I want to compute the difference $\Delta f(x_1,x_2) = f(x_1)-f(x_2)$ of a smooth function $f(x)$ at two points $x_1$ and $x_2$ which are close to each other. The magnitude of the expected result, $|\Delta f(x_1,x_2)|$, is often small compared to $|f(x_1)|$ and $|f(x_2)|$ and several significant digits are lost if I evaluate the expression $f(x_1)-f(x_2)$ straightforwardly. What is a recommended way to compute $\Delta f(x_1,x_2)$ accurately?

I have an analytical expression of $f(x)$. If a detail matters, it has a form like $f(x) = \exp(-a x^2) \sin(x)$.

Perhaps one method would be to use the Taylor expansion of $f(x)$, but I am afraid that the convergence might not be very fast. Also, writing down the higher order derivatives is tedious. I'd like to know something easier and more efficient.

  • 1
    $\begingroup$ Does your compute platform offer FMA (fused multiply-add)? If you program in C or C++, this should be available as the standard math function fma(). Support for FMA would allow efficient double-double computation to improve the accuracy of the computation. $\endgroup$ – njuffa Aug 16 '17 at 4:05
  • $\begingroup$ What do you need these differences for? Depending on your use case there might be more efficient methods. $\endgroup$ – H. Rittich Aug 16 '17 at 9:53
  • 2
    $\begingroup$ You also should consider using extended/arbitrary-precision floating point numbers for the operation. You can cast the result back to a fundamental floating point type at the end of the operation. It will result in a slower computation, but may be what you're looking for. $\endgroup$ – Tyler Olsen Aug 16 '17 at 12:34
  • 2
    $\begingroup$ @norio Sorry, while I have rudimentary Fortran reading skills, I don't write Fortran code, and I don't think Fortran 2008 provides an FMA built-in. There are double-double implementations in Fortran, I think (not sure), so maybe look for those. double-double is basically a computationally cheap approximation to proper quadruple-precision arithmetic, and use of FMA allows particularly efficient implementations. $\endgroup$ – njuffa Aug 16 '17 at 16:44
  • 1
    $\begingroup$ @norio In this case, I do not think there is an easy way around the cancellation error. I assume that you have to build the table only once. In this case, speed is probably not too critical and you can use high precision arithmetic, as already suggested. For Fortran implementations see here: stackoverflow.com/questions/15226606 $\endgroup$ – H. Rittich Aug 16 '17 at 20:20

Your Answer

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

Browse other questions tagged or ask your own question.