cubic solver

[Albert Chan]

Previous posts assumed depressed cubic. Here is how to depress a cubic.
This version, if cubic coefficients are integers, so does depressed cubic.

Let's define {a,b,c,d} = roots of polynomial, a*x^3 + b*x^2 + c*x + d
To force b to 0, we need substitution, x = y - b/A, where A=3a

But we don't want fractions. Let's scale the problem. Bonus: cubic coefficient = 1

{a, b, c, d} = {A, 3b, 3c, 3d} = {1, 3b, 3Ac, 3A²d} / A

Synthetic Division with offset b will clear 3b coefficient

Code:

-b>   1   3b    3Ac          3A²d
-b>   1   2b    3Ac-2b²      D = 3A²d-b*(3Ac-2b²)
-b>   1    b    C = 3Ac-3b²
      1    0

{a,b,c,d} = ({1, 0, C, D} - b) / A

Example, super-golden ratio, ψ³ = ψ²+1, we have

a, b, c, d = 1, -1, 0, -1

A = 3 * a = 3
C = 3 * (Ac-b²) = 3 * -1
D = 3A²d - b*(C+b²) = -29

α³, β³ = {1, D, -(C/3)³} = {1,-29,1} = (29 ± √(93²-4))/2 = (29 ± 3√93)/2

ψ = ({1, 0, -3, -29} - (-1)) / 3 = (α + β + 1) / 3

Cas> (cubic(3,29) .+ 1) / 3.

[1.46557123188, −0.232785615938+0.792551992515*i, −0.232785615938-0.792551992515*i]

This matched Michael Penn's ψ solution

(08-07-2022 07:55 AM)Thomas Klemm Wrote:  Recently I stumbled upon this video by Micheal Penn: What is the super-golden ratio??

It is defined by this equation:

\(
\begin{align}
\psi^3=\psi^2+1
\end{align}
\)

He solves it algebraically to find:

\(
\begin{align}
\psi=\frac{1}{3}\left(1+\sqrt[3]{\frac{1}{2}\left(29-3\sqrt{93}\right)}+\sqrt[3]{\frac{1}{2}\left(29+3\sqrt{93}\right)}\right)
\end{align}
\)

Code:

function cubic(a,b,c,d)         -- a*x^3 + b*x^2 + c*x + d = 0
    a = 3*a
    c = b*b-a*c
    d = 3*a*a*d - b*(b*b-3*c)   -- y coefs = 1, 0, -3*c, d
    d = quadratic(1, d, c^3)    -- x roots = ((y roots)-b)/a
    d = d:imag()==0 and cbrt(d:real()) or d^(1/3) -- beta
    c = I.abs(d)>0 and c/d or 0 -- alpha
    c, d = I(-.5)*(d+c), I(0,sqrt(.75))*(d-c)
    return (-c-c-b)/a, (c+d-b)/a, (c-d-b)/a
end    

function quadratic(a,b,c)       -- a*x^2 + b*x + c = 0
    b = -.5*b
    local y = I.sqrt(b*b-a*c)
    y = I.abs(b-y) > I.abs(b+y) and b-y or b+y
    return y/a, (I.abs(y)==0 and y or c/y)
end