Boost C++ Libraries: Ticket #5095: incorrect results from hypergeometric pdf

attachment set

Wed, 19 Jan 2011 17:06:17 GMT

attachment → hyper.cpp

example

Wed, 19 Jan 2011 17:26:05 GMT

As a workaround, if the lgamma backup version of hypergeometric_pdf_lanczos_imp is used instead, correct results are obtained.

In /usr/local/include/boost/math/distributions/detail/hypergeometric_pdf.hpp: 414c414 < result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy()); ---

result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, lanczos::undefined_lanczos(), forwarding_policy());

Wed, 19 Jan 2011 17:27:13 GMT

Properly formatted diff of boost/math/distributions/detail/hypergeometric_pdf.hpp

414c414
<       result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy());
---
>       result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, lanczos::undefined_lanczos(), forwarding_policy());

John Maddock — Thu, 20 Jan 2011 12:35:39 GMT

I've been unable to reproduce here either with VC10 or with gcc-4.4 on Ubuntu Linux, what compiler/platform are you on?

The output I see from your program is: 0.00640003 0.00639305 0.00638443 which I believe agrees with the online calculator.

anonymous — Thu, 20 Jan 2011 13:01:34 GMT

gcc version 4.4.3 (Ubuntu 4.4.3-4ubuntu5) on Ubuntu 10.4

Also, on Ubuntu 10.10 with boost version 1.40 and gcc version 4.4.5 (Ubuntu/Linaro 4.4.4-14ubuntu5) and OS X (10.6) with boost version 1.41 and gcc version 4.2.1 (Apple Inc. build 5664).

Perhaps for some reason your installation is already falling back on the lgamma version of hypergeometric_pdf_lanczos_imp?

anonymous — Thu, 20 Jan 2011 13:22:36 GMT

No it's calling for example:

hypergeometric_pdf_lanczos_imp<long double,boost::math::lanczos::lanczos13m53,boost::math::policies::policy<boost::math::policies::promote_float<0>,boost::math::policies::promote_double<0>,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy> >(double formal, unsigned int x, unsigned int r, unsigned int n, unsigned int N, double formal, double formal);

I notice that you're using an out of date Boost distribution, can you:

Please try the last release and if that still fails, then:
Let me have the program output when BOOST_MATH_INSTRUMENT is #defined when building?

Thanks, John.

attachment set

Thu, 20 Jan 2011 14:48:33 GMT

attachment → BMIdump

program output with BOOST_MATH_INSTRUMENT defined

Thu, 20 Jan 2011 14:49:13 GMT

The Ubuntu 10.04 box has boost 1.45 on it. I also tried it on two other machines with the version of boost that they already have. I have also grabbed and built the svn trunk on the OS X machine, and tried removing /usr/local/include/boost and reinstalling on the 10.4 machine. I still get the failure.

I have attached the output with BOOST_MATH_INSTRUMENT defined.

Thanks.

anonymous — Thu, 20 Jan 2011 16:46:43 GMT

Thanks for trying that, I really appreciate the help you're putting in to get to the bottom of this, unfortunately I still can't see what's gone wrong (other than something has), I'm also confused because there's really nothing in the code that should rely on the evenness of the random variable!

I'm attaching an updated hypergeometric_pdf.hpp with a lot more instrumentation, can I get you to try again? Just the failing case this time - otherwise the output is going to be huge! :(

Many thanks,

Still confused yours, John.

attachment set

anonymous — Thu, 20 Jan 2011 16:47:43 GMT

attachment → hypergeometric_pdf.hpp

anonymous — Thu, 20 Jan 2011 16:49:06 GMT

Forgot to say - if the program output is too large to attach here - zip it up and mail me direct at john at johnmaddock.co.uk

Thanks! John.

attachment set

Thu, 20 Jan 2011 16:54:40 GMT

attachment → BMIdumpmore.gz

even more instrumentation (bad value only, gzipped)

Thu, 20 Jan 2011 17:07:57 GMT

Attached. Here are some interesting results that may partially explain the reproducibility problem:

dkoes@quasar:~/tmp$ g++  hyper.cpp; ./a.out
1.11519e-09
dkoes@quasar:~/tmp$ g++ -m32 hyper.cpp; ./a.out
0.00639305
dkoes@quasar:~/tmp$ g++ -mfpmath=387 hyper.cpp; ./a.out
0.00639305
dkoes@quasar:~/tmp$ g++ -mfpmath=sse hyper.cpp; ./a.out
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse hyper.cpp; ./a.out
In file included from /usr/local/include/boost/config/no_tr1/cmath.hpp:21,
                 from /usr/local/include/boost/math/policies/error_handling.hpp:15,
                 from /usr/local/include/boost/math/distributions/detail/common_error_handling.hpp:12,
                 from /usr/local/include/boost/math/distributions/hypergeometric.hpp:12,
                 from hyper.cpp:1:
/usr/include/c++/4.4/cmath: In function ‘double std::abs(double)’:
/usr/include/c++/4.4/cmath:94: error: SSE register return with SSE disabled
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse2 hyper.cpp; ./a.out
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse3 hyper.cpp; ./a.out
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse4 hyper.cpp; ./a.out
1.11519e-09
dkoes@quasar:~/tmp$ g++ -DBOOST_MATH_SPECIAL_FUNCTIONS_LANCZOS_SSE2 hyper.cpp; ./a.out
1.11519e-09

anonymous — Thu, 20 Jan 2011 18:33:12 GMT

Thanks for that, some of the values in the "exponents" table are getting truncated, I'm attaching an updated header with some extra typecasts present to try and prevent this, can you give this a try?

Thanks, John.

attachment set

anonymous — Thu, 20 Jan 2011 18:33:49 GMT

attachment → hypergeometric_pdf.2.hpp

Thu, 20 Jan 2011 18:37:38 GMT

Bingo. That fixed it. Thanks!

John Maddock — Fri, 21 Jan 2011 12:03:51 GMT

(In [68347]) Added note about fixed bug. Refs #5095.

status changed; resolution set

anonymous — Fri, 21 Jan 2011 12:53:58 GMT

status new → closed
resolution → fixed

This was fixed in revision [68346].

John Maddock — Sat, 22 Jan 2011 17:31:32 GMT

(In [68367]) Merge fix for issue 5095. Refs #5095.