As reported in scipy/scipy#15101, there is a potential divide-by-zero in ibeta_power_terms L392-L393 for some values of n, p, and k, e.g., n=1000, p=0.01, k=996. This may not be the only instance of a possible divide-by-zero for these parameter choices.

Suggested patch would be to avoid the logarithm calculation on line 393 if p1 == 0. We are also curious why this underflow is not reported according to boost::math policy or why we can't reproduce any observable issue when run as a simple C++ application (GCC 9), e.g.:

#include <iostream>
#include "boost/math/distributions/binomial.hpp"

int main() {
  constexpr int n = 1000;
  constexpr double p = .01;
  constexpr int k = 996;
  boost::math::binomial_distribution<double> b(n, p);
  std::cout << std::setprecision(16) << "ans is " << boost::math::pdf(b, k) << std::endl;
  return 0;
}