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: | 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 (0)
Change History (5)
comment:1 Changed 5 years ago by johnmaddock
- Component changed from None to math
- Owner set to johnmaddock
comment:2 Changed 5 years ago by Stephen Montgomery-Smith <stephen@…>
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.
It looks like the mistake is in Figure 6 of the paper by Hull, Fairgrieve and Tang. It looks like this patch fixes it.