Trig algorithms on HP RPN?

[Albert Chan]

(11-30-2024 03:49 PM)Thomas Klemm Wrote:  

Code:

0-0727 0530  load constant 5        reg['C'] = 0x00785398163500

Should last digit for pi/4 a 4?

lua> pi/4
0.7853981633974483

(11-30-2024 03:49 PM)Thomas Klemm Wrote:  Honestly, I don't know why that wasn't used instead.

Just a guess, if angles are in radian, we can just use tan code to get atan

tan sum formula: tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
Flip for atan:          atan(x) = atan( y + (x-y) )

atan(x) = atan(y) + atan(z),     where z = (x-y)/(1+x*y)

if y = tan(x), (x-y) ≈ x^3/3 --> z shrink fast, cubic convergence!

But we can do better!

y can be set to anything. As long as y is close to x, z will shrink fast.
Below, y = tan(x2 = x/(1+x*x/3)), to get convergence rate of O(x^5)

Code:

function myatan(x, s)
    if abs(x)>1 then return myatan(-1/x, copysign(pi/2,x)) end
    s = s or 0
    print(s .. " + atan( " .. x .. " )")
    if 3+x*x==3 then return s + x end   -- atan(x) ≈ x
    local x2 = x / (1+x*x/3)
    local y = tan(x2)           -- y is close to x
    local z = (x-y) / (1+x*y)   -- atan(x) = atan(y) + atan(z)
    return myatan(z, s + x2)    -- atan(y) = x2
end

lua> myatan(0.1)
0 + atan( 0.1 )
0.09966777408637874 + atan( 8.784047832967151e-07 )
0.09966865249116204 + atan( -1.0587911840670585e-22 )
0.09966865249116204

lua> myatan(1)
0 + atan( 1 )
0.75 + atan( 0.03541295579818369 )
0.7853981584542247 + atan( 4.943223638603625e-09 )
0.7853981633974483