HP 17B Solver - another ARCTAN(Y/X) approximation

[Albert Chan]

The approximation formula were good, giving 6+ significant digits.
But, we can do better, by more argument reduction.

Let y = tan(θ):

tan(θ - atan(x)) = (y-x)/(1+y*x)         // tan(a ± b) = (tan(a) ± tan(b)) / (1 ∓ tan(a)*tan(b))
θ - atan(x) = atan((y-x)/(1+y*x))       // assume |θ - atan(x)| ≤ pi/2, so atan(tan()) round-trip

atan(x) = θ + atan((x-y)/(1+x*y))

We implicitly applied this when |x| > 1, with θ = pi/2:

atan(x) = θ + atan(x/y-1)/(1/y+x) = pi/2 + atan((0-1)/(0+x)) = pi/2 - atan(1/x)

We could do 1 more, with θ = pi/4:

atan(x) = pi/4 + atan((x-1)/(x+1))

However, this is more messy with unknown sign of x (we may need θ = -pi/4 instead)

I was wondering ... is there a half-"angle" atan formula ? atan(x) = 2*atan(h) + atan(ε)
This is the required taylor series for h, to keep ε small.

XCas> series(tan(atan(x)/2),x,0,12)

x/2 - x^3/8 + x^5/16 - 5*x^7/128 + 7*x^9/256 - 21*x^11/1024 + x^13*order_size(x)

Using Pade approximation, I found one that atan(ε) can be approximated by ε

XCas> h := x*(1/6 + 1/(x^2+4) + 1/(9*x^2+12))
XCas> series(h,x,0,12)

x/2 - x^3/8 + x^5/16 - 5*x^7/128 + 7*x^9/256 - 41*x^11/2048 + x^13*order_size(x)

This is ε, by applying sum-of-tangents formula twice:

XCas> corr(x,y) := factor(((x-y)/(1+x*y))
XCas> corr(corr(x,h),h)       // = ε

-x^11 / (11*x^10 + 220*x^8 + 1232*x^6 + 2816*x^4 + 2816*x^2 + 1024)

Example, with x = 1 (where max-rel-error occurs), it will split as:

atan(x) = 2*atan(h) + atan(ε)
atan(1) = 2*atan(29/70) + atan(-1/8119)

atan(h) - estimated_atan(h) ≈ 3.66e-12
atan(ε) - ε ≈ -ε^3/3 ≈ 6.23e-13, small enough to ignore.

This compared to without splitting:

atan(1) - estimated_atan(1) = pi/4 - 45824/58345 ≈ 9.6e-7

Code:

function my_atan2(y, x)
    local z, offset, p = y/x, 0, (y >= 0 and pi or -pi)
    if abs(z)>1 then z, offset = -x/y, p/2 end
    if x < 0 then offset = p - offset end -- Q2, Q3
    local s = z*z
    x = z * (1/6 + 1/(s+4) + 1/(9*s+12))  -- x/z <= 1/2
    y = (z-x)/(1+z*x)   -- atan(z) = 2*atan(x) + atan(y)
    y = (y-x)/(1+y*x)   -- atan(y) approximated by y
    s = x*x             -- atan(x) approximated by below
    x = x * (((15159/35*s + 4213)*s + 9867)*s + 6435) /
      ((((35*s + 1260)*s + 6930)*s + 12012)*s + 6435)
    return 2*x + y + offset
end