Let's say that we have the function

$$f(x) = 1- \frac{\sqrt{1 + x^2}}{1 + x^2/2}$$

for small $x$.

What I am asking is the following:

I am going to solve this function numerically for $10^{-10}<x<1$. To have the best precision, should I use the function $f(x)$ as it is and make the computer solve it or should I expand the function using Taylor series? In other words:

Since computers use certain arithmetics (64-bit double-precision etc.) it can only contain some of the digits and there will always be an error (loss of precision). What I am trying to understand is can avoid it (not entirely of course) using Taylor series. if yes; To what order should I expand it so that computer can store that kind of numbers?

  • $\begingroup$ What do you mean by "solve this function numerically"? $\endgroup$ – nicoguaro Sep 19 '19 at 17:25
  • $\begingroup$ I mean I am going to plot the function for $10^{-10}<x<1$ using let's say Matlab. $\endgroup$ – Ekrem Sep 19 '19 at 17:30
  • $\begingroup$ I don't see why evaluating that function would be troublesome for that interval. $\endgroup$ – nicoguaro Sep 19 '19 at 18:06
  • $\begingroup$ It's not a specific interval. I just wanted to say that x is small and I need numbers with less error. Since computers use certain arithmetics (64-bit double-precision etc.) it can only contain some of the digits and there will always be an error (loss of precision). What I am trying to understand is can avoid it (not entirely of course) using Taylor series. $\endgroup$ – Ekrem Sep 19 '19 at 18:53
  • $\begingroup$ That's not what you ask in your OP, I think that you should correct it accordingly. Furthermore, numerical software is written with this kind of things in mind, I doubt that you can just get better results using a naive approach. $\endgroup$ – nicoguaro Sep 19 '19 at 18:58

If we want to approximate $$f(x) := 1 - \frac{\sqrt{1+x^{2}}}{1+\frac{1}{2}x^2}$$ on $[10^{-10}, 1]$ as accurately as possible, and in particular with an accuracy close to the machine precision associated with a given IEEE-754 floating-point format, neither naive computation nor a Taylor series expansion seem suitable. There is an issue with subtractive cancellation that is particular severe for small arguments. We can certainly expand $f(x)$ around $0$, e.g. with the help of Wolfram Alpha: $$\frac{x^{4}}{8}-\frac{x^{6}}{8}+\frac{13x^{8}}{128}-\frac{5x^{10}}{64}+\frac{61x^{12}}{1024}+\mathcal{O}(x^{13})$$

While this would allow us to accurately compute $f(x)$ for arguments very small in magnitude, it is not practical for $x$ near $1$. And since $f(1) \approx 0.0572$, it is easy to see that we encounter the subtractive cancellation issue even there when using the naive approach. Sure enough, a quick experiment shows that the naive computation leads to errors in excess of $256$ ulps on $[\frac{1}{2}, 1]$, meaning the least significant $8$ bits of results for that interval cannot be relied upon.

#include <math.h> // for sqrt()

double func_naive (double x)
    return 1.0 - sqrt (1.0 + x * x) / (1.0 + x * x / 2);

I briefly examined whether we can somehow use the function sqrt1pm1, which accurately computes $\sqrt{x + 1}-1$, to minimize the impact from subtractive cancellation, but I couldn't find a way to do this.

Intuition indicates that if we want to generate results small in magnitude from arguments large in magnitude in a numerically stable fashion, we could use computation based on division. There are various ways to do this. Given that here we are dealing with fixed-precision IEEE-754 floating-point formats and a specific, fairly narrow interval, a rational minimax approximation seems to be the best approach. A rational approximation is basically the ratio of two polynomials, i.e. $r(x) = p(x) / q(x)$. A numerically robust way of evaluating each polynomial is Horner's scheme.

Software like Mathematica or Maple can be used to find such approximations, or one could build one's own tool based on the Remez algorithm. I chose the latter approach here. On platforms with hardware and software support for fused multiply-add (FMA) we can use it to achieve a slightly more accurate computation.

Based on testing with 100M random arguments in $[10^{-10}, 1]$ with IEEE-754 binary64 (double precision) computation, maximum errors around $5$ ulps are observed, which seems reasonable for this quickly constructed code. With careful tweaking of the coefficients the maximum error could likely be reduced a bit further.

#include <math.h> // for fma()

#define FMAD(a,b,c) fma(a,b,c)
#define FMAD(a,b,c) (a*b+c)
#endif // HAVE_FMA

