+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: HP Prime (/forum-5.html)
+--- Thread: cubic solver (/thread-20476.html)
cubic solver - Albert Chan - 09-08-2023 02:26 PM
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
RE: cubic solver - Albert Chan - 09-08-2023 02:44 PM
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)\)
RE: cubic solver - parisse - 09-08-2023 06:56 PM
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/2f228406dbe099aaf7abf8d6b3fdb6ea8fd94c9d
RE: cubic solver - Albert Chan - 06-29-2024 02:46 PM
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
RE: cubic solver - Albert Chan - 06-29-2024 03:29 PM
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;
#endIs there a way to print out discriminant on the same screen? (i.e. without pressing key to get back?)
RE: cubic solver - Albert Chan - 07-02-2024 02:34 PM
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]
RE: cubic solver - Albert Chan - 07-02-2024 05:12 PM
(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]
RE: cubic solver - Albert Chan - 07-03-2024 10:17 PM
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)
endRE: cubic solver - Albert Chan - 07-03-2024 11:48 PM
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;
#endXCas
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];
}Quarter Solver - Albert Chan - 07-04-2024 08:37 PM
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
endlua> 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)