Cancellation Problem in Borwein's quartic Algorithm for computing $\pi$

Question

In the center of Borwein's Algorithm from 1985 occurs a term $$T_n = 1-\sqrt[4]{1-y_n^4}$$ with $y_n \to 0_+$. For small $y_n$ you need 4 times the precision of result $T_n$ due to cancellation. With the decomposition $$\sqrt[4]{1-y_n^4} = \sqrt[4]{1-y_n^2}\sqrt[4]{1+y_n^2}$$ you still loose half of the precision! A further complex decomposition of the latter radicand would remain quadratic.

On the other hand, regarding the Taylor series, $T_n$ is quite flat, $O(y_n^4)$, and there is a fast linear convergent evaluation without loss.

Question: Is there a numerically stable evaluation algorithm for $T_n$ with at least quadratic convergence, e.g. similar to what we have for the quartic root

$$r_{\infty}=\sqrt[4]{x}: \ r_{n+1} = r_n \cdot \{\frac{5}{4} - \frac{r_n^4}{4x}\} $$

?

Answer 1

If you want more accurate than Taylor expansions, think about Padé approximants.

The simplest would be

$$1-\sqrt{1-y^4}\sim\frac{2 y^4}{8-3 y^4}+O(y^{12})$$ To give an idea $$\Phi=\int_0^{\frac 12} \Bigg[\Big[1-\sqrt{1-y^4} \Big]-\frac{2 y^4}{8-3 y^4} \Bigg]^2\,dy=5.22\times 10^{-12}$$ We can do much better, as for example $$1-\sqrt{1-y^4}\sim \frac{8 y^4 \left(2-y^4\right) }{64-56y^4+7y^8 }+O(y^{20})$$ $$\Phi=\int_0^{\frac 12} \Bigg[\Big[1-\sqrt{1-y^4} \Big]-\frac{8 y^4 \left(2-y^4\right) }{64-56y^4+7y^8 } \Bigg]^2\,dy=2.37\times 10^{-20}$$

Answer 2

Use the binomial formula $$ A^4-B^4=(A-B)(A+B)(A^2+B^2) $$ with $A=1$ and $B=\sqrt[4]{1-y_n^4}$ to get a formula without cancellation, $$ 1-\sqrt[4]{1-y^4}=\frac{y^4}{(1+\sqrt[4]{1-y^4})(1+\sqrt{1-y^4})}. $$ The parts in the resulting formula can be computed in the following steps as

y2 = y*y
y4 = y2*y2
w2 = sqrt(1-y4)
w4 = sqrt(w2)
T = y4/(1+w4)/(1+w2)