id	summary	reporter	owner	description	type	status	milestone	component	version	severity	resolution	keywords	cc
11395	Beta Binonial Distriction/ log beta / boost::math::beta not usable when the result would be very small	Oli Quinet <quinet.olivier@…>	John Maddock	"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(...))
"	Feature Requests	closed	To Be Determined	math	Boost 1.57.0	Problem	obsolete