using Polynomials, PyPlot, LinearAlgebra
# force IJulia to display as LaTeX rather than HTML
Base.showable(::MIME"text/html", ::Polynomial) = false

# some "random" matrix:
A = [ 0.325  -0.075   0.075  -0.075
      0.025   0.225  -0.025  -0.275
       0.15   -0.05    0.25   -0.05 
       -0.1    -0.1     0.1     0.4  ]

λ = eigvals(A)

ξ = range(0,0.6,length=100)
plot(ξ, [det(A - λ*I) for λ in ξ], "r-")
plot(ξ, 0ξ, "k--")
plot(λ, 0λ, "bo")
xlabel(L"\lambda")
ylabel(L"\det(A - \lambda I)")
title("characteristic polynomial")

A = [ 1 1
     -2 4 ]

eigvals(A)

A * [1, 1]

A * [1, 2]

λ, X = eigen(A)

λ

X

X[:,1] / X[1,1]  # first column, with first entry scaled to 1

X[:,2] / X[1,2] # second column, with second entry scaled to 1

A^5 * [1,1] # gives 2⁵ * [1,1] = [32, 32]

y = A^100.0 * [2,3]

2^100.0 * [1,1] + 3^100.0 * [1,2]  # same, but computed with eigenvectors

y[2] / y[1]

z = A^100.0 * [17,-5]

z[2]/z[1] # approximately 2 since z is nearly parallel to (1,2)

eigvals(A')

Y = eigvecs(A')
Y ./ Y[1,:]' # normalize so that the first components are 1, for easier comparison

eigvals([1 3
        -2 4])

(5 + sqrt(15)*im) / 2

w = prod([Polynomial([-n, 1.0]) for n = 1:10])

roots(w)

N = 10
w = prod([Polynomial([-n, 1.0]) for n = 1:N])
fig = figure()
#using Interact
#@manipulate for logR in -9:0.1:-1
for logR in -9:0.5:-1
    display(
    withfig(fig) do
        plot(1:N, zeros(10), "r*")
        R = exp10(logR)
        for i = 1:100
            r = roots(Polynomial(coeffs(w) .* (1 .+ R .* randn(N+1))))
            plot(real(r), imag(r), "b.")
        end
        xlabel("real part of roots")
        ylabel("imaginary part of roots")
        title("roots of \$(x-1)\\cdots(x-10)\$ with coeffs perturbed by R=$R")
    end
    )
end

function companion(p::Polynomial)
    c = coeffs(p)
    n = degree(p)
    c = c[1:n] / c[end]
    C = [ [ zeros(n-1)'; I ] -c ]'
    return C
end

p = Polynomial([-2, 1]) * Polynomial([-3, 1]) # (x - 2) * (x - 3)

C = companion(p)

eigvals(C)

 # (x - 2) * (x - 3) * (x - 4) * (x + 1)
p = Polynomial([-2, 1]) * Polynomial([-3, 1]) * Polynomial([-4, 1]) * Polynomial([1, 1])

C = companion(p)

eigvals(C)

@which roots(p)

A1000 = rand(1000,1000)
LinearAlgebra.BLAS.set_num_threads(1) # use 1 cpu for benchmarking
@time lu(A1000)
@time qr(A1000)
@time eigvals(A1000)
@time eigen(A1000);

@elapsed(eigen(A1000)) / @elapsed(lu(A1000))