HP50g simplifing a root
10-07-2020, 06:12 PM (This post was last modified: 10-10-2020 05:25 PM by Albert Chan.)
Post: #18
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: HP50g simplifing a root
(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
 « Next Oldest | Next Newest »

 Messages In This Thread HP50g simplifing a root - peacecalc - 09-29-2020, 09:22 PM RE: HP50g simplifing a root - Albert Chan - 09-29-2020, 11:47 PM RE: HP50g simplifing a root - Albert Chan - 09-30-2020, 02:22 AM RE: HP50g simplifing a root - Albert Chan - 09-30-2020, 10:50 PM RE: HP50g simplifing a root - Albert Chan - 10-01-2020, 07:31 AM RE: HP50g simplifing a root - peacecalc - 09-30-2020, 05:33 AM RE: HP50g simplifing a root - peacecalc - 10-01-2020, 02:20 PM RE: HP50g simplifing a root - Albert Chan - 10-01-2020, 05:22 PM RE: HP50g simplifing a root - peacecalc - 10-04-2020, 06:05 PM RE: HP50g simplifing a root - Albert Chan - 10-04-2020, 11:48 PM RE: HP50g simplifing a root - peacecalc - 10-04-2020, 07:36 PM RE: HP50g simplifing a root - peacecalc - 10-05-2020, 11:36 AM RE: HP50g simplifing a root - Albert Chan - 10-05-2020, 05:01 PM RE: HP50g simplifing a root - peacecalc - 10-06-2020, 05:25 AM RE: HP50g simplifing a root - Albert Chan - 10-06-2020, 09:40 AM RE: HP50g simplifing a root - Albert Chan - 10-06-2020, 12:06 PM RE: HP50g simplifing a root - Albert Chan - 10-06-2020, 04:13 PM RE: HP50g simplifing a root - Albert Chan - 10-07-2020 06:12 PM RE: HP50g simplifing a root - Albert Chan - 10-09-2020, 12:20 AM RE: HP50g simplifing a root - Albert Chan - 10-09-2020, 02:31 PM RE: HP50g simplifing a root - Albert Chan - 10-11-2020, 06:28 PM RE: HP50g simplifing a root - Albert Chan - 10-12-2020, 03:17 AM RE: HP50g simplifing a root - Albert Chan - 10-24-2020, 02:19 PM RE: HP50g simplifing a root - Albert Chan - 10-12-2020, 10:54 PM RE: HP50g simplifing a root - CMarangon - 10-12-2020, 11:45 PM RE: HP50g simplifing a root - grsbanks - 10-13-2020, 06:46 AM RE: HP50g simplifing a root - Albert Chan - 10-09-2020, 05:21 PM RE: HP50g simplifing a root - Albert Chan - 10-10-2020, 03:58 PM RE: HP50g simplifing a root - Albert Chan - 10-10-2020, 04:49 PM RE: HP50g simplifing a root - peacecalc - 10-12-2020, 08:49 PM RE: HP50g simplifing a root - peacecalc - 10-13-2020, 06:30 AM RE: HP50g simplifing a root - peacecalc - 10-13-2020, 06:36 AM

User(s) browsing this thread: 1 Guest(s)