Examples with Gamma org.hipparchus.special.Gamma used on opensource projects

Example 1 with Gamma

use of org.hipparchus.special.Gamma in project symja_android_library by axkr.

the class GammaJS method gamma.

/**
 * Incomplete gamma function.
 *
 * @param x
 * @param y
 * @return
 */
public static Complex gamma(Complex x, Complex y) {
    // patch lower end or evaluate exponential integral independently
    if (F.isZero(x)) {
        if (F.isZero(y)) {
            throw new ArgumentTypeException("Gamma function pole");
        }
        // if (Complex.equals(x, Complex.ZERO, Config.SPECIAL_FUNCTIONS_TOLERANCE)) {
        // // taylorSeries => (-EulerGamma - Log(y)) + x - 1/4 * x^2 + 1/18 * x^3 - 1/96 * x^4 + 1/600
        // * x^5
        // Complex result = y.log().add(ConstantDefinitions.EULER_GAMMA).negate();
        // double[] coeff = new double[] { 1.0, -0.25, 1.0 / 18.0, -1.0 / 96.0, 1.0 / 600.0 };
        // Complex yPow = y;
        // for (int i = 0; i < coeff.length; i++) {
        // result = result.add(yPow.multiply(coeff[i]));
        // yPow = yPow.multiply(y);
        // }
        // return result;
        Complex yLogNegate = y.log().negate();
        Complex result = // 
        expIntegralEi(y.negate()).negate().add(// 
        y.negate().log().subtract(y.reciprocal().negate().log()).multiply(0.5)).add(yLogNegate);
        if (F.isZero(y.getImaginary()) && y.getReal() > 0.0) {
            return new Complex(result.getReal());
        }
        return result;
    }
    double delta = 1e-5;
    if (x.norm() < delta) {
    // TODO implement case for abs value near 0
    // return taylorSeries( t => gamma(t,y), mul( x, delta/x.abs( ) ), 2.0)(x);
    }
    // dlmf.nist.gov/8.4.15
    double xRe = x.getReal();
    if (xRe < 0.0 && F.isNumIntValue(xRe) && F.isZero(x.getImaginary())) {
        // x is a negative integer
        final double n = -xRe;
        int iterationLimit = EvalEngine.get().getIterationLimit();
        final Complex t = // 
        y.negate().exp().multiply(ZetaJS.complexSummation(// 
        k -> new Complex(Math.pow(-1.0, k) * factorialInt(k)).divide(y.pow(k + 1)), 0.0, n - 1.0, iterationLimit));
        // dlmf.nist.gov/8.4.4
        final double plusMinusOne = Math.pow(-1.0, n);
        return gamma(Complex.ZERO, y).subtract(t).multiply(plusMinusOne / factorialInt(n));
    }
    return Arithmetic.lanczosApproxGamma(x).subtract(gamma(x, 0.0, y));
}