cubic solver

[Albert Chan]

This program solve x^3 = a*x + b, for x

Code:

cubic_ab(a,b) := 
BEGIN
LOCAL d;
a, b := a/3, b/2;    /* x^3 = (3a)*x + (2b) */
d := sqrt(b^2-a^3) * (-1)^(sign(re(b))==-1);
b := surd(b+d, 3);
return [a/b, b];     /* alpha, beta */
END;

Code:

cubic(a,b) :=
BEGIN
LOCAL w := exp(2*pi/3*i);
[[1,1],[w,conj(w)],[conj(w),w]] * cubic_ab(a,b);
END;

Michael Penn's example




Let x = sin(θ) + cos(θ), it reduced to cubic: x^3 = 3*x - 11/8

Cas> r := cubic(3, -11/8)
Cas> simplify(r)                 → [1/2,    (3*√5-1)/4,       (-3*√5-1)/4]
Cas> float(Ans)                  → [0.5, 1.42705098312, −1.92705098312]
Cas> float(r)                      → [1.42705098312, −1.92705098312, 0.5]

We have Cas simplify bug here! symbolic and approximate numbers don't match.
Luckily, the algorithm doesn't care which cube roots is used (no need for principle root).

But, this bug should be fixed.

Cas> b := cubic_ab(3, -11/8) [2]     → (1/256*(-√15*48*i-176))^(1/3)
Cas> polar(float(b))                         → 1., −0.77627903074
Cas> polar(simplify(b))                    → 1 , atan(√15)

Simplify bug had phase angle off by 2*pi/3:

Cas> atan(√15) - 2*pi/3.0               → −0.77627903074

Update: cubic_ab(a,b) code explained here

[Albert Chan]

Root simplify returned wrong (not principle) root bug existed for XCas 1.5.0

But, it is worse on XCas 1.9.0-31 (mingw32 win32)
Result does not match phase *and* magnitude.

XCas> b := (1/256*(-√15*48*i-176))^(1/3)
XCas> polar(float(b))

1.0, -0.77627903074

XCas> simplify(polar(b))                     // good

\(\displaystyle 1\;,\;\frac{-\pi +\arctan \left(\frac{3}{11} \sqrt{15}\right)}{3}\)

XCas> simplify(polar(simplify(b)))       // bad

\(256\;,\;-\pi +\arctan \left(\frac{3}{11} \sqrt{15}\right)\)

[parisse]

Yes, there is an additional check in recent versions of Xcas, and here it failed because an argument was not reduced mod 2*pi.
https://github.com/geogebra/giac/commit/...ea8fd94c9d

[Albert Chan]

Updated cubic solver (x^3=a*x+b), use quadratic solver (x^2=a*x+b) to get alpha^3

Code:

quadratic(a,b) := { local d;
a /= 2; /* x^2 = 2a*x + b */
print(d := normal(a*a+b));
if (d==0) return [a,a];
d := sqrt(d);
a -= d * (-1) ^ bool(abs(a+d) >= abs(a-d));
return [a, -b/a];
}

Code:

cubic(a,b) := {
a /= 3; /* x^3 = 3a*x + b */
b := quadratic(b,-a^3)[0];
if (b==0) return [0,0,0];
b := surd(b,3)*(-1)^[0,2/3,-2/3]; /* betas */
return b .+ a ./ b;               /* roots  */
}

(06-29-2024 12:54 AM)Albert Chan Wrote:  x^3 - 4*x^2 + 8*x - 8 = 0            // let x = y + 4/3
y^3 + 8/3*y - 56/27 = 0               // c = 8/3, d = 56/27, to match above formula
...
x1 = 4/3 + (4/3)     + (-2/3)     = 2
x2 = 4/3 + (4/3)*ω + (-2/3)/ω = 1 + i√3
x3 = 4/3 + (4/3)/ω + (-2/3)*ω = -1/3 − i√3

XCas> cubic(-8/3, 56/27) .+ 4/3.
16/9

[2.0, 1.0+1.73205080757*i, 1.0-1.73205080757*i]

Note: (green) cubic discriminant is really quadratic discriminant.

Update: both solver now work even if both a, b are 0's

[Albert Chan]

Previous post was for XCas. Here is for HP Prime Cas

Code:

#cas
quadratic(a,b) := 
BEGIN
LOCAL d;
a /= 2; /* x^2 = 2a*x + b */
print(d := normal(a*a+b));
IF d==0 THEN return [a,a] END;
d := sqrt(d);
a -= d * (-1) ^ bool(abs(a+d) >= abs(a-d));
return [a, -b/a];
END;

