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