cubic solver

[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