cubic(a,b) :=
BEGIN
a /= 3; /* x^3 = 3a*x + b */
b := quadratic(b,-a^3)[1];
IF b==0 THEN return [0,0,0] END;
b := surd(b,3)*(-1)^[0,2/3,-2/3]; 
return b .+ a ./ b;
END;
#end

Is there a way to print out discriminant on the same screen? (i.e. without pressing key to get back?)

[Albert Chan]

Here is solving cubic with trig identities

cos(3θ) = 4*cos(θ)^3 - 3*cos(θ)
sin(3θ) = 3*sin(θ) - 4*sin(θ)^3

Or, to make formula with same pattern

4*cos(θ)^3 = 3*cos(θ) + cos(-3θ)
4*sin(θ)^3  = 3*sin(θ)  +  sin(-3θ)

x^3 = a*x + b            // let x = k*y
k³*y³ = a*k*y + b        // multiply both side by 4/k³, assuming k≠0
4*y³ = 4a/k²*y + 4b/k³
4*y³ = 3*y + t             // where k = √(4a/3), t = 4b/k³

This matched both trig identities. We just need to fix when k=0 (i.e. a=0)
For cos version, replace sin/asin with cos/acos.

Code:

cubic_kt(a,b) := { a:=sqrt(4*a/3); return normal([a, b*4/a^3]) };

cubic_sin(a,b) := { /* solve for x^3 = a*x + b */
  if (a==0) return surd(b,3) * (-1)^[0,2/3,-2/3];
  a, b := cubic_kt(a,b);
  return a * sin(([0,2*pi,-2*pi] .- asin(b)) / 3);
};

x^3 - 4*x^2 + 8*x - 8 = 0      // x = y + 4/3
y^3 = -8/3*y + 56/27

XCas> cubic_sin(-8/3, 56/27) .+ 4/3.

[2.0, 1.0+1.73205080757*i, 1.0-1.73205080757*i]

[Albert Chan]

(07-02-2024 02:34 PM)Albert Chan Wrote:  For cos version, replace sin/asin with cos/acos.

For numerical calculations, sin version perhaps is better. (is there exceptions?)
But there is a use for cosine, if we want to get real (α, β) from real (k, t)

cos(real or purely imaginary) = real

y^3 = a*y+b,      k = √(4a/3),      t = 4b/k³,      θ = acos(t)/3

y = k* cos(2n/3*pi - θ)
   = k*(cos(2n/3*pi)*cos(θ) + sin(2n/3*pi)*sin(θ))
   = [cos(2n/3*pi), i*sin(2n/3*pi)] • [k*cos(θ), k*sin(θ)/i]

Compare this to radical (α, β) version

y = cis(2m/3*pi) * α + cis(-2m/3*pi) * β
   = [cos(2m/3*pi), i*sin(2m/3*pi)] • [α+β, α−β]

Assuming m=n (it doesn't matter if they are not the same), we have:

α+β = k*cos(θ)
α−β = -i*k*sin(θ)

2α = (α+β) + (α−β) = k*cis(-θ)
2β = (α+β) − (α−β) = k*cis(θ)

LHS = one of RHS pairs. (we don't know if m=n, or sort order of (α, β))

k/2*cis(±acos(t)/3) = {α, β} | {α*ω, β/ω} |{α/ω, β*ω} | {β, α} | {β*ω, α/ω} | {β/ω, α*ω}

Example, x^3 = x + 1

XCas> a, b := cubic_ab(1.,1.)       → [0.337726750973, 0.986991206271]
XCas> k, t := cubic_kt(1.,1.)        → [1.15470053838, 2.59807621135]
XCas> θ := acos(t)/3                   → 0.536211995193*i
XCas> k/2 * [cis(θ), cis(-θ)]         → [0.337726750973, 0.986991206271]

[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

[Albert Chan]

Based from previous post lua code (depressed cubic code not needed anymore)

HP Prime Cas

Code:

#cas
quadratic(a,b,c) := 
BEGIN
LOCAL d;
b /= -2;
d := normal(sqrt(b*b-a*c));
b -= d * (-1) ^ bool(abs(b+d) >= abs(b-d));
return [b/a, (b ? c/b : b)];
END;

cubic(a,b,c,d) :=
BEGIN
a *= 3;
c := b*b-a*c;
d := 3*a*a*d - b*(b*b-3*c);
d := quadratic(1, d, c^3)[1];
d := normal(surd(d,3));
c := d ? c/d : 0;
c, d := (d+c)/-2, sqrt(-3)/2*(d-c);
return [(-c-c-b)/a, (c+d-b)/a, (c-d-b)/a];
END;
#end

XCas

Code:

quadratic(a,b,c) := { local d;
b /= -2;
print(d := normal(b*b-a*c));
d := sqrt(d);
b -= d * (-1) ^ bool(abs(b+d) >= abs(b-d));
return [b/a, (b ? c/b : b)];
}

cubic(a,b,c,d) := {
a *= 3;
c := b*b-a*c;
d := 3*a*a*d - b*(b*b-3*c);
d := quadratic(1, d, c^3)[0];
d := normal(surd(d,3));
c := d ? c/d : 0;
c, d := (d+c)/-2, sqrt(-3)/2*(d-c);
return [(-c-c-b)/a, (c+d-b)/a, (c-d-b)/a];
}

[Albert Chan]

We can use same trick (post #8) to depress quartic, a*x^4 + b*x^3 + c*x^2 + d*x + e = 0

{a, b, c, d, e} = {(A = 4a), 4b, 4c, 4d, 4e} = {1, 4b, 4Ac, 4A²d, 4A³e} / A

Synthetic Division to force coefficient 4b to 0

Code:

-b>   1   4b    4Ac          4A²d               4A³e
-b>   1   3b    4Ac-3b²      4A²d-b*(4Ac-3b²)   4A³e-b*(4A²d-b*(4Ac-3b²)) = E
-b>   1   2b    4Ac-5b²      4A²d-b*(8Ac-8b²) = D
-b>   1    b    4Ac-6b² = C
      1    0

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

y^4 + C*y^2 + D*y + E = 0
(y^2+s)^2 = (2s-C)*y^2 - D*y + (s^2-E)      // Add an offset s, making LHS perfect square

Let k² = 2s-C ≠ 0, and force RHS as perfect square too.

{(k² = 2s-C), -D, (s²-E)} = {1, -D, (D²/4 = (s²-E)*k²)} / k²

D² = 4*(s²-E)*k² = ((k²+C)² - 4E) * k²

(k²)^3 + 2*C*(k²)^2 + (C²-4E)*(k²) - D² = 0

The reason we solve for k² instead of s is to easily check k²≠0
k² roots all zero ⇒ {C,D,E} all zeroes ⇒ y roots all zero too.

(y²+s)² = k² * (y-(D/2)/k²)²      // let t = D/2/k
(y²+s) = ±(k*y - t)

(y^2 - ky + s+t) * (y^2 + ky + s-t) = 0

(s±t), one of them may cause huge cancellation, we only use the good one.
The other can be derived by product of all roots = (s+t) * (s-t) = E

Last step is simply map y back to x: (x roots) = ((y roots) - b) / A

Code:

function quartic(a,b,c,d,e) -- a*x^4 + b*x^3 + c*x^2 + d*x + e = 0
    a = 4*a
    c, d, e = 4*a*c, 4*a^2*d, 4*a^3*e
    e = e - b*(d - b*(c - 3*b*b))
    d = d - b*(c+c - 8*b*b)
    c = c - 6*b*b           -- y coefs = 1, 0, c, d, e
    local k2 = {cubic(1, 2*c, c*c-4*e, -d*d)}
    table.sort(k2, function(a,b) return I.abs(a) > I.abs(b) end)
    if I.abs(k2[1])==0 then c=I(-b/a); return c,c,c,c end
    local k = I.sqrt(k2[1]) -- (y^2+s)^2 = (k*y-t)^2
    c, d = k2[1]+c, d/k     -- 2*s, 2*t
    c, d = (c+d)/2, (c-d)/2 -- quadratic constant term
    if I.abs(c) < I.abs(d) then c=d; k=-k end
    local y1, y2 = quadratic(1, -k,  c)
    local y3, y4 = quadratic(1, k, e/c)
    return (y1-b)/a, (y2-b)/a, (y3-b)/a, (y4-b)/a
end

lua> quartic(1,2,3,4,5)
(0.28781547955764775+1.4160930801719083*I)
(-1.287815479557648+0.8578967583284901*I)
(0.2878154795576481-1.4160930801719085*I)
(-1.287815479557648-0.8578967583284899*I)