Modify

Ticket #7290 (closed Bugs: fixed)

Opened 20 months ago

Last modified 17 months ago

complex acos is occasionally wrong

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

Description

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 http://www.boost.org/doc/libs/1_51_0/boost/math/complex/acos.hpp, 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

Change History

comment:1 Changed 20 months ago by johnmaddock

  • Owner set to johnmaddock
  • Component changed from None to math

comment:2 Changed 20 months 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 @@
             }
             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);
                }
             }

comment:3 Changed 20 months 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 20 months 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 17 months ago by johnmaddock

  • Status changed from new to closed
  • Resolution set to fixed

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

View

Add a comment

Modify Ticket

Change Properties
<Author field>
Action
as closed
The resolution will be deleted. Next status will be 'reopened'
Author


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

 
Note: See TracTickets for help on using tickets.