HP50g simplifing a root

[Albert Chan]

I was being stupid. Why solve cubic equation for a
This is the way to get cube roots simplified:

\(\sqrt[3]{A ± \sqrt{R}} = a ± \sqrt{r} \quad ⇒ \quad a = \large\frac{\sqrt[3]{A+\sqrt{R}}
\;+\; \sqrt[3]{A-\sqrt{R}}}{2} \)

Round a to closest halves, get r, and double check if it round-trip back to (A, R)

Code:

function simp_cbrt4(A,B,k)  -- simplify cbrt(A + B * sqrt(k))
    local R = B*B*k
    local a = sqrt(abs(R))
    if R<0 then a = (A*A-R)^(1/6) * cos(atan2(a,A)/3)
    else        a = (cbrt(A+a) + cbrt(A-a))/2
    end
    a = a + 0x3p50 - 0x3p50 -- round to halves
    local r = (A/a-a*a)/3
    local b = B/(3*a*a+r)
    if r == b*b*k then return a,b,k end
end

lua> simp_cbrt4(1859814842094, -59687820010, 415)
11589      -145         415
lua> simp_cbrt4(300940299,103940300,101)
99           100          101
lua> simp_cbrt(180, 23, 157)   -- (a,b) can be halves
1.5          0.5          157
lua> simp_cbrt4(-36, 20, -7)    -- work with complex roots too.
3             1             -7
lua> simp_cbrt4(81,30,-3)       -- this is simplest, see comment from previous post.
4.5          0.5          -3

[Albert Chan]

(10-09-2020 05:21 PM)Albert Chan Wrote:  \(\sqrt[3]{A ± \sqrt{R}} = a ± \sqrt{r} \quad ⇒ \quad a = \large\frac{\sqrt[3]{A+\sqrt{R}}
\;+\; \sqrt[3]{A-\sqrt{R}}}{2} \)

Slightly off topics. I was curious how to calculate a accurately (assumed R > 0)
Since cube root is odd function, we pull out the sign.

a = sign(A)/2 * (³√(|A|+√R) + ³√(|A|-√R))

Last term can be rewritten as ³√(A²-R) / ³√(|A|+√R)

Let c = ³√(A²-R), d = ³√(|A|+√R)², we have:

a = sign(A)/2 * (d + c) / √d

Since both c, d are cube roots, we try identity c³ + d³ = (c + d) (c² - c d + d²)

c³ + d³ = (A² - R) + (A² + R + 2|A|√R) = |2A| * (|A|+√R) = |2A| * d√d

Substitute back in, many terms get cancelled, a = A / (c²/d - c + d)

If ³√(A ± √R) can be simplified, c turns to integer (we had shown this earlier).
In fact, the whole denominator turns to integer = a² + 3r = 4a² - 3c

Code:

function calc_a(A,R)
    local c = cbrt(A*A-R)
    local d = cbrt(A*A+R + 2*abs(A)*sqrt(R))
    return A / (c*c/d - c + d)
end

lua> function test_a(A,R)
:           local q = sqrt(R)
:           return calc_a(A,R), (cbrt(A+q) + cbrt(A-q))/2
:      end
lua>
lua> test_a(26, 15^2 * 3)
2      1.9999999999999982
lua> test_a(9416, 4256^2 * 5)
11    11.000000000000005
lua> test_a(300940299,103940300^2 * 101)
99    99.00000000000006

[Albert Chan]

(10-10-2020 03:58 PM)Albert Chan Wrote:  Substitute back in, many terms get cancelled, a = A / (c²/d - c + d)

If ³√(A ± √R) can be simplified, c turns to integer (we had shown this earlier).
In fact, the whole denominator turns to integer = a² + 3r = 4a² - 3c

Proof is trivial:

if ³√(A ± √R) = a ± √r, then

c = ³√(A² - R) = a² - r
d = ³√(|A|+√R)² = (|a|+√r)²

c²/d + d - c = (|a|-√r)² + (|a|+√r)² - c = (2a² + 2r) - (a² - r) = a² + 3r

[Albert Chan]

(10-09-2020 02:31 PM)Albert Chan Wrote:  Lets rearrange the cubic to match form x³ + 3px - 2q = 0

c³ = A² - R
4*a³ = 3*c*a + A
a³ + 3*(-c/4)*a - 2*(A/8) = 0

→ Cubic discriminant = p³ + q² = (-c/4)³ + (A/8)² = (-c³ + A²) / 64 = R/64
→ If R < 0, we got 3 real roots.

