Roots of a cubic polynomial

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]

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

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

roots: [1,1,1]

Raku

Works with: Raku version 2023.11

Phix

Translation of: Wren

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 
    // provided 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 eval_poly(sequence coefs, atom x)
    atom res = coefs[1]
    for i=2 to length(coefs) do res = res*x + coefs[i] 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' above(!).
    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

function poly(sequence si)
    -- display helper
    string r = ""
    for t=length(si) to 1 by -1 do
        integer 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("%d",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 = repeat(0,3),
           errs = repeat(0,3)
    // find first root
    roots[1] = root_poly(cp)
    errs[1] = eval_poly(cp, roots[1])
    // divide out to get quadratic
    // aside: no idea why reverse helps...
    sequence num = reverse(cp),
             den = {-roots[1],1}
    num = reverse(poly_div(num,den))
    // find second root
    roots[2] = root_poly(reverse(num))
    errs[2] = eval_poly(cp, roots[2])
    // divide out to get linear
    den = {-roots[2], 1}
    num = reverse(poly_div(num,den))
    // find third root
    roots[3] = -num[1]
    errs[3] = 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 tests = {{{ 1, -1,  0},
                   { 0,  1,  1},
                   { 0,  0,  1}},
                  {{-2, -4,  2},
                   {-2,  1,  2},
                   { 4,  2,  5}}}
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,-5,6}
errors: {0,0,0}

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-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 polynomial after each root is found.

import "./matrix" for Matrix
import "./math" for Math
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])
    }
}

var eigenvalues = Fn.new { |m|
    var roots = List.filled(3, 0)
    var errs  = List.filled(3, 0)
    // 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[0] = Math.rootPoly(cp)
    errs[0] = Math.evalPoly(cp, roots[0])
    // divide out to get quadratic
    var num = Polynom.new(cp[-1..0])
    var den = Polynom.new([-roots[0], 1])
    num = (num/den)[0]
    // find second root
    roots[1] = Math.rootPoly(num.factors[-1..0])
    errs[1] = Math.evalPoly(cp, roots[1])
    // divide out to get linear
    den = Polynom.new([-roots[1], 1])
    num = (num/den)[0]
    // find third root
    roots[2] = -num.factors[0]
    errs[2] = Math.evalPoly(cp, roots[2])
    return [cp, roots, errs]
}

var mats = [
    Matrix.new ( [
        [ 1, -1,  0],
        [ 0,  1,  1],
        [ 0,  0,  1]
    ]),
    Matrix.new ([
        [-2, -4,  2],
        [-2,  1,  2],
        [ 4,  2,  5]
    ])
]

for (m in mats) {
    System.print("For matrix:")
    Fmt.mprint(m, 2, 0)
    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, -5, 6]
and the corresponding errors are: [0, 0, 0]