HP50g simplifing a root

[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