HP50g simplifing a root

[peacecalc]

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}} = ? \]

[Albert Chan]

(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)

[Albert Chan]

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]

[peacecalc]

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...

[Albert Chan]

(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[1]), 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[0] + x[1] * sqrt((n/x[0]-x[0]^2)/(3*x[1]^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}]\)

[Albert Chan]

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}]\)

[peacecalc]

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.

[Albert Chan]

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}}\)

[peacecalc]

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
        »
      »
    »
  »
»

[peacecalc]

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.

[Albert Chan]

(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
see https://www.hpmuseum.org/forum/thread-15...#pid136957

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>

Update: added a ≡ A (mod 3), see https://www.hpmuseum.org/forum/thread-15...#pid137252
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)

[peacecalc]

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.

[Albert Chan]

(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

[peacecalc]

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...?

[Albert Chan]

(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}]\)

[Albert Chan]

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

[Albert Chan]

(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)

[Albert Chan]

(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

[Albert Chan]

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.
see https://www.hpmuseum.org/forum/thread-15...#pid137252

---

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

[Albert Chan]

(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)