Roots of a cubic polynomial

From Rosetta Code
Roots of a cubic polynomial is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Find all eigenvalues of a real 3x3 matrix and estimate their errors.

See Wikipedia Cubic equation. [1] Example.

a*x**3 + b*x**2 + c*x + d = 0

matrix: [[1,-1,0], [0,1,-1],[0,0,1]]

polynomial: [a,b,c,d] =[1, -3,+3, -1]

roots: [1,1,1]


J

Using J's polynomial primitive for most of the heavy lifting here. Also using J's generalized determinant to find the characteristic polynomial (first turning each element of the matrix into a polynomial then finding the determinant of the result). See also Oh, No, Not Eigenvalues Again!.

Implementation:

NB. matrix characteristic polynomial
mcp=: {{ ;pdf/ .ppr y,each -=i.#y }}
pdf=: -/@,:L:0      NB. polynomial difference
ppr=: +//.@(*/)L:0  NB. polynomial product
proots=: 1 {:: p.   NB. polynomial roots

task=: {{
  poly=. mcp y
  roots=. proots poly 
  errors=. poly p."1 0 roots
  lines=. (,:'matrix: '),&.|:&": y
  lp=. 'coefficients: ',":poly
  lr=. 'roots:  ',":roots
  le=. 'errors: ',":errors
  lines,lp,lr,:le
}}

Task examples:

   task 1 _1,0 1 1,:0 0 1
matrix: 1 _1 0         
        0  1 1         
        0  0 1         
coefficients: 1 _3 3 _1
roots:  1 1 1          
errors: 0 0 0          
   
   task _2 _4 2,_2 1 2,:4 2 5
matrix: _2 _4 2          
        _2  1 2          
         4  2 5          
coefficients: _90 27 4 _1
roots:  6 _5 3           
errors: 0 0 0            
   
   task 1 _1,0 1 _1,:0 0 1
matrix: 1 _1  0        
        0  1 _1        
        0  0  1        
coefficients: 1 _3 3 _1
roots:  1 1 1          
errors: 0 0 0          
   
   task 2 0,0 _1,:0 0 _1
matrix: 2  0  0       
        0 _1  0       
        0  0 _1       
coefficients: 2 3 0 _1
roots:  2 _1 _1       
errors: 0 0 0         
   
   task 2 0,0 3 4,:0 4 9
matrix: 2 0 0             
        0 3 4             
        0 4 9             
coefficients: 22 _35 14 _1
roots:  11 2 1            
errors: 0 0 0             
   
   task 1 0,0,.+.1ad_45 1ad45
matrix: 1        0         0                              
        0 0.707107 _0.707107                              
        0 0.707107  0.707107                              
coefficients: 1 _2.41421 2.41421 _1                       
roots:  1 0.707107j0.707107 0.707107j_0.707107            
errors: 0 4.44089e_15j_1.77636e_15 4.44089e_15j1.77636e_15

We could instead use lapack. The lapack approach would be faster and more stable for very large matrices, but that's not an issue here. Also, the approach we use here better fits the current title of this task.

jq

rootPoly is adapted from Wren

The first few functions are taken from Polynomial_long_division#jq but are reproduced here for clarity.

In the following, a method for computing the roots of a polynomial of arbitrary degree is used, and so no restriction to polynomials of degree 3 is imposed.

Note that the polynomial of degree n, P = SIGMA (c[i] * x^i), is represented by the array [ c[0] ... c[n] ], and the estimated error in a root, r, is computed using the formula (dY / (dP/dX)(r)) except when the denominator or the expression as a whole is computed as 0, in which case a very small value (almost indistinguishable from 0.0) is used.

# Input: a possibly empty numeric array [c0, ... ] 
# Emit the canonical form of the polynomial: SIGMA c[i] * x^i
def canonical:
  if length == 0 then [0]
  elif .[-1] == 0 then .[:-1]|canonical
  else .
  end;

# string representation
def poly2s: "Polynomial(\(join(",")))";

# Polynomial division
# Output [ quotient, remainder]
def divrem($divisor):
   ($divisor|canonical) as $divisor
   | { curr: canonical}
   | .base = ((.curr|length) - ($divisor|length))
   | until( .base < 0;
          (.curr[-1] / $divisor[-1]) as $res
          | .result += [$res]
          | .curr |= .[0:-1]
          |  reduce range (0;$divisor|length-1) as $i (.;
                .curr[.base + $i] +=  (- $res * $divisor[$i])  )
          | .base += -1
        )
   | [(.result | reverse),  (.curr | canonical)] ;

# Evaluate a polynomial at $x
# Input: the array representation of a polynomial
def evalPoly($x):
  length as $c
  | . as $coeffs
  | reduce range(1;length+1) as $i (0; . * $x + $coeffs[$c - $i]) ;

# Differentiate a polynomial
# Input: the array representation of the polynomial to be differentiated
def diffPoly:
  (length - 1) as $c
  | if $c == 0 then [0]
    else . as $coefs
    | reduce range(0; $c) as $i ([]; .[$i] = ($i+1) * $coefs[$i+1])
    end;

# Emit a stream
# No check is made that the input represents a quadratic polynomial
def quadraticRoots:
  . as $coefs
  | ($coefs[1]*$coefs[1] - 4*$coefs[0]*$coefs[2]) as $d
  | select($d >= 0)
  | ($d|sqrt) as $s
  | -($s + $coefs[1]), ($s - $coefs[1])
  | . / (2 * $coefs[2]);

# Attempt to find a real root of a polynomial using Newton's method from
# an initial guess, a given tolerance, a maximum number of iterations
# and a given multiplicity (usually 1).
# If a root is found, it is returned, otherwise 'null' is returned.
# If the root is near an integer, check if the integer is in fact the root.
# If the polynomical degree is 0, return null.
# Input: [c0, c1 ... ]
def rootPoly($guess; $tol; $maxIter; $mult):
  . as $coefs
  | (length - 1) as $deg
  | if $deg == 0 then null
    elif $deg == 1 then -$coefs[0]/$coefs[1]
    elif $deg == 2 then first(quadraticRoots)
    elif evalPoly(0) == 0 then 0
    else 0.001 as $eps
    | diffPoly as $deriv
    | { x0: $guess, iter: 1, return: null }
    | until(.return;
        .x0 as $x0
        | ($deriv|evalPoly($x0)) as $den
        | if $den == 0
          then .x0 |= (if . >= 0 then . + $eps else . - $eps end)
          else ($coefs|evalPoly($x0)) as $num
          | (.x0 - ($num/$den) * $mult) as $x1
          | if (($x1 - .x0)|length) <= $tol  # abs
            then ($x1 | round) as $r
            | if (($r - $x1)|length) <= $eps and ($coefs|evalPoly($r)) == 0
              then .return = $r
              else .return = $x1
              end
            else .x0 = $x1
            end
          | if .iter == $maxIter then .return = true
            else .iter += 1
            end
          end)
    | if .return != true then .return
      else (.x0|round) as $x0
      | if ($coefs | evalPoly($x0)) == 0 then $x0
        else null
        end
      end
    end;

# Convenience versions of rootPoly/4:
def rootPoly($guess):
  rootPoly($guess; 1e-15; 100; 1);
def rootPoly:
  rootPoly(0.001; 1e-15; 100; 1);

# Emit a stream of real roots
def roots:
  if length==3 then quadraticRoots
  else rootPoly as $root
  | select($root)
  | $root,
    # divide by (x - $root)
    ( divrem( [- $root, 1] ) as [$div, $rem]
      | $div
      | roots )
  end ;

# Given $root is an estimated root of the input polynomial,
# estimate the error dx from deriv = (dY / $root)
# except that if $root or deiv is 0, then use a very small value:
# since 1E-324 == 1E-325 but 1E-323 != 1E-324, we choose 1E-323
def estimatedDeltaX($root):
  1E-323 as $tiny
  | evalPoly($root) as $dY
  | (diffPoly | evalPoly($root)) as $deriv
  | if $deriv == 0 then $tiny
    else (($dY / $deriv) | length) as $dx
    | if $dx == 0 then $tiny
      else $dx
      end
    end ;

def roots_with_errors:
  rootPoly as $root
  | select($root)
  | estimatedDeltaX($root) as $dx
  | [$root, $dx],
    # divide by (x - $root)
    ( divrem( [- $root, 1] ) as [$div, $rem]
      | $div
      | roots_with_errors );

def illustration($text; $poly):
  $text,
  ($poly | roots_with_errors),
  "";

"Each root and corresponding estimated error is presented in an array.",
illustration("X^3 = 1";    [-1, 0, 0, 1] ),
illustration("X^3 - 4X^2 - 27X + 90 = 0";  # [-5, 3 ,6]
  [90, -27, -4, 1] ),
illustration("(X - 1)^3";  [-1,3,-3,1] ),
illustration("X^2 = 2";    [-2,0,1] ),
illustration("X^3 = 2";    [-2,0,0,1] )
Output:
Each root and corresponding estimated error is presented in an array.
X^3 = 1
[1,1e-323]

X^3 - 4X^2 - 27X + 90 = 0
[3,1e-323]
[-5,1e-323]
[6,1e-323]

(X - 1) ^3
[1,1e-323]
[1,1e-323]
[1,1e-323]

X^2 = 2
[-1.4142135623730951,1.570092458683775e-16]
[1.4142135623730951,1e-323]

X^3 = 2
[1.2599210498948732,1e-323]


Julia

Julia has extensive matrix math support via the LAPACK library. The Horner function is from the Raku example.

using LinearAlgebra

""" Cubic polynomial using Horner's method. """
horner(xarr, a, b, c, d) = @. ((a * xarr + b) * xarr + c) * xarr + d

""" Get parameters of characteristic cubic polynomial (with x^3 term 1) from 3x3 matrix """
polynomial3(t) = (1, -tr(t), (tr(t)^2 - tr(t^2)) / 2, -det(t)) # a, b, c, d

# test values
const sin45 = 1 / sqrt(2)
const testmats = [
    [1 -1  0; 0 1 1; 0 0 1],
    [-2 -4 2; -2 1 2; 4 2 5],
    [1 -1 0; 0 1 -1; 0 0 1],
    [2 0 0; 0 -1 0; 0 0 -1],
    [2 0 0; 0 3 4; 0 4 9],
    [1 0 0; 0 sin45 -sin45; 0 sin45 sin45],
]

for mat in testmats
    display(mat)
    a, b, c, d = polynomial3(mat)
    print("The characteristic polynomial is: ")
    println("x³ ", b < 0 ? "- $(-b)" : " + $b", "x^² ", c < 0 ? "- $(-c)" : "+ $c", "x ",
       d < 0 ? "- $(-d)" : "+ $d")
    evals = eigvals(mat)
    println("The LAPACK library computed eigenvalues are: $evals")
    println("Errors are: ", horner(evals, a, b, c, d), "\n")
end
Output:
3×3 Matrix{Float64}:
 1.0  -1.0  0.0
 0.0   1.0  1.0
 0.0   0.0  1.0
The characteristic polynomial is: x³ - 3.0x^² + 3.0x - 1.0
The LAPACK library computed eigenvalues are: [1.0, 1.0, 1.0]
Errors are: [0.0, 0.0, 0.0]

3×3 Matrix{Float64}:
 -2.0  -4.0  2.0
 -2.0   1.0  2.0
  4.0   2.0  5.0
The characteristic polynomial is: x³ - 4.0x^² - 27.0x + 90.0
The LAPACK library computed eigenvalues are: [-5.000000000000005, 3.0, 6.0]
Errors are: [-4.831690603168681e-13, 0.0, 0.0]

3×3 Matrix{Float64}:
 1.0  -1.0   0.0
 0.0   1.0  -1.0
 0.0   0.0   1.0
The characteristic polynomial is: x³ - 3.0x^² + 3.0x - 1.0
The LAPACK library computed eigenvalues are: [1.0, 1.0, 1.0]
Errors are: [0.0, 0.0, 0.0]

3×3 Matrix{Float64}:
 2.0   0.0   0.0
 0.0  -1.0   0.0
 0.0   0.0  -1.0
The characteristic polynomial is: x³  + -0.0x^² - 3.0x - 2.0
The LAPACK library computed eigenvalues are: [-1.0, -1.0, 2.0]
Errors are: [0.0, 0.0, 0.0]

3×3 Matrix{Float64}:
 2.0  0.0  0.0
 0.0  3.0  4.0
 0.0  4.0  9.0
The characteristic polynomial is: x³ - 14.0x^² + 35.0x - 22.0
The LAPACK library computed eigenvalues are: [1.0, 2.0, 11.0]
Errors are: [0.0, 0.0, 0.0]

3×3 Matrix{Float64}:
 1.0  0.0        0.0
 0.0  0.707107  -0.707107
 0.0  0.707107   0.707107
The characteristic polynomial is: x³ - 2.414213562373095x^² + 2.414213562373095x - 0.9999999999999998
The LAPACK library computed eigenvalues are: ComplexF64[0.7071067811865475 - 0.7071067811865475im, 0.7071067811865475 + 0.7071067811865475im, 1.0 + 0.0im]
Errors are: ComplexF64[1.1102230246251565e-16 + 1.1102230246251565e-16im, 1.1102230246251565e-16 - 1.1102230246251565e-16im, 2.220446049250313e-16 + 0.0im]

Phix

Translation of: Wren

Note the six functions matrix_trace() ... root_poly() are straightforward translations of code from matrix.wren and math.wren, and are not tested beyond matching the output here. The previously "suspect use of reverse() in eigenvalues(), which kinda brute-forces it to work" is now in Wren too.

with javascript_semantics
function matrix_trace(sequence m)
    atom res = 0
    for i=1 to length(m) do
        res += m[i,i]
    end for
    return res
end function

function matrix_minor(sequence m, integer x, y)
    integer l = length(m)-1
    sequence res = repeat(repeat(0,l),l)
    for i=1 to l do
        integer ri = i+(i>=x)
        for j=1 to l do
            integer rj = j+(j>=y)
            res[i,j] = m[ri,rj]
        end for
    end for
    return res
end function

function matrix_determinant(sequence m)
    // Returns the determinant of the current instance 
    // provieded it's square, using Laplace expansion.
    integer l = length(m)
    if l=1 then return m[1,1] end if
    if l=2 then return m[2,2]*m[1,1]-m[1,2]*m[2,1] end if
    integer sgn = 1
    atom det = 0
    for i=1 to l do
        sequence mm = matrix_minor(m,1,i)
        det += sgn*m[1,i]*matrix_determinant(mm)
        sgn = -sgn
    end for
    return det
end function

function cofactors(sequence m)
    integer l = length(m)
    sequence res = repeat(repeat(0,l),l)
    for i=1 to l do
        for j=1 to l do
            sequence mm = matrix_minor(m,i,j)
            atom d = matrix_determinant(mm)
            res[i,j] = d*power(-1,i+j)
        end for
    end for
    return res
end function

function diff_poly(sequence coefs)
    // Returns the coefficients of a polynomial following differentiation.
    // Polynomials are represented as described in 'eval_poly' below.
    integer c = length(coefs)-1
    if c=0 then return {0} end if
    sequence deriv = repeat(0,c)
    for i=1 to c do deriv[i] = (c-i+1) * coefs[i] end for
    return deriv
end function

function root_poly(sequence coefs, atom guess=0.001, tol=1e-15, maxIter=100, mult=1)
    integer deg = length(coefs)-1
    if deg=0 then return null end if
    if deg=1 then return -coefs[2]/coefs[1] end if
    if deg=2 and coefs[2]*coefs[2] - 4*coefs[1]*coefs[3] < 0 then return null end if
    if eval_poly(coefs, 0)=0 then return 0 end if
    sequence deriv = diff_poly(coefs)
    atom eps = 0.001,
          x0 = guess
    for iter=1 to maxIter do
        atom den = eval_poly(deriv, x0)
        if den=0 then
            x0 = iff(x0>=0 ? x0 + eps : x0 - eps)
        else
            atom num = eval_poly(coefs, x0),
                  x1 = x0 - num/den * mult
            if abs(x1-x0)<=tol then
                atom r = round(x1)
                if abs(r-x1)<=eps and eval_poly(coefs, r)=0 then return r end if
                return x1
            end if
            x0 = x1
        end if
    end for
    x0 = round(x0)
    if eval_poly(coefs, x0)=0 then return x0 end if
    return null
end function

function poly_div(sequence num,den)
    sequence curr = trim_tail(num,0),
            right = trim_tail(den,0),
              res = {}
    integer base = length(curr)-length(right)
    while base>=0 do
        atom r = curr[-1] / right[-1]
        res &= r
        curr = curr[1..-2]
        for i=1 to length(right)-1 do
            curr[base+i] -= r * right[i]
        end for
        base -= 1
    end while
    return res
end function

include complex.e

// assumes real coefficients, highest to lowest order
function quadratic(atom a, b, c)
    atom d = b*b - 4*a*c
    if d=0 then                 // single real root
        atom root = -b/(a*2)
        return {root, root}
    elsif d>0 then              // two real roots
        atom sr = sqrt(d)
        d = iff(b<0 ? sr - b : -sr - b)
        return {d/(2*a), 2*c/d}
    end if
    // two complex roots
    atom den = 1/(2*a)
    complex t1 = complex_new(-b*den, 0),
            t2 = complex_new(0, sqrt(-d)*den)
    return {complex_add(t1,t2), complex_sub(t1,t2)}
end function

// same assumption as quadratic but 'x' can be real or complex.
function eval_poly(sequence coefs, object x)
    object res = coefs[1]
    for i=2 to length(coefs) do
        if complex(x) then
            res = complex_add(complex_mul(res,x),coefs[i])
        else
            res = res*x + coefs[i]
        end if
    end for
    return res
end function

function poly(sequence si)
    -- display helper
    string r = ""
    for t=length(si) to 1 by -1 do
        atom sit = si[t]
        if sit!=0 then
            if sit=1 and t>1 then
                r &= iff(r=""? "":" + ")
            elsif sit=-1 and t>1 then
                r &= iff(r=""?"-":" - ")
            else
                if r!="" then
                    r &= iff(sit<0?" - ":" + ")
                    sit = abs(sit)
                end if
                r &= sprintf("%g",sit)
            end if
            r &= iff(t>1?"x"&iff(t>2?sprintf("^%d",t-1):""):"")
        end if
    end for
    if r="" then r="0" end if
    return r
end function

function eigenvalues(sequence m)
    // find the characteristic polynomial
    sequence cp = {1,
                  -matrix_trace(m),
                   matrix_trace(cofactors(m)),
                  -matrix_determinant(m)},
          roots = {},
           errs = {}
    // find first root
    roots &= root_poly(cp)
    errs &= eval_poly(cp, roots[1])
    // divide out to get quadratic
    sequence den = {-roots[1],1},
             qdr = poly_div(reverse(cp),den)
    // find second and third roots
    roots &= quadratic(qdr[1],qdr[2],qdr[3])
    errs = append(errs,eval_poly(cp, roots[2]))
    errs = append(errs,eval_poly(cp, roots[3]))
    string pm = ppf(m,{pp_Nest,1,pp_Indent,8,pp_IntFmt,"%2d"})
    return {pm, poly(reverse(cp)), roots, errs}
end function

constant v = 1/sqrt(2)
constant tests = {{{ 1, -1,  0},
                   { 0,  1,  1},
                   { 0,  0,  1}},
                  {{-2, -4,  2},
                   {-2,  1,  2},
                   { 4,  2,  5}},
                  {{ 1, -1,  0},
                   { 0,  1, -1},
                   { 0,  0,  1}},
                  {{ 2,  0,  0},
                   { 0, -1,  0},
                   { 0,  0, -1}},
                  {{ 2,  0,  0},
                   { 0,  3,  4},
                   { 0,  4,  9}},
                  {{ 1,  0,  0},
                   { 0,  v, -v},
                   { 0,  v,  v}}}
constant fmt = """
Matrix: %s
characteristic polynomial: %s
eigenvalues: %v
errors: %v

"""

for m in tests do
    printf(1,fmt,eigenvalues(m))
end for
Output:
Matrix: {{ 1,-1, 0},
         { 0, 1, 1},
         { 0, 0, 1}}

characteristic polynomial: x^3 - 3x^2 + 3x - 1
eigenvalues: {1,1,1}
errors: {0,0,0}

Matrix: {{-2,-4, 2},
         {-2, 1, 2},
         { 4, 2, 5}}

characteristic polynomial: x^3 - 4x^2 - 27x + 90
eigenvalues: {3,6,-5}
errors: {0,0,0}

Matrix: {{ 1,-1, 0},
         { 0, 1,-1},
         { 0, 0, 1}}

characteristic polynomial: x^3 - 3x^2 + 3x - 1
eigenvalues: {1,1,1}
errors: {0,0,0}

Matrix: {{ 2, 0, 0},
         { 0,-1, 0},
         { 0, 0,-1}}

characteristic polynomial: x^3 - 3x - 2
eigenvalues: {-1,2,-1}
errors: {0,0,0}

Matrix: {{ 2, 0, 0},
         { 0, 3, 4},
         { 0, 4, 9}}

characteristic polynomial: x^3 - 14x^2 + 35x - 22
eigenvalues: {1,11,2}
errors: {0,0,0}

Matrix: {{ 1, 0, 0},
         { 0,0.7071067812,-0.7071067812},
         { 0,0.7071067812,0.7071067812}}

characteristic polynomial: x^3 - 2.41421x^2 + 2.41421x - 1
eigenvalues: {1.0,{0.7071067812,0.7071067812},{0.7071067812,-0.7071067812}}
errors: {2.22044605e-16,{2.22044605e-16,-2.775557561e-16},{2.22044605e-16,2.775557561e-16}}

Translation of: Raku

Note "wren2" gives eigenvalues of {6,3,-5}...

with javascript_semantics

function polynomial(sequence t)
    atom a = 1,                                 // create characteristic polynomial
         b = -(t[1,1] + t[2,2] + t[3,3]),                   // = -trace
         c = ( t[1,1]*t[2,2] + t[2,2]*t[3,3] + t[3,3]*t[1,1] )
             -(t[2,3]*t[3,2] + t[3,1]*t[1,3] + t[1,2]*t[2,1] ),
         d = - t[1,1] * t[2,2] * t[3,3]
             - t[1,2] * t[2,3] * t[3,1]
             - t[1,3] * t[2,1] * t[3,2]
             + t[1,1] * t[3,2] * t[2,3]
             + t[1,2] * t[2,1] * t[3,3]
             + t[1,3] * t[2,2] * t[3,1];                    // = -determinant
    return {a, b, c, d}
end function

function reduction(sequence poly)
    atom {a, b, c, d} = poly,
         b3 = power(b,3),
         b2 = power(b,2),
         a3 = power(a,3),
         a2 = power(a,2),
         delta = 18*a*b*c*d - 4*b3*d + b2*power(c,2) - 4*a*power(c,3) - 27*a2*power(d,2),
         p = (3*a*c - b2) / (3*a2),
         q = (2*b3 - 9*a*b*c + 27*a2*d) / (27*a3),
         d0 = b2 - 3*a*c,
         d1 = 2*b3 - 9*a*b*c + 27*a2*d;
    return { delta, p, q, d0, d1 }
end function

include complex.e

function horner(object x, atom a, b, c, d ) // cubic polynomial using horner's method
    if complex(x) then
        complexn r = 0
        for f in {a,b,c,d} do
            r = complex_add(complex_mul(r,x),f)
        end for
        return r
    end if
    return ((a*x + b)*x + c)*x + d
end function

function solution(atom a, b, /*c*/, /*d*/, delta, p, q, d0, d1)
    sequence x
    if abs(delta)<1e-12 then        // multiple real roots
        if abs(p)<1e-12 then        // triple equal real roots
            x = {0,0,0}
        else                        // double real root
            atom x23 = -3*q/(2*p)
            x = {3*q/p, x23, x23}
        end if
    elsif delta > 0 then            // three distinct real roots
        x = {}
        for i=1 to 3 do
            x &= 2*sqrt(-p/3) * cos(arccos(sqrt(-3/p)*3*q/(2*p))/3 - 2*PI*(i-1)/3)
        end for
    else // delta < 0               // one real root and two complex conjugate roots
        complexn g
        if d0=0 and d1<0 then
            g = -power(-d1,1/3)
        elsif d0<0 and d1=0 then
            g = -sqrt(-d0)
        else
            complex s = complex_sqrt(power(d1,2) - 4*power(d0,3)),
                   g1 = complex_power(complex_div(complex_add(d1,s),2),1/3),
                   g2 = complex_power(complex_div(complex_sub(d1,s),2),1/3)
            g = iff(g1={0,0} ? g2 : g1)
        end if
        complex z = complex_mul(g,complex_div(complex_add(-1,complex_sqrt(-3)),2)),
                x1 = complex_mul(-1/3,complex_add(g,complex_div(d0,g))),
                x2 = complex_mul(-1/3,complex_add(z,complex_div(d0,z)))
        x = {x1,x2,complex_conjugate(x2)}
        for i=1 to 3 do
            x[i] = complex_sub(x[i],(1/3)*b/a)
        end for
        return x
    end if
    for i=1 to 3 do
        x[i] -= (1/3)*b/a
    end for
    return x
end function

function spectrum(sequence m)
    sequence poly = polynomial(m)
    atom {a, b, c, d} = poly,
         {delta, p, q, d0, d1} = reduction(poly);
    sequence s = solution(a, b, c, d, delta, p, q, d0, d1 ),
             e = apply(true,horner,{s,a,b,c,d})
    return {s, e}
end function

constant r = 1/sqrt(2),
         m = {{"wren1",    {{ 1, -1,  0}, { 0,  1,  1}, { 0,  0,  1}}},
              {"wren2",    {{-2, -4,  2}, {-2,  1,  2}, { 4,  2,  5}}},
              {"triple",   {{ 1, -1,  0}, { 0,  1, -1}, { 0,  0,  1}}},
              {"double",   {{ 2,  0,  0}, { 0, -1,  0}, { 0,  0, -1}}},
              {"distinct", {{ 2,  0,  0}, { 0,  3,  4}, { 0,  4,  9}}},
              {"rotation", {{ 1,  0,  0}, { 0,  r, -r}, { 0,  r,  r}}}}
for i,pair in m do
    {string desc, sequence t} = pair
    printf(1," %d:  %8s matrix: %v\n",{i,desc,t})
    sequence poly = polynomial(t),
           reduct = reduction(poly),
            {s,e} = spectrum(t)
    printf(1,"          polynomial: %v\n",{poly})
    printf(1,"           reduction: %v\n",{reduct})
    printf(1,"         eigenvalues: %v\n",{s})
    printf(1,"              errors: %v\n",{e})
end for
Output:
 1:     wren1 matrix: {{1,-1,0},{0,1,1},{0,0,1}}
          polynomial: {1,-3,3,-1}
           reduction: {0,0,0,0,0}
         eigenvalues: {1,1,1}
              errors: {0,0,0}
 2:     wren2 matrix: {{-2,-4,2},{-2,1,2},{4,2,5}}
          polynomial: {1,-4,-27,90}
           reduction: {69696,-32.33333333,49.25925926,97,1330}
         eigenvalues: {6,3.0,-5.0}
              errors: {0,-4.263256414e-14,-1.705302566e-13}
 3:    triple matrix: {{1,-1,0},{0,1,-1},{0,0,1}}
          polynomial: {1,-3,3,-1}
           reduction: {0,0,0,0,0}
         eigenvalues: {1,1,1}
              errors: {0,0,0}
 4:    double matrix: {{2,0,0},{0,-1,0},{0,0,-1}}
          polynomial: {1,0,-3,-2}
           reduction: {0,-3,-2,9,-54}
         eigenvalues: {2,-1,-1}
              errors: {0,0,0}
 5:  distinct matrix: {{2,0,0},{0,3,4},{0,4,9}}
          polynomial: {1,-14,35,-22}
           reduction: {8100,-30.33333333,-61.92592593,91,-1672}
         eigenvalues: {11,2.0,1.0}
              errors: {0,1.065814103e-14,-1.77635684e-14}
 6:  rotation matrix: {{1,0,0},{0,0.7071067812,-0.7071067812},{0,0.7071067812,0.7071067812}}
          polynomial: {1,-2.414213562,2.414213562,-1.0}
           reduction: {-0.686291501,0.4714045208,-0.0994922778,-1.414213562,-2.686291501}
         eigenvalues: {{1.0,0},{0.7071067812,-0.7071067812},{0.7071067812,0.7071067812}}
              errors: {{-2.22044605e-16,0},{-3.330669073e-16,-3.885780587e-16},{-3.330669073e-16,3.885780587e-16}}

Complex results on the last entry are just shown as {real,imag} pairs.

Raku

my \r = 1/sqrt(2);
my Pair @m = [[
      triple       => [[ 1, -1,  0], [ 0,  1, -1], [ 0,  0,  1] ],
      double       => [[ 2,  0,  0], [ 0, -1,  0], [ 0,  0, -1] ],
      distinct     => [[ 2,  0,  0], [ 0,  3,  4], [ 0,  4,  9] ],
      rotation     => [[ 1,  0,  0], [ 0,  r, -r], [ 0,  r,  r] ],
   ]];
for ^@m {
   my Pair $pair = @m[$_];                  print "\n  {1+$_}:  {$pair.key} ";
   my @t = $pair.value.[];                  say 'matrix: 	  ', @t.raku;
   my @poly = polynomial( @t );             say 'polynomial:  ', @poly.raku;
   my @reduction = reduction( |@poly );     say 'reduction :  ', @reduction.raku;
   my ( $s, $e ) = spectrum( @t );          say 'eigenvalues: ', @$s.raku;
                                            say 'errors: 	  ', @$e.raku;
}
sub horner( \x, \a, \b, \c, \d ) {          # cubic polynomial using horner's method
    ((a * x + b) * x + c) * x + d
}
sub polynomial( @t ) {
    my \a = 1;                                  # create characteristic polynomial
    my \b = -(@t[0;0] + @t[1;1] + @t[2;2]);                     # = -trace
    my \c = ( @t[0;0]*@t[1;1] + @t[1;1]*@t[2;2] + @t[2;2]*@t[0;0] )
            -(@t[1;2]*@t[2;1] + @t[2;0]*@t[0;2] + @t[0;1]*@t[1;0] );
    my \d = - @t[0;0] * @t[1;1] * @t[2;2]
            - @t[0;1] * @t[1;2] * @t[2;0]
            - @t[0;2] * @t[1;0] * @t[2;1]
            + @t[0;0] * @t[2;1] * @t[1;2]
            + @t[0;1] * @t[1;0] * @t[2;2]
            + @t[0;2] * @t[1;1] * @t[2;0];                  # = -determinant
    return [ a, b, c, d]>>.narrow;
}
sub reduction( \a, \b, \c, \d ) {
    my \delta = 18 * a * b * c * d - 4 * b **3 * d + b**2 * c**2 - 4 * a * c**3 - 27 * a**2 * d**2;
    my \p = (3 * a * c - b * b) / (3 * a * a);
    my \q = (2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d) / (27 * a ** 3);
    my \d0 = b*b - 3 * a * c;
    my \d1 = 2*b**3 - 9*a*b*c + 27 * a**2 * d;
    return [ delta, p, q, d0, d1 ];
}
sub solution( \a, \b, \c, \d, \delta, \p, \q, \d0, \d1 ) {
    my @x;
    if abs(delta) =~= 0 {           # say " multiple real roots ", p.raku ;
        if abs(p) =~= 0 {           # say " triple equal real roots: ";;
            @x[$_] = 0 for ^3;
        }
        else {                      # say ' double real root:';
            @x[0] = 3 * q / p;
            @x[1] = -3 * q /(2 * p);
            @x[2] = @x[1];
        }
    }
    elsif delta > 0 {               # say " three distinct real roots:";
        @x[$_] = 2*sqrt(-p/3) * cos( acos( sqrt( -3 / p ) * 3 * q /( 2 * p ) )/3 - 2*pi * $_/ 3 ) for ^3;
    }
    else {      # delta < 0;        say " one real root and two complex conjugate roots:";
        my $g = do {
            if d0 == 0 and d1 < 0 {
                -(-d1)**;
            }
            elsif d0 < 0 and d1 == 0 {
                -sqrt(-d0);
            }
            else {
                my \s = Complex( d1**2 - 4 * d0**3 ).sqrt;
                my \g1 = (( d1 + s )/2)**;
                my \g2 = (( d1 - s )/2)**;
                g1 == 0 ?? g2 !! g1
            }
        }
        my Complex $z = $g * ( -1 + Complex(-3).sqrt )/2;
        @x[0] = - * ( $g + d0 / $g );
        @x[1] = - * ( $z + d0 / $z );
        @x[2] = @x[1].conj;
    }
    @x[$_] -=  * b / a for ^3;
    return @x;
}
sub spectrum( @mat ) {
    my ( \a, \b, \c, \d )         = polynomial( @mat );
    my (\delta, \p, \q, \d0, \d1) = reduction( a, b, c, d );
    my @s = solution( a, b, c, d, delta, p, q, d0, d1 )>>.narrow;
    my @e = @s.map( { horner( $_, a, b, c, d ) })>>.narrow;
    return ( @s, @e );
}
Output:
 1:  triple matrix: [[1, -1, 0], [0, 1, -1], [0, 0, 1]]
polynomial:  [1, -3, 3, -1]
reduction :  [0, 0.0, 0.0, 0, 0]
eigenvalues: [1, 1, 1]
errors:      [0, 0, 0]

  2:  double matrix: [[2, 0, 0], [0, -1, 0], [0, 0, -1]]
polynomial:  [1, 0, -3, -2]
reduction :  [0, -3.0, -2.0, 9, -54]
eigenvalues: [2, -1, -1]
errors:      [0, 0, 0]

  3:  distinct matrix: [[2, 0, 0], [0, 3, 4], [0, 4, 9]]
polynomial:  [1, -14, 35, -22]
reduction :  [8100, <-91/3>, <-1672/27>, 91, -1672]
eigenvalues: [11, 2, 0.9999999999999991e0]
errors:      [0, 0, -1.0658141036401503e-14]

  4: rotation matrix: [[1, 0, 0], [0, 0.7071067811865475e0, -0.7071067811865475e0], [0, 0.7071067811865475e0, 0.7071067811865475e0]]
polynomial:  [1, -2.414213562373095e0, 2.414213562373095e0, -0.9999999999999998e0]
reduction :  [-0.6862915010152442e0, 0.4714045207910316e0, -0.0994922778153789e0, -1.414213562373095e0, -2.68629150101523e0]
eigenvalues: [0.9999999999999992e0, <0.7071067811865478-0.7071067811865472i>, <0.7071067811865478+0.7071067811865472i>]
errors:      [0, <-3.3306690738754696e-16-3.885780586188048e-16i>, <-3.3306690738754696e-16+3.885780586188048e-16i>]

RPL

Works with: HP version 49
« JORDAN 4 ROLL DROP NIP
» 'TASK' STO
[[1 -1 0][0 1 -1][0 0 1]] TASK
Output:
2: 'X^3-3*X^2+3*X-1'
1: [1 1 1]

If the characteristic polynomial is not needed, the EGVL function directly returns the eigenvalues as a vector.

Wren

Library: Wren-matrix
Library: Wren-math
Library: Wren-complex
Library: Wren-fmt

The eigenvalues of a 3 x 3 matrix will be the roots of its characteristic polynomial.

We borrow code from the Polynomial_long_division#Wren task to divide out this cubic polynomial after the first root is found and then code from the Roots_of_a_quadratic_function#Wren task to solve the resulting quadratic polynomial which may have complex roots. Note that the former code assumes that polynomials are presented starting from the lowest order term and all other functions used here from the highest order term. It is therefore necessary to reverse the order of terms when switching between the two.

import "./matrix" for Matrix
import "./math" for Math
import "./complex" for Complex
import "./fmt" for Fmt

class Polynom {
    // assumes factors start from lowest order term
    construct new(factors) {
        _factors = factors.toList
    }

    factors { _factors.toList }

    /(divisor) {
        var curr = canonical().factors
        var right = divisor.canonical().factors
        var result = []
        var base = curr.count - right.count
        while (base >= 0) {
            var res = curr[-1] / right[-1]
            result.add(res)
            curr = curr[0...-1]
            for (i in 0...right.count-1) {
                curr[base + i] = curr[base + i] - res * right[i]
            }
            base = base - 1
        }
        var quot = Polynom.new(result[-1..0])
        var rem = Polynom.new(curr).canonical()
        return [quot, rem]
    }

    canonical() {
        if (_factors[-1] != 0) return this
        var newLen = factors.count
        while (newLen > 0) {
            if (_factors[newLen-1] != 0) return Polynom.new(_factors[0...newLen])
            newLen = newLen - 1
        }
        return Polynom.new(_factors[0..0])
    }
}

// assumes real coefficients, highest to lowest order
var quadratic = Fn.new { |a, b, c|
    var d = b*b - 4*a*c
    if (d == 0) {
        // single real root
        var root = -b/(a*2)
        return [root, root]
    }
    if (d > 0) {
        // two real roots
        var sr = d.sqrt
        d = (b < 0) ? sr - b : -sr - b
        return [d/(2*a), 2*c/d]
    }
    // two complex roots
    var den = 1 / (2*a)
    var t1 = Complex.new(-b*den, 0)
    var t2 = Complex.new(0, (-d).sqrt * den)
    return [t1+t2, t1-t2]
}

// same assumption as quadratic but 'x' can be real or complex.
var evalPoly = Fn.new { |coefs, x|
    var c = coefs.count
    if (c == 1) return coefs[0]
    var sum = coefs[0]
    for (i in 1...c) sum = x * sum + coefs[i]
    return sum
}

var eigenvalues = Fn.new { |m|
    var roots = []
    var errs  = []
    // find the characteristic polynomial
    var cp = List.filled(4, 0)
    cp[0] = 1
    cp[1] = -m.trace
    cp[2] = m.cofactors.trace
    cp[3] = -m.det
    // find first root
    roots.add(Math.rootPoly(cp))
    errs.add(evalPoly.call(cp, roots[0]))
    // divide out to get quadratic
    var num = Polynom.new(cp[-1..0])
    var den = Polynom.new([-roots[0], 1])
    var qdr = (num/den)[0]
    // find second and third roots
    var coefs = qdr.factors[-1..0]
    roots = roots + quadratic.call(coefs[0], coefs[1], coefs[2])
    errs.add(evalPoly.call(cp, roots[1]))
    errs.add(evalPoly.call(cp, roots[2]))
    return [cp, roots, errs]
}

var v = 0.70710678118655
var arrs = [
    [
        [ 1, -1,  0],
        [ 0,  1,  1],
        [ 0,  0,  1]
    ],
    [
        [-2, -4,  2],
        [-2,  1,  2],
        [ 4,  2,  5]
    ],
    [
        [ 1, -1,  0],
        [ 0,  1, -1],
        [ 0,  0,  1]
    ],
    [
        [ 2,  0,  0],
        [ 0, -1,  0],
        [ 0,  0, -1]
    ],
    [
        [ 2,  0,  0],
        [ 0,  3,  4],
        [ 0,  4,  9]
    ],
    [
        [ 1,  0,  0],
        [ 0,  v, -v],
        [ 0,  v,  v]
    ]
]

for (i in 0...arrs.count) {
    var m = Matrix.new(arrs[i])
    System.print("For matrix:")
    if (i < arrs.count-1) Fmt.mprint(m, 2, 0) else Fmt.mprint(m, 17, 14)
    var eigs = eigenvalues.call(m)
    Fmt.print("\nwhose characteristic polynomial is:")
    Fmt.pprint("$n", eigs[0], "", "x")
    Fmt.print("\nIts eigenvalues are: $n", eigs[1])
    Fmt.print("and the corresponding errors are: $n\n", eigs[2])
}
Output:
For matrix:
| 1 -1  0|
| 0  1  1|
| 0  0  1|

whose characteristic polynomial is:
x³ - 3x² + 3x - 1

Its eigenvalues are: [1, 1, 1]
and the corresponding errors are: [0, 0, 0]

For matrix:
|-2 -4  2|
|-2  1  2|
| 4  2  5|

whose characteristic polynomial is:
x³ - 4x² - 27x + 90

Its eigenvalues are: [3, 6, -5]
and the corresponding errors are: [0, 0, 0]

For matrix:
| 1 -1  0|
| 0  1 -1|
| 0  0  1|

whose characteristic polynomial is:
x³ - 3x² + 3x - 1

Its eigenvalues are: [1, 1, 1]
and the corresponding errors are: [0, 0, 0]

For matrix:
| 2  0  0|
| 0 -1  0|
| 0  0 -1|

whose characteristic polynomial is:
x³ - 3x - 2

Its eigenvalues are: [-1, 2, -1]
and the corresponding errors are: [0, 0, 0]

For matrix:
| 2  0  0|
| 0  3  4|
| 0  4  9|

whose characteristic polynomial is:
x³ - 14x² + 35x - 22

Its eigenvalues are: [1, 11, 2]
and the corresponding errors are: [0, 0, 0]

For matrix:
| 1.00000000000000  0.00000000000000  0.00000000000000|
| 0.00000000000000  0.70710678118655 -0.70710678118655|
| 0.00000000000000  0.70710678118655  0.70710678118655|

whose characteristic polynomial is:
x³ - 2.4142135623731x² + 2.4142135623731x - 1

Its eigenvalues are: [1, 0.70710678118655 + 0.70710678118655i, 0.70710678118655 - 0.70710678118655i]
and the corresponding errors are: [0, 0 + 0i, 0 + 0i]