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

Solution 1:

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. 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)

Solution 2:

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}$$