// Accurately compute 1-sqrt(1+x*x)/(1+x*x/2) on [1e-10, 1]
// HAVE_FMA=1: max ulp=4.55664 @ x=2.2644149793184315e-6; max rel err=6.36228e-16 @ x=3.5041674402633207e-4
// HAVE_FMA=0: max ulp=5.62062 @ x=6.4687640650198819e-1; max rel err=6.82920e-16 @ x=8.2367484833186755e-1
double func_rat_approx (double x)
    double p, q;
    double x2 = x * x;

    p =              - 3.2490268983475442E-5;
    p = FMAD (p, x2, + 2.2558281208246997E-3);
    p = FMAD (p, x2, + 7.9333170773742845E-1);
    p = FMAD (p, x2, + 1.1376929119278506E+1);
    p = FMAD (p, x2, + 4.8037227288592518E+1);
    p = FMAD (p, x2, + 7.1940878446934676E+1);
    p = FMAD (p, x2, + 3.4917444046359641E+1);

    q =          x2  + 2.2273111565876125E+1;
    q = FMAD (q, x2, + 1.6832675743004904E+2);
    q = FMAD (q, x2, + 5.8322456893291758E+2);
    q = FMAD (q, x2, + 1.0122010119537713E+3);
    q = FMAD (q, x2, + 8.5486657994635436E+2);
    q = FMAD (q, x2, + 2.7933955237087713E+2);

    return (p / q) * x2 * x2;
| cite | improve this answer | |
  • $\begingroup$ Thanks a lot for the explanation. Could you briefly explain the code? I didn't quite understand it. $\endgroup$ – Ekrem Sep 19 '19 at 23:25
  • 1
    $\begingroup$ @Ekrem: I'd be happy to improve the answer if you could you briefly explain what part of the code is unclear. A rational approximation is the ratio of two polynomials: $r(x) =p(x)/q(x)$. Here, both polynomials are evaluated using the Horner scheme. $\endgroup$ – njuffa Sep 19 '19 at 23:29
  • $\begingroup$ I don't know much c++. In fact, I know coding only at a basic level. I searched about this fma() function but I couldn't understand why did you use it. I kind of understood the math behind the Horner scheme but I don't see how you define $p$ and $q$ i.e. where those red numbers come from? Also why the function returns $(p/q) * x^{4}$? Sorry for having so much question. Thank you so much for your time. $\endgroup$ – Ekrem Sep 20 '19 at 20:51
  • 1
    $\begingroup$ @Ekrem FWIW, the code I posted is C, not C++. The "red numbers" are the coefficients of the polynomials, generated by my proprietary program that generates minimax approximations using the Remez algorithm. You can generate such approximations with widely-used commercial tools like Maple or Mathematica (which I cannot afford and therefore have never used). As to why the approximation takes the form $(p/q)*x^{4}$, see the Taylor series expansion further up in the answer. FMA provides two operations with a single rounding and gives some protection against subtractive cancellation; see Wikipedia $\endgroup$ – njuffa Sep 20 '19 at 21:55
  • $\begingroup$ Thank you. This helps a lot. $\endgroup$ – Ekrem Sep 20 '19 at 22:17

You can easily approximate $f(x)$ up to fourth order accuracy by using Taylor expansion as:

$$f(x) = 1 - \frac{\sqrt{1+x^{2}}}{1+\frac{x^{2}}{2}}$$

Taylor expansions:

$$\sqrt{1+x^{2}} = 1 + \frac{x^{2}}{2} - \frac{x^{4}}{8} + O(x^{5})$$

$$\frac{1}{1+\frac{x^{2}}{2}} = 1 - \frac{x^{2}}{2} + \frac{x^{4}}{4} + O(x^{6})$$


$$f(x) = 1 - \frac{\sqrt{1+x^{2}}}{1+\frac{x^{2}}{2}} = 1 - \frac{1 + \frac{x^{2}}{2} - \frac{x^{4}}{8}}{1+\frac{x^{2}}{2}} = 1 - \frac{1}{1+\frac{x^{2}}{2}} - \frac{x^{2}}{2} (\frac{1}{1+\frac{x^{2}}{2}}) + \frac{x^{4}}{8}(\frac{1}{1+\frac{x^{2}}{2}})$$


$$f(x) = 1 - 1 + \frac{x^{2}}{2} - \frac{x^{4}}{4} - \frac{x^{2}}{2}(1 - \frac{x^{2}}{2}+\frac{x^{4}}{4}) + \frac{x^{4}}{8}(1-\frac{x^{2}}{2}+\frac{x^{4}}{8})$$

Now do simplification and eliminate terms that have order of higher than four:

$$f(x) = \frac{x^{2}}{2} - \frac{x^{4}}{4} - \frac{x^{2}}{2} + \frac{x^{4}}{4} + \frac{x^{4}}{8}$$

So final form is:

$$f(x) = \frac{x^{4}}{8} + O(x^{5})$$

Check the results from Wolfram Alpha here: https://www.wolframalpha.com/input/?i=1-%5Csqrt%281%2Bx%5E2%29%2F%281%2Bx%5E2%2F2%29

If you want it up to eighth order accuracy, it's the solution from Wolfram Alpha:

$$f(x) = \frac{x^{4}}{8} - \frac{x^{6}}{8} + \frac{13 x^{8}}{128} + O(x^{9})$$

| cite | improve this answer | |
  • $\begingroup$ Thanks for the explanation but what I tried to ask was "If I compute the function numerically for let's say 10^(-10)<x<1, then should I manipulate the function with Taylor series to have a better precision?". I am going to edit the question to be more clear. $\endgroup$ – Ekrem Sep 19 '19 at 16:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.