Eigenvalues of a Matrix Given a square matrix, output the matrix's eigenvalues. Each eigenvalue should be repeated a number of times equal to its algebraic multiplicity.

The eigenvalues of a matrix A are scalar values λ such that, for some column vector v, A*v = λ*v. They are also the solutions to the characteristic polynomial of A: det(A - λ*I) = 0 (where I is the identity matrix with the same dimensions as A).

Outputs must be accurate to 3 significant digits. All inputs and outputs will be within the representable range of numerical values for your chosen language.

Builtins are acceptable, but you are encouraged to include solutions that do not use builtins.

Test Cases

In these test cases, I represents the imaginary unit. Complex numbers are written in the form a + b*I. All outputs have 3 significant digits of precision.

[[42.0]] -&gt; [42.0]
[[1.0, 0.0], [0.0, 1.0]] -&gt; [1.00, 1.00]
[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]] -&gt; [16.1, -1.12, -1.24e-15]
[[1.2, 3.4, 5.6, 7.8], [6.3, 0.9, -5.4, -2.3], [-12.0, -9.7, 7.3, 5.9], [-2.5, 7.9, 5.3, 4.4]] -&gt; [7.20 + 5.54*I, 7.20 - 5.54*I, -4.35, 3.75]
[[-3.22 - 9.07*I, 0.193 + 9.11*I, 5.59 + 1.33*I, -3.0 - 6.51*I, -3.73 - 6.42*I], [8.49 - 3.46*I, -1.12 + 6.39*I, -8.25 - 0.455*I, 9.37 - 6.43*I, -6.82 + 8.34*I], [-5.26 + 8.07*I, -6.68 + 3.72*I, -3.21 - 5.63*I, 9.31 + 3.86*I, 4.11 - 8.82*I], [-1.24 + 9.04*I, 8.87 - 0.0352*I, 8.35 + 4.5*I, -9.62 - 2.21*I, 1.76 - 5.72*I], [7.0 - 4.79*I, 9.3 - 2.31*I, -2.41 - 7.3*I, -7.77 - 6.85*I, -9.32 + 2.71*I]] -&gt; [5.18 + 16.7*I, -24.9 - 2.01*I, -5.59 - 13.8*I, 0.0438 - 10.6*I, -1.26 + 1.82*I]
[[-30.6 - 73.3*I, 1.03 - 15.6*I, -83.4 + 72.5*I, 24.1 + 69.6*I, 52.3 + 2.68*I, 23.8 + 98.0*I, 96.8 + 49.7*I, -26.2 - 5.87*I, -52.4 + 98.2*I, 78.1 + 6.69*I], [-59.7 - 66.9*I, -26.3 + 65.0*I, 5.71 + 4.75*I, 91.9 + 82.5*I, -94.6 + 51.8*I, 61.7 + 82.3*I, 54.8 - 27.8*I, 45.7 + 59.2*I, -28.3 + 78.1*I, -59.9 - 54.5*I], [-36.0 + 22.9*I, -51.7 + 10.8*I, -46.6 - 88.0*I, -52.8 - 32.0*I, -75.7 - 23.4*I, 96.2 - 71.2*I, -15.3 - 32.7*I, 26.9 + 6.31*I, -59.2 + 25.8*I, -0.836 - 98.3*I], [-65.2 - 90.6*I, 65.6 - 24.1*I, 72.5 + 33.9*I, 1.47 - 93.8*I, -0.143 + 39.0*I, -3.71 - 30.1*I, 60.1 - 42.4*I, 55.6 + 5.65*I, 48.2 - 53.0*I, -3.9 - 33.0*I], [7.04 + 0.0326*I, -12.8 - 50.4*I, 70.1 - 30.3*I, 42.7 - 76.3*I, -3.24 - 64.1*I, 97.3 + 66.8*I, -11.0 + 16.5*I, -40.6 - 90.7*I, 71.5 - 26.2*I, 83.1 - 49.4*I], [-59.5 + 8.08*I, 74.6 + 29.1*I, -65.8 + 26.3*I, -76.7 - 83.2*I, 26.2 + 99.0*I, -54.8 + 33.3*I, 2.79 - 16.6*I, -85.2 - 3.64*I, 98.4 - 12.4*I, -27.6 - 62.3*I], [82.6 - 95.3*I, 55.8 - 73.6*I, -49.9 + 42.1*I, 53.4 + 16.5*I, 80.2 - 43.6*I, -43.3 - 3.9*I, -2.26 - 58.3*I, -19.9 + 98.1*I, 47.2 + 62.4*I, -63.3 - 54.0*I], [-88.7 + 57.7*I, 55.6 + 70.9*I, 84.1 - 52.8*I, 71.3 - 29.8*I, -3.74 - 19.6*I, 29.7 + 1.18*I, -70.6 - 10.5*I, 37.6 + 99.9*I, 87.0 + 19.0*I, -26.1 - 82.0*I], [69.5 - 47.1*I, 11.3 - 59.0*I, -84.3 - 35.1*I, -3.61 - 35.7*I, 88.0 + 88.1*I, -47.5 + 0.956*I, 14.1 + 89.8*I, 51.3 + 0.14*I, -78.5 - 66.5*I, 2.12 - 53.2*I], [0.599 - 71.2*I, 21.7 + 10.8*I, 19.9 - 97.1*I, 20.5 + 37.4*I, 24.7 + 40.6*I, -82.7 - 29.1*I, 77.9 + 12.5*I, 94.1 - 87.4*I, 78.6 - 89.6*I, 82.6 - 69.6*I]] -&gt; [262. - 180.*I, 179. + 117.*I, 10.3 + 214.*I, 102. - 145.*I, -36.5 + 97.7*I, -82.2 + 89.8*I, -241. - 104.*I, -119. - 26.0*I, -140. - 218.*I, -56.0 - 160.*I]
Haskell, 576 554 532 507 bytes



