Question for Trig Gurus

[Albert Chan]

For 5+ digits accuracy, code is simple, using only 2 square roots.

Code:

function atand(x) -- = deg(atan(x)), |x| <= 1, 5+ digits accuracy
    local g = sqrt(1+x*x)
    local a = (1+g)/2
    a = 1 + 12*a + 32*sqrt(a*g)
    return x/a * 8100/pi
end

Assume we have signed zero, atand(x) = sign(x)*90 - atand(1/x)
--> maximum error (both absolute and relative) at x = ±1

lua> r3 = sqrt(3)
lua> atand(1/(2+r3)), atand(1/r3), atand(1)
15.000000161153318      30.000021200229554      45.00037955691671

Update: for |x| = 1/√3 .. √3, we can map it within ±1/√3, getting 6+ digits accuracy

atand(x) = sign(x)*45 - atand((1-x*x)/(2*x)) / 2

---

Code based on sequence a_k = (x/2^k) / tan(atan(x)/2^k)
see An algorithm for computing Logarithms and ArcTangents, by B. C. Carlson

a0 = x/x = 1, g0 = sqrt(1+x*x),
a1 = (x/2) / tan(atan(x)/2) = (a0+g0)/2, g1 = sqrt(a1*g0)
a2 = (x/4) / tan(atan(x)/4) = (a1+g1)/2

We use Richardson extrapolation to extrapolate for a_∞ = x / atan(x)
Since we don't care intermediate values, we can apply weight directly.
We just need to know what weight to use ...

CAS> [1]
CAS> 4^len(Ans)*Ans - extend([0],Ans)
[4, -1]
[64, -20, 1]
[4096, -1344, 84, -1]
[1048576, -348160, 22848, -340, 1]
...

45*a_∞ ≈ (a0 - 20*a1 + 64*a2) = (a0 - 20*a1 + 32*(a1+g1)) = (a0 + 12*a1 + 32*g1)