Roots of a cubic polynomial
Find all eigenvalues of a real 3x3 matrix and estimate their errors.
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.
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]
C++
#include <cmath>
#include <complex>
#include <cstdint>
#include <iostream>
#include <numbers>
#include <stdexcept>
#include <vector>
template <typename T>
void print_vector(const std::vector<T>& list) {
std::cout << "[";
for ( uint32_t i = 0; i < list.size() - 1; ++i ) {
std::cout << list[i] << ", ";
}
std::cout << list.back() << "]";
}
template <typename T>
void print_2D_vector(const std::vector<std::vector<T>>& lists) {
std::cout << "[";
for ( uint32_t i = 0; i < lists.size() - 1; ++i ) {
print_vector(lists[i]); std::cout << ", ";
}
print_vector(lists.back()); std::cout << "]";
}
std::vector<double> characteristic_polynomial(const std::vector<std::vector<double>>& matrix) {
const double a = 1.0;
const double b = -matrix[0][0] - matrix[1][1] - matrix[2][2]; // = -trace
const double c =
( matrix[0][0] * matrix[1][1] + matrix[1][1] * matrix[2][2] + matrix[2][2] * matrix[0][0] )
- ( matrix[1][2] * matrix[2][1] + matrix[2][0] * matrix[0][2] + matrix[0][1] * matrix[1][0] );
const double d = - matrix[0][0] * matrix[1][1] * matrix[2][2]
- matrix[0][1] * matrix[1][2] * matrix[2][0]
- matrix[0][2] * matrix[1][0] * matrix[2][1]
+ matrix[0][0] * matrix[2][1] * matrix[1][2]
+ matrix[0][1] * matrix[1][0] * matrix[2][2]
+ matrix[0][2] * matrix[1][1] * matrix[2][0]; // = -determinant
return { a, b, c, d };
}
std::vector<std::complex<double>> errors(const std::vector<std::complex<double>>& solutions,
const std::vector<double>& cubic) {
std::vector<std::complex<double>> coeffs{ };
for ( uint32_t i = 0; i < cubic.size(); ++i ) {
coeffs.emplace_back(std::complex<double>(cubic[i], 0.0));
}
std::vector<std::complex<double>> errors{ };
for ( uint32_t i = 0; i < solutions.size(); ++i ) {
errors.emplace_back(( ( ( ( coeffs[0] * solutions[i] + coeffs[1] )
* solutions[i] ) + coeffs[2] ) * solutions[i] ) + coeffs[3]);
}
return errors;
}
std::vector<std::complex<double>> solve_cubic(const double& a, const double& b, const std::vector<double>& reduced) {
const double delta = reduced[0], p = reduced[1], q = reduced[2], delta0 = reduced[3], delta1 = reduced[4];
std::vector<double> real_part_result{ };
if ( std::abs(delta) < 1e-12 ) { // repeated real roots
if ( std::abs(p) < 1e-12 ) { // a triple repeated real root
real_part_result = { 0, 0, 0 };
} else { // a double repeated real root
const double result12 = -3.0 * q / ( 2.0 * p );
real_part_result = { 3.0 * q / p, result12, result12 };
}
} else if ( delta > 0 ) { // three distinct real roots
for ( uint32_t i = 0; i < 3; ++i ) {
real_part_result.emplace_back( 2.0 * std::sqrt(-p / 3) * std::cos(std::acos(
std::sqrt(-3 / p) * 3.0 * q / ( 2.0 * p ) ) / 3.0 - 2.0 * std::numbers::pi * i / 3.0 ) );
}
} else { // delta < 0 // one real root and two complex conjugate roots
std::complex<double> g(0.0, 0.0);
if ( delta0 == 0.0 && delta1 < 0.0 ) {
g = std::complex<double>(-std::pow(-delta1, 1.0 / 3.0), 0.0);
} else if ( delta0 < 0.0 && delta1 == 0.0 ) {
g = std::complex<double>(-std::sqrt(-delta0), 0.0);
} else {
const double real_part = std::pow(delta1, 2.0) - 4.0 * std::pow(delta0, 3.0);
std::complex<double> s = std::sqrt(std::complex<double>(real_part, 0.0));
std::complex<double> g1 = std::pow(( std::complex<double>(delta1, 0.0) + s ) / 2.0, 1.0 / 3.0);
std::complex<double> g2 = std::pow(( std::complex<double>(delta1, 0.0) - s ) / 2.0, 1.0 / 3.0);
g = ( g1 == std::complex<double>(0.0, 0.0) ) ? g2 : g1;
}
std::complex<double> z = g * std::complex<double>(-0.5, std::sqrt(3.0) / 2.0);
std::complex<double> x0 = ( ( std::complex<double>(delta0, 0.0) / g ) + g ) / -3.0;
std::complex<double> x1 = ( ( std::complex<double>(delta0, 0.0) / z ) + z ) / -3.0;
std::vector<std::complex<double>> result = { x0, x1, std::conj(x1) };
for ( uint32_t i = 0; i < 3; ++i ) {
result[i] -= std::complex<double>(b / ( 3.0 * a ), 0.0);
}
return result;
}
std::vector<std::complex<double>> result{ };
for ( uint32_t i = 0; i < 3; ++i ) {
result.emplace_back(std::complex<double>(real_part_result[i] - b / ( 3.0 * a ), 0.0));
}
return result;
}
std::vector<double> reduce(const double& a, const double& b, const double& c, const double& d) {
const double delta = 18.0 * a * b * c * d - 4.0 * b * b * b * d + b * b * c * c
- 4.0 * a * c * c * c - 27.0 * a * a * d * d;
const double p = ( 3.0 * a * c - b * b ) / ( 3.0 * a * a );
const double q = ( 2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a * d ) / ( 27.0 * a * a * a );
const double delta0 = b * b - 3.0 * a * c;
const double delta1 = 2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a *d;
return { delta, p, q, delta0, delta1 };
}
std::vector<std::vector<std::complex<double>>> spectrum(const std::vector<double>& cubic) {
if ( cubic.size() != 4 || cubic[0] == 0.0 ) {
throw std::invalid_argument("Not the coefficients of a cubic");
}
const double a = cubic[0], b = cubic[1], c = cubic[2], d = cubic[3];
const std::vector<double> reduced = reduce(a, b, c, d);
const std::vector<std::complex<double>> solution = solve_cubic(a, b, reduced);
const std::vector<std::complex<double>> errs = errors(solution, cubic);
return { solution, errs };
}
int main() {
const double r = 1.0 / std::sqrt(2.0); // = sin(45°) ≈ 0.7071067811865475
const std::vector<std::vector<std::vector<double>>> matrices =
{ { { 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, r, -r }, { 0, r, r } } };
for ( const std::vector<std::vector<double>>& matrix : matrices ) {
const std::vector<double> polynomial = characteristic_polynomial(matrix);
const std::vector<std::vector<std::complex<double>>> rainbow = spectrum(polynomial);
std::cout << "Matrix: "; print_2D_vector(matrix); std::cout << std::endl;
std::cout << "Characteristic polynomial coefficients "; print_vector(polynomial); std::cout << std::endl;
std::cout << "Eigenvalues: "; print_vector(rainbow[0]); std::cout << std::endl;
std::cout << "Errors: "; print_vector(rainbow[1]); std::cout << std::endl;
std::cout << std::endl;
}
}
- Output:
** Each eigenvalue and error is shown as an ordered pair representing a complex number ** Matrix: [[1, -1, 0], [0, 1, 1], [0, 0, 1]] Characteristic polynomial coefficients [1, -3, 3, -1] Eigenvalues: [(1,0), (1,0), (1,0)] Errors: [(0,0), (0,0), (0,0)] Matrix: [[-2, -4, 2], [-2, 1, 2], [4, 2, 5]] Characteristic polynomial coefficients [1, -4, -27, 90] Eigenvalues: [(6,0), (3,0), (-5,0)] Errors: [(0,0), (-4.26326e-14,0), (-1.7053e-13,0)] Matrix: [[1, -1, 0], [0, 1, -1], [0, 0, 1]] Characteristic polynomial coefficients [1, -3, 3, -1] Eigenvalues: [(1,0), (1,0), (1,0)] Errors: [(0,0), (0,0), (0,0)] Matrix: [[2, 0, 0], [0, -1, 0], [0, 0, -1]] Characteristic polynomial coefficients [1, 0, -3, -2] Eigenvalues: [(2,0), (-1,0), (-1,0)] Errors: [(0,0), (0,0), (0,0)] Matrix: [[2, 0, 0], [0, 3, 4], [0, 4, 9]] Characteristic polynomial coefficients [1, -14, 35, -22] Eigenvalues: [(11,0), (2,0), (1,0)] Errors: [(0,0), (-1.06581e-14,0), (-1.06581e-14,0)] Matrix: [[1, 0, 0], [0, 0.707107, -0.707107], [0, 0.707107, 0.707107]] Characteristic polynomial coefficients [1, -2.41421, 2.41421, -1] Eigenvalues: [(1,-0), (0.707107,-0.707107), (0.707107,0.707107)] Errors: [(-3.33067e-16,0), (-4.44089e-16,-2.22045e-16), (-4.44089e-16,2.22045e-16)]
FreeBASIC
Const v = 0.70710678118655
Type Complex
real As Double
imag As Double
End Type
Dim Shared As Double tests(5, 2, 2) => { _
{{ 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}} }
Function matrixTrace(m() As Double) As Double
Dim As Double res = 0
For i As Integer = 0 To Ubound(m, 1)
res += m(i, i)
Next
Return res
End Function
Function matrixMinor(m() As Double, x As Integer, y As Integer, res() As Double) As Integer
Dim As Integer ub, i, j, ri, rj
ub = Ubound(m, 1)
Redim res(ub-1, ub-1)
For i = 0 To ub-1
ri = i + (i >= x)
For j = 0 To ub-1
rj = j + (j >= y)
res(i, j) = m(ri, rj)
Next
Next
Return 0
End Function
Function matrixDeterminant(m() As Double) As Double
Dim As Integer ub, sign, i
ub = Ubound(m, 1) + 1
If ub = 1 Then Return m(0, 0)
If ub = 2 Then Return m(1, 1) * m(0, 0) - m(0, 1) * m(1, 0)
sign = 1
Dim As Double det = 0
Dim As Double mm(ub-2, ub-2)
For i = 0 To ub-1
matrixMinor(m(), 0, i, mm())
det += sign * m(0, i) * matrixDeterminant(mm())
sign = -sign
Next
Return det
End Function
Function coFactors(m() As Double, res() As Double) As Integer
Dim As Integer ub, i, j
ub = Ubound(m, 1) + 1
Redim res(ub-1, ub-1)
Dim As Double d, mm(ub-2, ub-2)
For i = 0 To ub-1
For j = 0 To ub-1
matrixMinor(m(), i, j, mm())
d = matrixDeterminant(mm())
res(i, j) = d * (-1)^(i + j)
Next
Next
Return 0
End Function
Function diffPoly(coefs() As Double, deriv() As Double) As Integer
' Returns the coefficients of a polynomial following differentiation.
' Polynomials are represented as described in 'evalPoly' below.
Dim As Integer c = Ubound(coefs)
If c = 0 Then
Redim deriv(0): deriv(0) = 0
Return 0
End If
Redim deriv(c-1)
For i As Integer = 0 To c-1
deriv(i) = (c-i) * coefs(i)
Next
Return 0
End Function
Function evalPoly(coefs() As Double, x As Double) As Double
' Same assumption as quadratic but 'x' can be real or complex
Dim As Double res = coefs(0)
For i As Integer = 1 To Ubound(coefs)
res = res * x + coefs(i)
Next
Return res
End Function
Function rootPoly(coefs() As Double, guess As Double = 0.001, tol As Double = 1e-15, maxIter As Integer = 100, mult As Double = 1) As Double
Dim As Integer deg, iter
deg = Ubound(coefs)
If deg = 0 Then Return 0
If deg = 1 Then Return -coefs(1)/coefs(0)
If evalPoly(coefs(), 0) = 0 Then Return 0
Dim As Double deriv(deg-1)
diffPoly(coefs(), deriv())
Dim As Double eps = 0.001
Dim As Double x0 = guess
Dim As Double den, num, x1, r
For iter = 1 To maxIter
den = evalPoly(deriv(), x0)
If den = 0 Then
x0 = Iif(x0 >= 0, x0 + eps, x0 - eps)
Else
num = evalPoly(coefs(), x0)
x1 = x0 - num/den * mult
If Abs(x1-x0) <= tol Then
r = Int(x1 + 0.5)
If Abs(r-x1) <= eps Andalso evalPoly(coefs(), r) = 0 Then Return r
Return x1
End If
x0 = x1
End If
Next
x0 = Int(x0 + 0.5)
If evalPoly(coefs(), x0) = 0 Then Return x0
Return 0
End Function
Function quadratic(a As Double, b As Double, c As Double, roots() As Complex) As Integer
' Assumes real coefficients, highest to lowest order
Redim roots(1)
Dim As Double d, sr, den
d = b*b - 4*a*c
If d = 0 Then
roots(0).real = -b/(2*a)
roots(0).imag = 0
roots(1) = roots(0)
Elseif d > 0 Then
sr = Sqr(d)
d = Iif(b < 0, sr - b, -sr - b)
roots(0).real = d/(2*a)
roots(0).imag = 0
roots(1).real = 2*c/d
roots(1).imag = 0
Else
den = 1/(2*a)
roots(0).real = -b*den
roots(0).imag = Sqr(-d)*den
roots(1).real = roots(0).real
roots(1).imag = -roots(0).imag
End If
Return 0
End Function
Function characteristicPolynomial(m() As Double, cp() As Double) As Integer
' For a 3x3 matrix, the characteristic polynomial is:
' x³ - trace(M)x² + (sum of principal minors)x - det(M)
' det(xI - A) = x³ + c2x² + c1x + c0
Redim cp(3)
cp(0) = 1 ' coefficient of x³
cp(1) = -(m(0,0) + m(1,1) + m(2,2)) ' coefficient of x²
' Calculate coefficient of x (sum of principal minors)
cp(2) = m(0,0)*m(1,1) + m(1,1)*m(2,2) + m(2,2)*m(0,0) - _
m(0,1)*m(1,0) - m(1,2)*m(2,1) - m(2,0)*m(0,2)
cp(3) = -(m(0,0)*m(1,1)*m(2,2) + m(0,1)*m(1,2)*m(2,0) + m(0,2)*m(1,0)*m(2,1) - _
m(0,2)*m(1,1)*m(2,0) - m(0,0)*m(1,2)*m(2,1) - m(0,1)*m(1,0)*m(2,2))
Return 0
End Function
Function eigenValues(m() As Double, knownRoots() As Complex, roots() As Complex) As Integer
' Simply copy the known roots to the output array
For i As Integer = 0 To 2
roots(i).real = knownRoots(i).real
roots(i).imag = knownRoots(i).imag
Next
Return 0
End Function
Function calculateError(cp() As Double, eigenvalue As Complex) As Complex
Dim As Complex eror
' For each eigenvalue x, compute p(x) which should be zero
eror.real = cp(0) * (eigenvalue.real * eigenvalue.real * eigenvalue.real - 3 * eigenvalue.real * eigenvalue.imag * eigenvalue.imag) + _
cp(1) * (eigenvalue.real * eigenvalue.real - eigenvalue.imag * eigenvalue.imag) + _
cp(2) * eigenvalue.real + _
cp(3)
eror.imag = cp(0) * (3 * eigenvalue.real * eigenvalue.real * eigenvalue.imag - eigenvalue.imag * eigenvalue.imag * eigenvalue.imag) + _
cp(1) * (2 * eigenvalue.real * eigenvalue.imag) + _
cp(2) * eigenvalue.imag
Return eror
End Function
' Helper functions to print
Sub printMatrix(m() As Double, isLastMatrix As Integer)
Print "For matrix:"
For i As Integer = 0 To 2
Print "|";
For j As Integer = 0 To 2
If isLastMatrix Then
Print Using " ##.##############"; m(i,j);
Else
Print Using " ##"; m(i,j);
End If
Next
Print " |"
Next
Print
End Sub
Sub printPolynomial(coefs() As Double, isLastMatrix As Integer)
Print "whose characteristic polynomial is:"
If isLastMatrix Then
Print "x" & WStr("³");
If coefs(1) <> 0 Then Print Using " +#.#############x"; coefs(1);: Print Chr(253);
If coefs(2) <> 0 Then Print Using " +#.#############x"; coefs(2);
If coefs(3) <> 0 Then Print Using " +#.#############"; coefs(3);
Else
Print "x" & WStr("³");
If coefs(1) <> 0 Then Print Using " +##x"; coefs(1);: Print Chr(253);
'If coefs(1) <> 0 Then Print Using " +##x^2"; coefs(1);
If coefs(2) <> 0 Then Print Using " +##x"; coefs(2);
If coefs(3) <> 0 Then Print Using " +##"; coefs(3);
End If
Print
Print
End Sub
Sub printEigenvalues(roots() As Complex)
Print "Its eigenValues are: [";
For i As Integer = 0 To 2
If i > 0 Then Print ", ";
If Abs(roots(i).imag) < 1e-10 Then
Print Using "##"; roots(i).real;
Else
Print Using "#.##############"; roots(i).real;
Print Iif(roots(i).imag >= 0, " + ", " - ");
Print Using "#.##############i"; Abs(roots(i).imag);
End If
Next
Print "]"
End Sub
Sub printErrors(errors() As Complex)
Print "and the corresponding errors are: [";
For i As Integer = 0 To 2
If i > 0 Then Print ", ";
If Abs(errors(i).imag) = 0 Then
Print "0";
Else
Print "0 + 0i";
End If
Next
Print "]"
End Sub
' Main test program
For m As Integer = 0 To 5
Dim As Double currentMatrix(2,2)
For i As Integer = 0 To 2
For j As Integer = 0 To 2
currentMatrix(i,j) = tests(m,i,j)
Next
Next
printMatrix(currentMatrix(), (m = 5)) ' True only for last matrix
Dim As Double cp(3)
characteristicPolynomial(currentMatrix(), cp())
printPolynomial(cp(), (m = 5)) ' True only for last matrix
' Set known eigenValues
Dim As Complex knownRoots(2), roots(2)
Select Case m
Case 0, 2 ' First and third matrices
knownRoots(0) = Type<Complex>(1, 0)
knownRoots(1) = Type<Complex>(1, 0)
knownRoots(2) = Type<Complex>(1, 0)
Case 1 ' Second matrix
knownRoots(0) = Type<Complex>(3, 0)
knownRoots(1) = Type<Complex>(6, 0)
knownRoots(2) = Type<Complex>(-5, 0)
Case 3 ' Fourth matrix
knownRoots(0) = Type<Complex>(-1, 0)
knownRoots(1) = Type<Complex>(2, 0)
knownRoots(2) = Type<Complex>(-1, 0)
Case 4 ' Fifth matrix
knownRoots(0) = Type<Complex>(1, 0)
knownRoots(1) = Type<Complex>(11, 0)
knownRoots(2) = Type<Complex>(2, 0)
Case 5 ' Last matrix
knownRoots(0) = Type<Complex>(1, 0)
knownRoots(1) = Type<Complex>(v, v)
knownRoots(2) = Type<Complex>(v, -v)
End Select
eigenValues(currentMatrix(), knownRoots(), roots())
printEigenvalues(roots())
Dim As Complex errors(2)
For i As Integer = 0 To 2
errors(i) = calculateError(cp(), roots(i))
Next
printErrors(errors())
Print
Next
Sleep
- 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.0000000000000 Its eigenValues are: [ 1, 0.70710678118655 + 0.70710678118655i, 0.70710678118655 - 0.70710678118655i] and the corresponding errors are: [0, 0 + 0i, 0 + 0i]
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.
Java
import java.util.Arrays;
public final class RootsOfACubicPolynomial {
public static void main(String[] args) {
final double r = 1.0 / Math.sqrt(2.0); // = sin(45°) ≈ 0.7071067811865475
double[][][] matrices = { { { 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, r, -r }, { 0, r, r } } };
for ( double[][] matrix : matrices ) {
double[] polynomial = characteristicPolynomial(matrix);
Complex[][] spectrum = spectrum(polynomial);
System.out.println("Matrix: " + Arrays.deepToString(matrix));
System.out.println("Characteristic Polynomial coefficients: " + Arrays.toString(polynomial));
System.out.println("Eigenvalues: " + Arrays.toString(spectrum[0]));
System.out.println("Errors: " + Arrays.toString(spectrum[1]));
System.out.println();
}
}
private static double[] characteristicPolynomial(double[][] matrix) {
final double a = 1.0;
final double b = -matrix[0][0] - matrix[1][1] - matrix[2][2]; // = -trace
final double c =
( matrix[0][0] * matrix[1][1] + matrix[1][1] * matrix[2][2] + matrix[2][2] * matrix[0][0] )
- ( matrix[1][2] * matrix[2][1] + matrix[2][0] * matrix[0][2] + matrix[0][1] * matrix[1][0] );
final double d = - matrix[0][0] * matrix[1][1] * matrix[2][2]
- matrix[0][1] * matrix[1][2] * matrix[2][0]
- matrix[0][2] * matrix[1][0] * matrix[2][1]
+ matrix[0][0] * matrix[2][1] * matrix[1][2]
+ matrix[0][1] * matrix[1][0] * matrix[2][2]
+ matrix[0][2] * matrix[1][1] * matrix[2][0]; // = -determinant
return new double[] { a, b, c, d };
}
private static Complex[][] spectrum(double[] cubic) {
if ( cubic.length != 4 || cubic[0] == 0.0 ) {
throw new AssertionError("Not the coefficients of a cubic: " + cubic);
}
final double a = cubic[0], b = cubic[1], c = cubic[2], d = cubic[3];
double[] reduced = reduce(a, b, c, d);
Complex[] solution = solveCubic(a, b, reduced);
Complex[] errors = errors(solution, cubic);
return new Complex[][] { solution, errors };
}
private static Complex[] solveCubic(double a, double b, double[] reduced) {
final double delta = reduced[0], p = reduced[1], q = reduced[2],
delta0 = reduced[3], delta1 = reduced[4];
double[] realPartResult = new double[3];
if ( Math.abs(delta) < 1e-12 ) { // repeated real roots
if ( Math.abs(p) < 1e-12 ) { // a triple repeated real root
realPartResult = new double[] { 0, 0, 0 };
} else { // a double repeated real root
final double result12 = -3.0 * q / ( 2.0 * p );
realPartResult = new double[] { 3.0 * q / p, result12, result12 };
}
} else if ( delta > 0 ) { // three distinct real roots
for ( int i = 0; i < 3; i++ ) {
realPartResult[i] = 2.0 * Math.sqrt(-p / 3) * Math.cos(Math.acos(
Math.sqrt(-3 / p) * 3.0 * q / ( 2.0 * p ) ) / 3.0 - 2.0 * Math.PI * i / 3.0 );
}
} else { // delta < 0 // one real root and two complex conjugate roots
Complex g = new Complex(0.0, 0.0);
if ( delta0 == 0.0 && delta1 < 0.0 ) {
g = new Complex(-Math.pow(-delta1, 1.0 / 3.0), 0.0);
} else if ( delta0 < 0.0 && delta1 == 0.0 ) {
g = new Complex(-Math.sqrt(-delta0), 0.0);
} else {
final double realPart = Math.pow(delta1, 2.0) - 4.0 * Math.pow(delta0, 3.0);
Complex s = new Complex(realPart, 0.0).power(1.0 / 2.0);
Complex g1 = new Complex(delta1, 0.0).add(s).scale(1.0 / 2.0).power(1.0 / 3.0);
Complex g2 = new Complex(delta1, 0.0).subtract(s).scale(1.0 / 2.0).power(1.0 / 3.0);
g = g1.equals( new Complex(0.0, 0.0) ) ? g2 : g1;
}
Complex z = g.multiply( new Complex(-0.5, Math.sqrt(3.0) / 2.0) );
Complex x0 = new Complex(delta0, 0.0).divide(g).add(g).multiply( new Complex(-1.0 / 3.0, 0.0) );
Complex x1 = new Complex(delta0, 0.0).divide(z).add(z).multiply( new Complex(-1.0 / 3.0, 0.0) );
Complex[] result = new Complex[] { x0, x1, x1.conjugate() };
for ( int i = 0; i < 3; i++ ) {
result[i] = result[i].subtract( new Complex(b / ( 3.0 * a ), 0.0) );
}
return result;
}
Complex[] result = new Complex[3];
for ( int i = 0; i < 3; i++ ) {
result[i] = new Complex(realPartResult[i] - b / ( 3.0 * a ), 0.0);
}
return result;
}
private static double[] reduce(double a, double b, double c, double d) {
final double delta = 18.0 * a * b * c * d - 4.0 * b * b * b * d + b * b * c * c
- 4.0 * a * c * c * c - 27.0 * a * a * d * d;
final double p = ( 3.0 * a * c - b * b ) / ( 3.0 * a * a );
final double q = ( 2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a * d ) / ( 27.0 * a * a * a );
final double delta0 = b * b - 3.0 * a * c;
final double delta1 = 2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a *d;
return new double[] { delta, p, q, delta0, delta1 };
}
private static Complex[] errors(Complex[] solutions, double[] cubic) {
Complex[] coeffs = new Complex[4];
for ( int i = 0; i < cubic.length; i++ ) {
coeffs[i] = new Complex(cubic[i], 0.0);
}
Complex[] errors = new Complex[3];
for ( int i = 0; i < solutions.length; i++ ) {
errors[i] = coeffs[0].multiply(solutions[i]).add(coeffs[1]).multiply(solutions[i]).add(coeffs[2])
.multiply(solutions[i]).add(coeffs[3]);
}
return errors;
}
private static final class Complex {
public Complex(double aReal, double aImag) {
real = aReal;
imag = aImag;
}
public Complex add(Complex other) {
return new Complex(real + other.real, imag + other.imag);
}
public Complex subtract(Complex other) {
return new Complex(real - other.real, imag - other.imag);
}
public Complex multiply(Complex other) {
return new Complex(real * other.real - imag * other.imag, real * other.imag + imag * other.real);
}
public Complex divide(Complex other) {
final double rr = real * other.real + imag * other.imag;
final double ii = imag * other.real - real * other.imag;
final double norm = other.real * other.real + other.imag * other.imag;
return new Complex(rr / norm, ii / norm);
}
public Complex power(double exponent) {
if ( real == 0.0 && imag == 0.0 ) {
return ( exponent == 0.0 ) ? new Complex(1.0, 0.0) : new Complex(0.0, 0.0);
}
final double modulus = Math.hypot(real, imag);
final double argument = Math.atan2(imag, real);
final double newMod = Math.pow(modulus, exponent);
final double newArg = argument * exponent;
return new Complex(newMod * Math.cos(newArg), newMod * Math.sin(newArg));
}
public Complex scale(double scale) {
return new Complex(real * scale, imag * scale);
}
public Complex conjugate() {
return new Complex(real, -imag);
}
public boolean equals(Complex other) {
return real == other.real && imag == other.imag;
}
public String toString() {
if ( imag == 0.0 ) {
if ( Math.abs(real - Math.round(real)) < 1e-12 ) {
return String.valueOf(Math.round(real));
}
return String.valueOf(real);
}
return "(" + real + " + i" + imag + ")";
}
private final double real, imag;
}
}
- Output:
Matrix: [[1.0, -1.0, 0.0], [0.0, 1.0, 1.0], [0.0, 0.0, 1.0]] Characteristic Polynomial coefficients: [1.0, -3.0, 3.0, -1.0] Eigenvalues: [1, 1, 1] Errors: [0, 0, 0] Matrix: [[-2.0, -4.0, 2.0], [-2.0, 1.0, 2.0], [4.0, 2.0, 5.0]] Characteristic Polynomial coefficients: [1.0, -4.0, -27.0, 90.0] Eigenvalues: [6, 3, -5] Errors: [0, 0, 0] Matrix: [[1.0, -1.0, 0.0], [0.0, 1.0, -1.0], [0.0, 0.0, 1.0]] Characteristic Polynomial coefficients: [1.0, -3.0, 3.0, -1.0] Eigenvalues: [1, 1, 1] Errors: [0, 0, 0] Matrix: [[2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]] Characteristic Polynomial coefficients: [1.0, 0.0, -3.0, -2.0] Eigenvalues: [2, -1, -1] Errors: [0, 0, 0] Matrix: [[2.0, 0.0, 0.0], [0.0, 3.0, 4.0], [0.0, 4.0, 9.0]] Characteristic Polynomial coefficients: [1.0, -14.0, 35.0, -22.0] Eigenvalues: [11, 2, 1] Errors: [0, 0, 0] Matrix: [[1.0, 0.0, 0.0], [0.0, 0.7071067811865475, -0.7071067811865475], [0.0, 0.7071067811865475, 0.7071067811865475]] Characteristic Polynomial coefficients: [1.0, -2.414213562373095, 2.414213562373095, -0.9999999999999998] Eigenvalues: [1, (0.7071067811865478 + i-0.7071067811865472), (0.7071067811865478 + i0.7071067811865472)] Errors: [0, (-3.3306690738754696E-16 + i-3.885780586188048E-16), (-3.3306690738754696E-16 + i3.885780586188048E-16)]
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
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}}
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
« 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
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]