No built-ins!

import Data.Complex
s=sum
l=length
m=magnitude
i=fromIntegral
(&amp;)=zip
t=zipWith
(x!a)b=x*a+b
a#b=[[s$t(*)x y|y&lt;-foldr(t(:))([]&lt;$b)b]|x&lt;-a]
f a|let c=[1..l a];g(u,d)k|m&lt;-[t(+)a b|(a,b)&lt;-a#u&amp;[[s[d|x==y]|y&lt;-c]|x&lt;-c]]=(m,-s[s[b|(n,b)&lt;-c&amp;a,n==m]|(a,m)&lt;-a#m&amp;c]/i k)=snd&lt;$&gt;scanl g(0&lt;$c&lt;$c,1)c
p?x|let f=foldl1(x!);c=l p-1;n=i c;q p=init$t(*)p$i&lt;$&gt;[c,c-1..];o=f(q p)/f p;a|d&lt;-sqrt$(n-1)*(n*(o^2-f(q$q p)/f p)-o^2)=n/last(o-d:[o+d|m(o-d)&lt;m(o+d)])=last$p?(x-a):[x|m a&lt;1e-9]
z[a,b]=[-b/a]
z p=p?0:z(init$scanl1(p?0!)p)


Try it online!

Many thanks @ØrjanJohansen for a total of -47 bytes!

Explanation

First this computes the characteristic polynomial with the Faddeev–LeVerrier algorithm which is the function f. Then the function z computes all roots of that polynomial by iterating g which implements Laguerre's Method for finding a root, once a root is found it is removed and g gets called again until the polynomial has degree 1 which is trivially solved by z[a,b]=[-b/a].

Ungolfed

I re-inlined the functions sum,length,magnitude, fromIntegral, zipWith and (&amp;) as well as the little helper (!). The function faddeevLeVerrier corresponds to f, roots to z and g to laguerre respectively.

-- Transpose a matrix/list
transpose a = foldr (zipWith(:)) (replicate (length a) []) a

-- Straight forward implementation for matrix-matrix multiplication
(#) :: [[Complex Double]] -&gt; [[Complex Double]] -&gt; [[Complex Double]]
a # b = [[sum $ zipWith (*) x y | y &lt;- transpose b]|x&lt;-a]


-- Faddeev-LeVerrier algorithm
faddeevLeVerrier :: [[Complex Double]] -&gt; [Complex Double]
faddeevLeVerrier a = snd &lt;$&gt; scanl go (zero,1) [1..n]
  where n = length a
        zero = replicate n (replicate n 0)
        trace m = sum [sum [b|(n,b)&lt;-zip [1..n] a,n==m]|(m,a)&lt;-zip [1..n] m]
        diag d = [[sum[d|x==y]|y&lt;-[1..n]]|x&lt;-[1..n]]
        add as bs = [[x+y | (x,y) &lt;- zip a b] | (b,a) &lt;- zip as bs]
        go (u,d) k = (m, -trace (a#m) / fromIntegral k)
          where m = add (diag d) (a#u)


-- Compute roots by succesively removing newly computed roots
roots :: [Complex Double] -&gt; [Complex Double]
roots [a,b] = [-b/a]
roots   p   = root : roots (removeRoot p)
  where root = laguerre p 0
        removeRoot = init . scanl1 (\a b -&gt; root*a + b)

-- Compute a root of a polynomial p with an initial guess x
laguerre :: [Complex Double] -&gt; Complex Double -&gt; Complex Double
laguerre p x = if magnitude a &lt; 1e-9 then x else laguerre p new_x
  where evaluate = foldl1 (\a b -&gt; x*a+b)
        order' = length p - 1
        order  = fromIntegral $ length p - 1
        derivative p = init $ zipWith (*) p $ map fromIntegral [order',order'-1..]
        g  = evaluate (derivative p) / evaluate p
        h  = (g ** 2 - evaluate (derivative (derivative p)) / evaluate p)
        d  = sqrt $ (order-1) * (order*h - g**2)
        ga = g - d
        gb = g + d
        s = if magnitude ga &lt; magnitude gb then gb else ga
        a = order /s
        new_x = x - a