extreme_value_distribution pdf fails for large deviates

[Paul McClellan <paulm@…>]

extreme_value_distribution<double> pdf function raises an error when the normalized deviate (a-x)/b is too large positive.

For example, when a = 0.0, b = 1.0, and x = -1.0e3 the pdf computation:

  result = exp((a-x)/b) * exp(-exp((a-x)/b)) / b;

results in the term exp((a-x)/b) = exp(-x) overflowing to +inf, and the term exp(-exp((a-x)/b)) / b = exp(-exp(-x)) underflowing to 0.0 and the product of the terms being indeterminate.

My workaround is:

  // evmax distribution object:
  boost::math::extreme_value evmax(dL, dS);
  // BUGFIX Set pdf(x;a,b) to zero for (a-x)/b > 37.
  // This avoids an indeterminant computation.
  return (dL-dX)/dS > 37.0 ? 0.0 : pdf(evmax, dX);

[John Maddock]

Fixed in ​https://github.com/boostorg/math/commit/545961746cb0d67ea20e27de98c459b3ec0ade71