Boost C++ Libraries: Ticket #7290: complex acos is occasionally wrong

component changed; owner set

John Maddock — Sun, 26 Aug 2012 17:02:52 GMT

owner set to John Maddock
component None → math

Sun, 26 Aug 2012 18:31:37 GMT

It looks like the mistake is in Figure 6 of the paper by Hull, Fairgrieve and Tang. It looks like this patch fixes it.

diff -u boost-trunk/boost/math/complex/acos.hpp boost-trunk-new/boost/math/complex/acos.hpp
--- boost-trunk/boost/math/complex/acos.hpp     2012-08-26 18:24:02.000000000 +0000
+++ boost-trunk-new/boost/math/complex/acos.hpp 2012-08-26 18:25:48.000000000 +0000
@@ -172,14 +172,15 @@
             }
             else
             {
-               real = 0;
                if(((std::numeric_limits<T>::max)() / xp1) > xm1)
                {
                   // xp1 * xm1 won't overflow:
+                  real = y / std::sqrt(xm1*xp1);
                   imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
                }
                else
                {
+                  real = y / x;
                   imag = log_two + std::log(x);
                }
             }

Sun, 26 Aug 2012 19:02:08 GMT

Looking at page 327 of the paper by Hull et al, they say that the real part can be set to zero, because the absolute error will be negligible compared to the absolute error in the imaginary part. But everywhere else in the paper, they manage to control the relative error of the real and imaginary parts separately.

Looking at page 304, it seems that their goal is only to control the relative error of the real and imaginary parts together. So this wasn't really a mistake in their paper.

But I think the "higher goal" of controlling the relative error of the real and imaginary parts separately is worthy, so I still commend my patch to you.

Mon, 27 Aug 2012 01:22:30 GMT

Another thing - through experiment I have found that setting a_crossover to 10 instead of 1.5 works a little better. This is for asin as well as acos.

status changed; resolution set

John Maddock — Thu, 29 Nov 2012 13:27:36 GMT

status new → closed
resolution → fixed

(In [81624]) Lots of small patches. Update and regenerate docs. Fixes #7183. Fixes #7290. Fixes #7291. Fixes #7649. Refs #7492.