Ticket #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
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.