docs: @Reference(authors = "J. M. Bernando", //
        title = "Algorithm AS 103: Psi (Digamma) Function", //
        booktitle = "Statistical Algorithms", //
        url = "https://doi.org/10.2307/2347257", //
        bibkey = "doi:10.2307/2347257")
    public static double digamma(double x) {
      if(!(x > 0)) {
        return Double.NaN;
      }
      // Method of equation 5:
      if(x <= 1e-5) {
        return -EULERS_CONST - 1. / x;
      }
      // Method of equation 4:
      else if(x > 49.) {
        final double ix2 = 1. / (x * x);
        // Partial series expansion
        return FastMath.log(x) - 0.5 / x - ix2 * ((1.0 / 12.) + ix2 * (1.0 / 120. - ix2 / 252.));
        // + O(x^8) error
      }
      else {
        // Stirling expansion
        return digamma(x + 1.) - 1. / x;
      }
    }