3

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

?

Carl Christian
  • 13,666
Sam Ginrich
  • 161

2 Answers2

6

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)
Lutz Lehmann
  • 131,652
  • Don't you think that the problem is the square root ? I really would appreciate your opinion on this point. Cheers and (+1). – Claude Leibovici Dec 07 '21 at 02:07
  • 1
    Yes, with the square root you lose a bit or two. But the real cancellation is in the outer difference, the whole cancellation issue is more of the type $1-(1-x)$. In the transformed expression in the denominator, in the sum with $1$, these bits remain outside the mantissa range, so nothing is lost in the limits of the floating point type. – Lutz Lehmann Dec 07 '21 at 07:47
  • Thanks for answering. Having started with computers 60+ years ago, I am always interested by diffrent points of view about this kind of topics. – Claude Leibovici Dec 07 '21 at 07:53
  • Regard sqrt(1-y^4) = 1- y^4/2 + O(y^6), no loss! – Sam Ginrich Dec 07 '21 at 13:58
  • (+1) Given that the $y_n$ are known to be in $[0,1]$, it seems to me (and some quick experiments appear to confirm), that accuracy could be slightly improved by making the last step T = y4 / (((w4 * w2 + w2) + w4) + 1). – njuffa Sep 06 '22 at 11:52
3

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

Claude Leibovici
  • 289,558
  • This comes brilliant in the asymptotic case, when O(y^20) is under the carpet. – Sam Ginrich Dec 06 '21 at 19:10
  • @SamGinrich. I wanted to avoid the square root – Claude Leibovici Dec 07 '21 at 02:08
  • The Borwein Algorithm requires full precision from the beginning and stability is crucial. – Sam Ginrich Dec 07 '21 at 12:51
  • sqrt complexity: Guess, for "really large" numbers division by iterated multiplications is cheaper than other methods, so division is in the same quadratic convergence class as sqrt, where sqrt starts with an additional inversion of the radicand. According to the number of required mulitplications, the effort ratio sqrt:division is roughly (2+3):2. Hence, Lutz Lehmann's approach has effort 12 and your Padé Approximation effort 2. ´I'll combine both solutions. – Sam Ginrich Dec 07 '21 at 18:32