# HP Forums

Full Version: HP50g simplifing a root
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Pages: 1 2
Hello folks,

I was wondering why my good old hp50g is not able to simplify:

$\left(26 - 15\cdot\sqrt{3} \right)^{\frac{1}{3}} = 2 - \sqrt{3}$

Hmm.. it is indeed not trivial: the machine had to calculate:

$\left(26 - 15\cdot\sqrt{3} \right) = (a - b)^{3} = a^3 -3a^2b + 3ab^2 - b^3$

We first see that b must contain the root of three as a factor. And 26 must contain the cube of a:

So we get: $26 = a^3 + 3ab^2$ Why this? the two terms on the right side may not contain a root. So the other two term with a root should be: $- 15\sqrt{3} = - 3a^2b - b^3$

After a little guess we get: $26 = 8 + 18 = 2^3 + 18$ and $18 = 3ab^2$ We set a = 2 and for b we get the root of 3 out of the last equation.

Now it was clear why the calculator isn't able to transform the terms, there is no simple recipe for this. Look at the more complicate expression for simplification:

$\left( 9416 - 4256\sqrt{5} \right)^{\frac{1}{3}} = ?$
(09-29-2020 09:22 PM)peacecalc Wrote: [ -> ]$- 15\sqrt{3} = - 3a^2b - b^3$

This is enough to solve it all Let b = c√3, we have -15√3 = -3a²(c√3) - (c√3)³

c * (a²+c²) = 5 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ⇒ c = 1, a = ±2

(±2 - √3)³ = (±8) - 3(4)(√3) + 3(±2)(3) - (3√3) = ±26 - 15√3

→ ³√(26 - 15√3) = 2 - √3

Quote:Look at the more complicate expression for simplification:
$\left( 9416 - 4256\sqrt{5} \right)^{\frac{1}{3}} = ?$

(a + b√5)³ = a³ + 3a²b (√5) + 3ab²(5) + b³ (5√5)

Collect radical free terms: a³ + 15 ab² = a * (a²+15b²) = 9416 = (2³)(11)(107)

(a²+15b²) > 0 ⇒ a > 0

15b² = 9416/a - a², thus RHS must be divisible by 3, and ends in 0 or 5

a=1: 9416/1-1² = 9415, not divisible by 3
a=2: 9416/2-2² = 4704, not end in 0 or 5
a=4: 9416/4-4² = 2338, not end in 0 or 5
a=8: 9416/8-8² = 1113, not end in 0 or 5

a=11: 9416/11-11² = 735 = 15*49 → |b| = 7

This work:

(11 ± 7√5)³ = 1331 + 3(121)(±7√5) + 3(11)(49*5) + (±343)(5√5) = 9416 ± 4256√5

→ ³√(9416 - 4256√5) = 11 - 7√5

---

Comment: both examples assumed the root has form a + b√k, all symbols integer.
Assumption might not hold, they might be rationals (and, more likely, irrationals).

For first example, had we solve the radical free terms, this also produce 26

(1 + (±5/3)(√3))³ = 1 + (3)(±5√3/3) + (3)(25/3) + (125/9)(√3) = 26 ± (170/9)(√3)
For general case, to solve for all a, b:

﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ³√(n + m√k) = a + b√k

(a + b√k)³ = a³ + 3a²b√k + 3ab²k + b³k√k

n = a³ + 3ab²k ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → n/a - a² = 3b²k
m = 3a²b + b³k ﻿ ﻿ ﻿ ﻿ ﻿ → 3m/b - 9a² = 3b²k

Equate the 2 to eliminate k, we have b = 3ma / (n + 8a³)

Substitute b to n = ... equation, we have a cubic equation, as function of A = a³

XCas> f(n, m, k) := horner([64, -48*n, 27*k*m^2-15*n^2, -n^3], A)

For first example, ³√(26 - 15√3):

XCas> expand(f(26,-15,3)) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 64*A^3 - 1248*A^2 + 8085*A - 17576
XCas> proot(ans()) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ →﻿ [5.75-1.125*i , 5.75+1.125*i , 8.0]

