cubic solver

[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]