Opened 5 years ago

Closed 5 years ago

#7290 closed Bugs (fixed)

complex acos is occasionally wrong

Reported by: Stephen Montgomery-Smith <stephen@…> Owned by: John Maddock
Milestone: To Be Determined Component: math
Version: Boost 1.52.0 Severity: Problem
Keywords: Cc:


I am finding that the acos function is getting some of them wrong. For example, your program evaluates acos(1.00000002785941 + I*5.72464869028403e-200) as 0 - I*0.000236048349018331 whereas it should be 2.42520172707401e-196 - I*0.000236048349018331.

Looking at, I am somewhat sure that the mistake is in the last line of this code fragment:

This is the Hull et al exception handling code from Fig 6 of their paper: if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1))) {

if(x < one) {

real = std::acos(x); imag = y / std::sqrt(xp1*(one-x));

} else {

real = 0;

For asin, setting real = half_pi does just fine. But for acos, real should be something extremely small, but definitely not 0.

Attachments (0)

Change History (5)

comment:1 Changed 5 years ago by John Maddock

Component: Nonemath
Owner: set to John Maddock

comment:2 Changed 5 years ago by Stephen Montgomery-Smith <stephen@…>

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 @@
-               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));
+                  real = y / x;
                   imag = log_two + std::log(x);

comment:3 Changed 5 years ago by Stephen Montgomery-Smith <stephen@…>

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.

comment:4 Changed 5 years ago by Stephen Montgomery-Smith <stephen@…>

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.

comment:5 Changed 5 years ago by John Maddock

Resolution: fixed
Status: newclosed

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

Modify Ticket

Change Properties
Set your email in Preferences
as closed The owner will remain John Maddock.
The resolution will be deleted.

Add Comment

E-mail address and name can be saved in the Preferences.

Note: See TracTickets for help on using tickets.