nd4j-common/src/main/java/org/nd4j/linalg/util/BigDecimalMath.java
static public BigDecimal pochhammer(final BigDecimal x, final int n) {
        /* reduce to interval near 1.0 with the functional relation, Abramowitz-Stegun 6.1.33
         */
        if (n < 0) {
            throw new ProviderException("Unimplemented pochhammer with negative index " + n);
        } else if (n == 0) {
            return BigDecimal.ONE;
        } else {
            /* internally two safety digits
             */
            BigDecimal xhighpr = scalePrec(x, 2);
            BigDecimal resul = xhighpr;


            double xUlpDbl = x.ulp().doubleValue();


            double xDbl = x.doubleValue();
            /* relative error of the result is the sum of the relative errors of the factors
             */


            double eps = 0.5 * xUlpDbl / Math.abs(xDbl);


            for (int i = 1; i < n; i++) {
                eps += 0.5 * xUlpDbl / Math.abs(xDbl + i);


                resul = resul.multiply(xhighpr.add(new BigDecimal(i)));
                final MathContext mcloc = new MathContext(4 + err2prec(eps));
                resul = resul.round(mcloc);


            }
            return resul.round(new MathContext(err2prec(eps)));


        }
    }
static public BigDecimal pochhammer(final BigDecimal x, final int n) {
        /* reduce to interval near 1.0 with the functional relation, Abramowitz-Stegun 6.1.33
         */
        if (n < 0) {
            throw new ProviderException("Unimplemented pochhammer with negative index " + n);
        } else if (n == 0) {
            return BigDecimal.ONE;
        } else {
            /* internally two safety digits
             */
            BigDecimal xhighpr = scalePrec(x, 2);
            BigDecimal resul = xhighpr;


            double xUlpDbl = x.ulp().doubleValue();


            double xDbl = x.doubleValue();
            /* relative error of the result is the sum of the relative errors of the factors
             */


            double eps = 0.5 * xUlpDbl / Math.abs(xDbl);


            for (int i = 1; i < n; i++) {
                eps += 0.5 * xUlpDbl / Math.abs(xDbl + i);


                resul = resul.multiply(xhighpr.add(new BigDecimal(i)));
                final MathContext mcloc = new MathContext(4 + err2prec(eps));
                resul = resul.round(mcloc);


            }
            return resul.round(new MathContext(err2prec(eps)));


        }
    }
Pochhammer’s function.

@param x The main argument.
@param n The non-negative index.
@return (x)_n = x(x+1)(x+2)*...*(x+n-1).