Above solved for A = a³. Each A produce 3 roots of a, thus a totaled 3×3 = 9 roots.
For this example, the real root is integer, a=2.

XCas> b := 3*m*a/(n + 8*a^3)
XCas> subst([a, b], [n,m,a] = [26,-15,2]) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → [2, -1]
Hello Albert,

thank you for responses!

You wrote:
Comment: both examples assumed the root has form a + b√k, all symbols integer.
Assumption might not hold, they might be rationals (and, more likely, irrationals).

That ist true, but when we don't make the assumption that the searched a's or b's are integers, the numbers of possibilites raised ad infinitum.

The problem I described raised for me to program the HP50g to solve the reduced equation: y^3 +py +q = 0. If the HP50g calculates in exact mode you get very wired expressions and even the roots are integers or simple roots, you cannot see this in the results (nor EVAL or SIMPLIFY helps). My work around is to call:

EVAL 10 FIX ->NUM ->PIQ (I tried this for a example which have the roots 5, -3 and root of 2. It is scary for getting exact values...
(09-30-2020 02:22 AM)Albert Chan Wrote: [ -> ](a + b√k)³ = a³ + 3a²b√k + 3ab²k + b³k√k

n = a³ + 3ab²k ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → n/a - a² = 3b²k
m = 3a²b + b³k ﻿ ﻿ ﻿ ﻿ ﻿ → 3m/b - 9a² = 3b²k

Equate the 2 to eliminate k, we have b = 3ma / (n + 8a³)

Instead of solving cubic, we can simply filter all (a,b), keeping only integers.
Since a divides n, just check divisors of n:

XCas> find_ab(n,m) := remove(x -> frac(x), map(divisors(n).*sign(n) , a->[a, 3*m*a/(n+8*a^3)]))
XCas> find_abk(n,m) := map(find_ab(n,m), x -> x + x * sqrt((n/x-x^2)/(3*x^2)))

XCas> find_abk(26,-15) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$[2 - \sqrt{3}]$$
XCas> find_abk(9416, -4256) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$[11 - 7 \sqrt{5}]$$

XCas> simplify((99+100*sqrt(101))^3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$300940299 + 103940300\sqrt{101}$$
XCas> find_abk(300940299, 103940300) ﻿ ﻿ ﻿ ﻿ ﻿ → $$[99 + 100 \sqrt{101}]$$
Here is an equivalent version, faster and shorter.

sign(b√k) = sign(b) = sign(m)

$$|b| \sqrt{k} = |b| \sqrt{{n/a-a^2 \over 3b^2}} = \sqrt{{n/a-a^2 \over 3}}$$

XCas> find_a(n,m) := remove(a -> remain(3*m, n/a+8*a*a), divisors(n) .* sign(n))
XCas> find_ar(n,m) := map(a -> a + sign(m)*sqrt((n/a-a*a)/3) , find_a(n,m))

XCas> find_ar(300940299, 103940300) ﻿ ﻿ ﻿ ﻿ ﻿ → $$[99 + 100 \sqrt{101}]$$
Hello Albert,

thank you for sharing your knowledge, it is very enlighting for me!.

I will try to program this with the help of the hp50g CAS, which is a subset of XCAS. If I succeed, I will share it, too.
Timings suggest it may be faster to solve the cubic, than factoring for divisors.
Again, we assumed n, m, k all integers.

Code:
find_cbrt(n,m,k):= {    local a;      // solved for A=a^3, then take cube root for a   a := fsolve(horner([64,-48*n,-15*n*n+27*m*m*k,-n^3],x), x=n);   a := exact(surd(a,3));   if (remain(n,a)) return surd(n+m*sqrt(k),3);   return a + sign(m) * sqrt((n/a-a*a)/3); }

XCas> find_cbrt(300940299, 103940300, 101) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$99 + 100 \sqrt{101}$$

If simplified form search failed, it returned the cube root as-is.

XCas> find_cbrt(300940299, 103940300, 102) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$\left(300940299 + 103940300\cdot \sqrt{102}\right)^{\frac{1}{3}}$$
Hello Albert,

I tried two way, the first was to probe every small a from 1 until third root of A. From your example a has to tested from 1 to 670, what a luck, that a = 99. My second approach is to calculate every divisor of A (with DIVI) as a possible candidate for small a. Then I'm calculate the difference A - third power of a. Negative results can be sorted out. The rest is divided by 3*k*a and be square rooted. If the result is a integer and not a algebraic expression would be a solution for smal b. and we get:

$\left( A\pm B\cdot\sqrt{k} \right)^{1/3} = a \pm b\cdot\sqrt{k}$.

For me the second approach makes mor sense, because the numbers of possibilities can be reduce very fast.

Code:
 « \-> E   « E 3 ^ DUP 'E' STO OBJ-> ->STR     IF "+" ==     THEN 1     ELSE -1     END SWAP DROP SWAP OBJ-> DROP2 DUP 2 ^ EVAL \-> A SG B RT K     « A DIVIS SORT \-> L       « A L 3 ^ -         «  \-> B           «             IF B 0 >             THEN B             END           »         » DOLIST          \-> L3         « L 1.0000 L3 SIZE SUB 'L' STO L3 L / 3 K * / \SQRT 'L3' STO            L3           «  \-> B             « B TYPE               IF 28.0000 ==               THEN B               ELSE 0               END             »           » DOLIST            'L3' STO { }                       1.0000 L3 SIZE           FOR I L3 I GET DUP             IF 0 ‹             THEN + I +             ELSE DROP             END           NEXT 'L3' STO                      IF L3 { } <>           THEN              1.0000 L3 SIZE             FOR I L3 I GET SG * L L3 I 1.0000 + GET GET SWAP RT * +             EVAL              2.0000 STEP           ELSE A B SG RT * * + 3 INV ^ EVAL           END         »       »     »   » »
Hello all,

an example for finding:

$\left( 1859814842094 - 59687820010\cdot\sqrt{415}\right)^{1/3} = 11589 - 145\cdot\sqrt{415}$

with a real HP50g in 27 sec. 1859814842094 has 96 divisors, after A - a^3 < 0 remain 25 divisors and the last one shows the integer 145 after dividing with 3*415*a and extracting root.
(10-04-2020 06:05 PM)peacecalc Wrote: [ -> ]I tried two way, the first was to probe every small a from 1 until third root of A. From your example a has to tested from 1 to 670, what a luck, that a = 99. My second approach is to calculate every divisor of A (with DIVI) as a possible candidate for small a. Then I'm calculate the difference A - third power of a. Negative results can be sorted out. The rest is divided by 3*k*a and be square rooted. If the result is a integer and not a algebraic expression would be a solution for smal b. and we get:

$\left( A\pm B\cdot\sqrt{k} \right)^{1/3} = a \pm b\cdot\sqrt{k}$

Based on above description ... not really.

You are solving ³√(A ± B√k) = a ± √r

r = b²k = (A-a³)/(3a) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // note: no B's !

(I might be wrong, your code may have used B somewhere)

Here is my lua code, that checked b actually get to B

Code:
cbrt = require'mathx'.cbrt function simp_cbrt1(A, B, k)         -- simplify cbrt(A + B * sqrt(k))     local N = abs(A)     for a = 1+(N-1)%3, cbrt(N), 3 do -- a mod 3 = A mod 3         if N % a ~= 0 then continue end         local r = (N/a - a*a)/3         if r % k ~= 0 then continue end         local b = B / (3*a*a+r)         if r == b*b*k then return (A<0 and -a or a), b, k end     end end

lua> A, B, k = 1859814842094, -59687820010, 415
lua> simp_cbrt1(A,B,k)
11589 ﻿ ﻿ ﻿ ﻿ ﻿ -145 ﻿ ﻿ ﻿ ﻿ ﻿ 415
lua> simp_cbrt1(A,B+1,k)﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -- nothing returned
lua>

Update2: fixed a bug if A is negative
Update3: inspired by simp_cbrt2(), calculate b without square root (i.e. not from r = b²k)
Hello Albert,

thank you for critics. Of course I forgot to implement the testing procedure being shure that I found a correct solution! The checking is absolut necessary.

Code:
 « \-> E   « E 3 ^ DUP 'E' STO OBJ-> ->STR     IF "+" ==     THEN 1     ELSE -1     END SWAP DROP SWAP OBJ-> DROP2 DUP 2 ^ EVAL \-> A SG B RT K     « A DIVIS SORT \-> L       « A L 3 ^ -         «  \-> B           «             IF B 0 >             THEN B             END           »         » DOLIST          \-> L3         « L 1.0000 L3 SIZE SUB 'L' STO L3 L / 3 K * / \SQRT 'L3' STO            L3           «  \-> B             « B TYPE               IF 28.0000 ==               THEN B               ELSE 0               END             »           » DOLIST            'L3' STO { }                       1.0000 L3 SIZE           FOR I L3 I GET DUP             IF 0 ‹             THEN + I +             ELSE DROP             END           NEXT 'L3' STO                      IF L3 { } <>           THEN              1.0000 L3 SIZE             FOR I               L3 I GET               L L3 I 1.0000 + GET GET               \-> a b               «                 B b 2 ^ K * a  2 ^ 3 * + b *                 IF ==                 THEN a b SG RT * * + EVAL                 ELSE E 3 INV ^ EVAL                  END              »                      2.0000 STEP           ELSE E 3 INV ^ EVAL           END         »       »     »   » »

Code is corrected, now the program tests if the found "a" leads to B as input.
(10-05-2020 11:36 AM)peacecalc Wrote: [ -> ]thank you for critics. Of course I forgot to implement the testing procedure being shure that I found a correct solution!
The checking is absolut necessary.

This checking is the reason my previous posts had different arguments requirement.

find_ar(n, m) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // note, no k's

It returned simplified a + b√k, but the user *must* check √k matches.

find_cbrt(n, m, k) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // we required k to build the cubic

It solved the cubic, checked a is integer, but assumed b is also integer.
Here is a proof that shows b is indeed an integer.

(09-30-2020 02:22 AM)Albert Chan Wrote: [ -> ]For general case, to solve for all a, b:

﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ³√(n + m√k) = a + b√k

(a + b√k)³ = a³ + 3a²b√k + 3ab²k + b³k√k

n = a³ + 3ab²k ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → n/a - a² = 3b²k
m = 3a²b + b³k ﻿ ﻿ ﻿ ﻿ ﻿ → 3m/b - 9a² = 3b²k

Since a|n, n/a - a² = 3b²k = integer

Assume b ≠ integer, but b = c/3, where c = integer.

m = 3a²b + b³k = a²c + c³k/27 = integer

Assuming √k is fully reduced, k can not have factor of 27 = 3³

(c³k/27 = integer) ⇒ (3|c) ⇒ (b = c/3 = integer)

Assumption were wrong, b must be integer. ﻿ ﻿ ﻿ ﻿ ﻿ QED
Hello Albert,

The HP50g reduces every squareroot if it contains a single integer fully. So the input and the output has to had the same squareroot k. My program takes that from the input and doesn't change it. So what do I overlook?

For calculating A and B I use k as a given value, hmm is that wrong? Are there possibly other solutions with another k as given? But when this would be possibly then we have a lot of more solutinon, too much...?
(10-06-2020 05:25 AM)peacecalc Wrote: [ -> ]For calculating A and B I use k as a given value, hmm is that wrong? Are there possibly other solutions with another k as given? But when this would be possibly then we have a lot of more solutinon, too much...?

It is possible for the same A,B, but different k's, we can simplify the cube roots.
(Actually, I am not so sure. Seems pretty rare for this to happen ...)

I know that for same A,k, but different B's, it is possible (at least, if we included imaginary numbers)

Below searched for all A, such that ³√(A ± √R) = a ± √r, ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ where R = B²k, r = b²k

r = b²k = (A/a - a*a)/3 = (x-y)/3, where x=A/a, y=a*a
B = b * (3a² + b²k)
R = B²k = r * (3a² + r)^2 = (x-y)/3 * (3y + (x-y)/3)² = (x-y)/27*(x+8*y)²

XCas> getR(x,y) := (x+8*y)^2*(x-y)/27
XCas> find_all_a(A) := remove(a-> denom(getR(A/a,a*a))!=1, divisors(A).*sign(A))
XCas> find_all_abk(A) := map(a -> a + sqrt((A/a-a*a)/3), find_all_a(A))

For A=90, k=7, we have 2 cube roots able to simplify, B = 34, 101i
Note: because of "±", we also have the "-" solutions.

XCas> s := find_all_abk(90)﻿ ﻿ → $$[\;3+\sqrt7,\quad\quad 6+i\sqrt7,\quad\;\; 15+i\sqrt{73}, \quad\quad 30+i\sqrt{299}]$$
XCas> simplify(s .^ 3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$[90+34\sqrt7,\; 90+101i\sqrt7,\; 90+602i\sqrt{73}, \;90+2401i \sqrt{299}]$$
I discovered a necessary condition to simplify cubic: ³√(A ± √R) = a ± √r If we found a, we can get r = (A/a-a*a)/3

Simplify equations from find_cbrt(): 4a³ - A = (3a) ³√(A²-R)

This implied ³√(A²-R) must be integer !

Example, ³√(9416 - 4256√5)

A² - R = 9416² - 4256² * 5 = -1906624
c = ³√(A² - R) = -124 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // cube root may be simplified

4a³ = A + 3*c*a
a = ³√(A/4 + (3c/4)*a)

With my Casio FX 115MS:

5 =
³√(2354 - 93 Ans
= ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 12.36167494
= ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 10.63945257
= ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 11.09160665
= ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 10.97648027 // locking to 11, try it
11 = ↑
= ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 11 ﻿ ﻿ ﻿ ﻿ ﻿ // = a
9416 / Ans - Ans Ans
= ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 735 ﻿ ﻿ ﻿ // = 3r
√(Ans / 3 / 5
= ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 7 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // = b

Add the sign, ³√(9416 - 4256√5) = 11 - 7√5
(10-06-2020 12:06 PM)Albert Chan Wrote: [ -> ]Simplify equations from find_cbrt(): 4a³ - A = (3a) ³√(A²-R)

This implied ³√(A²-R) must be integer !

There is a much more elegant way to get above relation !

³√(A ± √R) = a ± √r
³√(A +√R) ³√(A -√R) = (a + √r) (a - √r)
³√(A² - R) = a² - r = a² - (A/a - a*a)/3 = (4a² - A/a)/3

Because RHS is integer, LHS must too.

We have a|A, and 3|(4a² - A/a), which simplified to a≡A (mod 3)
I've updated my lua code now to check a in steps of 3 Update: Above assumed a is integer; it might be halves.
Including halves, we have 2a ≡ 2a*(A/a) ≡ 2A (mod 3)
(10-06-2020 04:13 PM)Albert Chan Wrote: [ -> ]³√(A +√R) ³√(A -√R) = (a + √r) (a - √r)
³√(A² - R) = a² - r = a² - (A/a - a*a)/3 = (4a² - A/a)/3

We can solve above cubics directly.
However, with limited precision, it is hard to get the correct LHS (is it really an integer ?)
One trick is assume that it is (with error under ± 1/2), i.e. we use c = round(³√(A² - R))

We solve the cubics, 4a³ = (3c)*a + A, then confirm (a, r) will, in reverse, produce (A, R)
(Actually, it calculated 3r both ways, using either A or B; both must match)

Below uses Kahan's Solve a Real Cubic Equation method

Code:
cbrt = require'mathx'.cbrt function cubic_guess(c,d)   -- guess for x^3 = c*x + d     if c <= 0 then return cbrt(d) end     c = 1.324718 * max(cbrt(abs(d)), sqrt(c))     return d < 0 and -c or c end     function simp_cbrt2(A,B,k)  -- simplify cbrt(A + B * sqrt(k))     local H = 0x3p50        -- for rounding to halves     local c = cbrt(A*A - B*B*k) + 2*H-2*H     local a = cubic_guess(0.75*c, 0.25*A)     repeat                  -- Newton's method         local a0, x = a, 4*a*a         a = a - (a*(x-3*c)-A) / (3*(x-c))     until a+H/2 == a0+H/2     a = a+H-H     local b = B/(4*a*a-c)   -- B = b*(4*a*a-(a*a-b*b*k))     if 3*b*b*k == A/a-a*a then return a,b,k end end

lua> simp_cbrt2(1859814842094, -59687820010, 415)
11589 ﻿ ﻿ -145 ﻿ ﻿ ﻿ 415
lua> simp_cbrt2(300940299,103940300,101)
99 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 100 ﻿ ﻿ ﻿ ﻿ ﻿ 101
lua> simp_cbrt2(9416,-4256,5)
11 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -7 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 5
lua> simp_cbrt2(26,-15,3)
2 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -1 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 3
lua> simp_cbrt2(26,-15+1,3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -- no returns, as expected
lua>

For first example, this version run 55X speed of previous simp_cbrt1()
To show why, lets add a debug line, "print(a)", right below keyword "repeat"

lua> simp_cbrt2(1859814842094, -59687820010, 415) -- debug version
12856.226745378446
11738.132585006717
11591.443457378273
11589.000672078266
11589 ﻿ ﻿ -145 ﻿ ﻿ ﻿ 415
Unfortunately, our divisors based cube root simplify routines were flawed !

(10-01-2020 05:22 PM)Albert Chan Wrote: [ -> ]
Code:
find_cbrt(n,m,k):= {    local a;      // solved for A=a^3, then take cube root for a   a := fsolve(horner([64,-48*n,-15*n*n+27*m*m*k,-n^3],x), x=n);   a := exact(surd(a,3));   if (remain(n,a)) return surd(n+m*sqrt(k),3);   return a + sign(m) * sqrt((n/a-a*a)/3); }

Real rational root for a is a divisor of n, *but* possibly divided by 4.
Above XCas version is correct. It does not assume a is integer.
Updated previous post simp_cbrt2() code to handle "quarter" bug.

XCas> find_cbrt(81,30,-3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$\frac{9}{2} + \frac{i}{2}\cdot\sqrt{3}$$
lua> simp_cbrt2(81,30,-3) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 4.5﻿ ﻿ ﻿ 0.5﻿ ﻿ ﻿ -3

Update: roots impossible as quarters, but possible as halves.
Updated previous post simp_cbrt2() code to handle "halves" bug.

---

I added a simple third version, going for converged (a, r), rational or not.
This may be convenient if we do not have √R factored to B√k
Code:
function simp_cbrt3(A,R)    -- simplify cbrt(A +/- sqrt(R))     local c = cbrt(A*A - R)     local a = cubic_guess(0.75*c, 0.25*A)     repeat                  -- Newton's method         local x = 4*a*a-c         local eps = (a*(x-2*c)-A) / (3*x)         a = a - eps     until a == a + eps*eps     return a, R/(4*a*a-c)^2 end

lua> simp_cbrt3(10, 108) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ -- ³√(10 ± √108) = 1 ± √3
1 ﻿ ﻿ ﻿ ﻿ 3

lua> A, R = 123, 456
lua> a, r = simp_cbrt3(A, R) ﻿ ﻿ ﻿ ﻿ ﻿ -- not possible to simplify, (A²-R) not cubes
lua> a, r
4.956345549114864 ﻿ ﻿ ﻿ ﻿ 0.0837701441234425
lua> cbrt(A + sqrt(R)) , a + sqrt(r)
5.245776273336468 ﻿ ﻿ ﻿ ﻿ 5.245776273336468
(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{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)
Pages: 1 2
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :