docs: @Reference(title = "Fast and accurate computation of binomial probabilities", //
        authors = "C. Loader", booktitle = "", //
        url = "http://projects.scipy.org/scipy/raw-attachment/ticket/620/loader2000Fast.pdf", //
        bibkey = "web/Loader00")
    public static double pmf(double x, int n, double p) {
      // Invalid values
      if(x < 0 || x > n) {
        return 0.;
      }
      // Extreme probabilities
      if(p <= 0.) {
        return x == 0 ? 1. : 0.;
      }
      if(p >= 1.) {
        return x == n ? 1. : 0.;
      }
      final double q = 1 - p;
      // FIXME: check for x to be integer, return 0 otherwise?
  
      // Extreme values of x
      if(x == 0) {
        if(p < .1) {
          return FastMath.exp(-devianceTerm(n, n * q) - n * p);
        }
        else {
          return FastMath.exp(n * FastMath.log(q));
        }
      }
      if(x == n) {
        if(p > .9) {
          return FastMath.exp(-devianceTerm(n, n * p) - n * q);
        }
        else {
          return FastMath.exp(n * FastMath.log(p));
        }
      }
      final double lc = stirlingError(n) - stirlingError(x) - stirlingError(n - x) - devianceTerm(x, n * p) - devianceTerm(n - x, n * q);
      final double f = (MathUtil.TWOPI * x * (n - x)) / n;
      return FastMath.exp(lc) / FastMath.sqrt(f);
    }