[Boost-bugs] [Boost C++ Libraries] #11395: boost::math::beta not working for large a, b values

Subject: [Boost-bugs] [Boost C++ Libraries] #11395: boost::math::beta not working for large a, b values
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2015-06-12 07:37:31

#11395: boost::math::beta not working for large a,b values
-------------------------------------------+-------------------------
 Reporter:  Oli Quinet <quinet.olivier@…>  |      Owner:  johnmaddock
     Type:  Patches                        |     Status:  new
Milestone:  To Be Determined               |  Component:  math
  Version:  Boost 1.57.0                   |   Severity:  Problem
 Keywords:                                 |
-------------------------------------------+-------------------------
 Because the computation inside beta_imp is not done in the logarithm scale
 it returns zero instead of the correct value.  Proposed corrected
 implementation:

 // Lanczos calculation:
 T agh = a + Lanczos::g() - T(0.5);
 T bgh = b + Lanczos::g() - T(0.5);
 T cgh = c + Lanczos::g() - T(0.5);
 result = Lanczos::lanczos_sum_expG_scaled(a) *
 Lanczos::lanczos_sum_expG_scaled(b) / Lanczos::lanczos_sum_expG_scaled(c);

 // If a and b were originally less than 1 we need to scale the result:
 result *= prefix;
 result = log(result);

 result += (.5*(1-log(bgh)));

 T ambh = a - T(0.5) - b;
 if((fabs(b * ambh) < (cgh * 100)) && (a > 100))
 {
   // Special case where the base of the power term is close to 1
   // compute (1+x)^y instead:
   result += (ambh * boost::math::log1p(-b / cgh, pol));
   result += (b*(log(agh)+log(bgh)-2*log(cgh)));
 }
 else
 {
   result += (a-T(0.5))*(log(agh)-log(cgh))+(b*(log(bgh)-log(cgh)));
 }

 return exp(result);

 Also, it would be good to also propose a 'l_beta_imp' function returning
 the logarithm of beta_imp, i.e., not doing the 'exp(result)'.  The same
 can be done for 'binomial_coefficient' by introducing a
 'l_binomial_coefficient' function.

 The main interest is the combination of those functions without problems,
 e.g.:
    binomial_coefficient(...) * beta(...) / beta(...) can be converted to:
    exp(l_binomial_coefficient(...) + l_beta(...) - l_beta(...))

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/11395>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.