HP50g simplifing a root

[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