Post Reply 
cubic solver
09-08-2023, 02:26 PM (This post was last modified: 03-03-2024 11:58 PM by Albert Chan.)
Post: #1
cubic solver
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
Find all posts by this user
Quote this message in a reply
09-08-2023, 02:44 PM
Post: #2
RE: cubic solver
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)\)
Find all posts by this user
Quote this message in a reply
09-08-2023, 06:56 PM
Post: #3
RE: cubic solver
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
Find all posts by this user
Quote this message in a reply
06-29-2024, 02:46 PM (This post was last modified: 07-02-2024 11:06 AM by Albert Chan.)
Post: #4
RE: cubic solver
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
Find all posts by this user
Quote this message in a reply
06-29-2024, 03:29 PM (This post was last modified: 07-02-2024 11:12 AM by Albert Chan.)
Post: #5
RE: cubic solver
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?)
Find all posts by this user
Quote this message in a reply
07-02-2024, 02:34 PM
Post: #6
RE: cubic solver
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]
Find all posts by this user
Quote this message in a reply
07-02-2024, 05:12 PM
Post: #7
RE: cubic solver
(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]
Find all posts by this user
Quote this message in a reply
07-03-2024, 10:17 PM (This post was last modified: 07-05-2024 09:38 PM by Albert Chan.)
Post: #8
RE: cubic solver
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
Find all posts by this user
Quote this message in a reply
07-03-2024, 11:48 PM
Post: #9
RE: cubic solver
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];
}
Find all posts by this user
Quote this message in a reply
07-04-2024, 08:37 PM (This post was last modified: 07-06-2024 03:20 PM by Albert Chan.)
Post: #10
Quarter Solver
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)
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 3 Guest(s)