Instead of cubic discriminant, we can show this with more familiar quadratic discriminant.

Let f(x) = 4x³ - 3*c*x - A, so that f(a) = 0

f(x) = 4 * (x-a) * (x² + a*x + (a² - 3*c/4))

Quadratic discriminant = a² - 4*(a² - 3*c/4) = -3*(a² - c) = -3*r
→ if r < 0 (due to R < 0), f got 3 real roots.

We can solve the quadratics, but unnecessary.
If we have 1 solution, where (a ± √r)³ = A ± √R, we can get the other 2.

(a ± √r)³ = ((a ± √r) × w)³ = ((a ± √r) / w)³,     where w = (-1)^(2/3)

Example, simplify ³√(36 + 20i√7)

>>> from mpmath import *
>>> A, B, k = 36, 20, -7
>>> z = A + B*sqrt(k)
>>> q = cbrt(z); print q
(3.79128784747792 + 1.27520055582102j)

>>> c = cbrt(A*A-B*B*k); print c
16.0

q looks hopeless to simplify, but c is integer, thus possible.
Perhaps the "nice" root is not the principle root ?

>>> w = exp(2j/3*pi)
>>> print q*w, q/w
(-3.0 + 2.64575131106459j)      (-0.791287847477919 - 3.92095186688561j)

>>> a = -3
>>> f = lambda x: 4*x**3 - 3*c*x - A
>>> print f(a)
0.0
>>> print a, (A/a - a*a)/3
-3       -7.0

We have ³√(36 + 20i√7) = (-3 + i√7) / (-1)^(2/3)

However, converted to form a ± √r, expression look messier than original ³√(A + √R)

[Albert Chan]

(10-11-2020 06:28 PM)Albert Chan Wrote:  Example, simplify ³√(36 + 20i√7)

>>> from mpmath import *
>>> A, B, k = 36, 20, -7
...
>>> c = cbrt(A*A-B*B*k); print c
16.0

Here is a novel way to solve for r instead, then get a.
(Note: b can be halves, just like a)

R/r = (B/b)² = (3c + 4r)² = integer, where r = b²k

For this case, we hit the jackpot. If b=±1, 3c+4r = 3c+4k = B

r = k = -7
a = A / (c+4r) = 36/ (16-28) = -3

3*arg(a±√r) = 3*atan2(±√7, -3) ≈ ±2.3098 pi

We need to match above angle to arg(A+√R), which is on the first quadrant.
+2.3098 pi - 2 pi = 0.3098 pi matches.

³√(36 + 20i√7) = (-3 + i√7) / (-1)^(2/3)

This is just for fun. No need to solve the cubics.
Directly calculate ³√(A+√R) (previous post) is preferred.

Comment: Solving for b instead, we can take square root of both side

B/b = 3c + 4r = 3c + 4kb² = integer

This removed the ambiguity of b = ±1 to b = B/(3c + 4kb²) = 1
Above example need to match angles anyway, so I skipped this.

[peacecalc]

Hello Albert,

I'm really stunnig about your investigations in this topic. Wow, in future I should be more careful in arising such topics! ;-)

[Albert Chan]

(10-09-2020 02:31 PM)Albert Chan Wrote:  Lets rearrange the cubic to match form x³ + 3px - 2q = 0

c³ = A² - R
4*a³ = 3*c*a + A
a³ + 3*(-c/4)*a - 2*(A/8) = 0

→ Cubic discriminant = p³ + q² = (-c/4)³ + (A/8)² = (-c³ + A²) / 64 = R/64
→ If R < 0, we got 3 real roots.

see https://proofwiki.org/wiki/Cardano%27s_Formula

I just realized my direct calculation for a is really using Cardano's formula

a = ³√(q + √Δ) + ³√(q - √Δ)
   = ³√((A/8) + √(R/64)) + ³√((A/8) - √(R/64))
   = (³√(A+√R) + ³√(A-√R)) / 2

[CMarangon]

Hello!

Well, I conclude that HP50 does not have the formulas below, into memory.
Maybe, it can be included, in next ROM update.

Formula:
(A ± B⋅√k)^(1/3) = a ± b⋅√k

See also post #22, from Albert Chan, above.


Carlos (BR)

(10-09-2020 02:31 PM)Albert Chan Wrote:  

(10-09-2020 12:20 AM)Albert Chan Wrote:  Unfortunately, our divisors based cube root simplify routines were flawed !
Real rational root for a is a divisor of n, *but* possibly divided by 4 ...

