arctangent approximation with explicit error bound

Question

I need to be able to generate rational approximations of arctangent evaluated at rational numbers to arbitrarily small error margins. In particular, I need to be able to determine when I've done enough iterations of the approximation algorithm to guarantee that the result is within the error margin. I would also like a method that doesn't involve any operations or values that themselves would have to be approximated, like root extractions, to spare myself the complexity of having to figure out how accurate each of those intermediary approximations would have to be in order for the approximation of the arctangent value to be valid to the desired accuracy; this probably means arithmetic operations on rational numbers only.

The Taylor series representation meets both of these requirements - since it's an alternating series, I'd know I've reached the desired accuracy once the magnitude of the last term added was smaller than the error margin - but I'd prefer a different method because the Taylor series converges so slowly. Are there faster methods available that meet my requirements? (6) from this paper would be suitably easy to calculate, but I don't know how to derive an error bound for it.

Answer 1

Using $$\tan^{-1}(x)=\sum_{n=0}^\infty (-1)^n\frac{x^{2n+1}}{2n+1}$$ you need to find $$\frac{x^{2p+3}}{2p+3} \leq \epsilon$$ in order that the summation up to $p$ meets your requirement.

Excluding the case $x=1$, the solution is given by $$p \geq-\frac{3}{2}-\frac{1}{2 \log (x)}W\left(-\frac{\log (x)}{\epsilon }\right)$$ where appears Lambert function.

Since the argument is large (positive for $x < 1$), you can use the expansion $$W(t)=L_1-L_2+\frac{L_2}{L_1}+\frac{L_2(L_2-2)}{2L_1^2}+\frac{L_2(6-9L_2+2L_2^2)}{6L_1^3}+\cdots$$ where $L_1=\log(t)$ and $L_2=\log(L_1)$.

If you set $\epsilon=10^k$, you should notice that $p$ is almost linear with $k$ .

Answer 2

Assume you did the usual reductions so that $x=\frac ab\in (0,1)$. Then as you know, the angle is the argument of the complex numbers and is the same for all points on the ray through $(b,a)$ $$ \arctan\frac ab=\arg(1+i\frac ab)=\arg(b+ia) $$ Now the aim is to reduce this point to somewhere close to the real axis. For that purpose one might try powers of $(1+\frac1{bm})^k$ where $k$ is close to $ma$, for instance for $x=\frac34$ and $m=2$ we get

k      z = (8+i)^k        4*Im(z)/Re(z)    (4+3i)/z (scaled)

4     3713 +    2016i    2.17182870994      836 +   123j 
5    27688 +   19841i    2.8663681017      6811 +   148j 
6   201663 +  186416i    3.69757466665    54636 -  5627j 
7  1426888 + 1692991i    4.74596744804   431461 - 99652j 

From that we read off the formulas \begin{align} \arctan \frac34 &= 5\arctan\frac18+\arctan\frac{148}{6811}\\[.5em] \arctan \frac34 &= 6\arctan\frac18-\arctan\frac{5627}{54636} \end{align} where the arcus tangent series on the right converge faster than the original on the left.

In practice one would chose a larger $m$ and do binary search or similar to rapidly find the correct power of $k$.