Another problem is we had assumed there is only 1 real root for a.
Lets rearrange the cubic to match form x³ + 3px - 2q = 0

c³ = A² - R
4*a³ = 3*c*a + A
a³ + 3*(-c/4)*a - 2*(A/8) = 0

→ Cubic discriminant = p³ + q² = (-c/4)³ + (A/8)² = (-c³ + A²) / 64 = R/64
→ If R < 0, we got 3 real roots.

see https://proofwiki.org/wiki/Cardano%27s_Formula

Quote:XCas> find_cbrt(81,30,-3)       → \(\frac{9}{2} + \frac{i}{2}\cdot\sqrt{3}\)

XCas> find_all_a(A,R) := solve(surd(A*A-R,3) = a*a - (A/a-a*a)/3 , a)

XCas> find_all_a(81, 30^2 * -3)       → [-3, -3/2, 9/2]
XCas> simplify( [-3+2i*sqrt(3), -3/2-5i/2*sqrt(3), 9/2+i/2*sqrt(3)] .^ 3)

\(\qquad\quad [81 + 30i\sqrt3\;,\; 81 + 30i\sqrt3\;,\; 81 + 30i\sqrt3 \)

If we consider integer as simplest, then \(\sqrt[3]{81 + 30i\sqrt3} = -3+2i\sqrt3\)

Comment: I was wrong on above example.

x³ = y does not imply ³√y = x. We need to consider branch-cut
In other words, simplified form should match numerical evaluated values.

>>> (81 + 30j * 3**0.5) ** (1/3)
(4.5+0.86602540378443837j)

[peacecalc]

Hello Albert,

Quote:a = ³√(q + √Δ) + ³√(q - √Δ)
   = ³√((A/8) + √(R/64)) + ³√((A/8) - √(R/64))
   = (³√(A+√R) + ³√(A-√R)) / 2

This was the origin of my thoughts for simplyfing, because I construct examples to test cardanos formular i keyed in my HP50g and I got that kind of "garbage" for simple solutions i. e. 1 minus root of two. So my question arises, why my calc cannot simplify this...the rest is history in nearly 30 posts...

[peacecalc]

Hello CMaragon,

don't wait for a rom update, I'm feared you don't live long enough, because the HP50g is no longer produced and so the last rom version is 2.15. It's a pity, but for HP calc department is "prime time".

[grsbanks]

(10-12-2020 11:45 PM)CMarangon Wrote:  Well, I conclude that HP50 does not have the formulas below, into memory.
Maybe, it can be included, in next ROM update.

Sadly the 50g was discontinued 5 years ago. There will be no ROM update

[Albert Chan]

(10-11-2020 06:28 PM)Albert Chan Wrote:  We have ³√(36 + 20i√7) = (-3 + i√7) / (-1)^(2/3)

Mathematica gives many "simplifed" forms, here are 2 of them.

(-3+i√7) * (-1-i√3)/2 = (√21+3)/2 + i/2*(3√3-√7) = (√21+3)/2 + i/2*√(34-6√21)

I was curious, given √(34-6√21), how to "simplify" to (3√3-√7) ?

(10-09-2020 05:21 PM)Albert Chan Wrote:  \(\sqrt[3]{A ± \sqrt{R}} = a ± \sqrt{r} \quad ⇒ \quad a = \large\frac{\sqrt[3]{A+\sqrt{R}}
\;+\; \sqrt[3]{A-\sqrt{R}}}{2} \)

We can use the same trick (if LHS is real, for both ±)

\( \sqrt{A ± B\sqrt{c d}} = a\sqrt{c} \,±\, b\sqrt{d} \quad ⇒ \quad
a\sqrt{c} = \large\frac{\sqrt{A+B\sqrt{c d}}\;+\;\sqrt{A-B\sqrt{c d}}}{2} \)

lua> A, B, k = 34, -6, 21
lua> C = B * sqrt(k)
lua> mean = (sqrt(A+C) + sqrt(A-C))/2
lua> mean^2
27

27 = 3*3*3 = a*a*c. We try a=c=3, d=7, so just solve for b:

(3√3 + b√7)² = (27+7b²) + 6b√21 = 34 - 6√21

We have b = -1, thus √(34 - 6√21) = 3√3 - √7

It is nice that XCas do this automatically

XCas> simplify(√(34-6*√(21)))      → 3√3 - √7