Search in sources :

Example 1 with ResultException

use of org.matheclipse.core.eval.exception.ResultException in project symja_android_library by axkr.

the class HypergeometricJS method hypergeometric2F1.

public static Complex hypergeometric2F1(Complex a, Complex b, Complex c, Complex x, double tolerance) {
    if (F.isFuzzyEquals(a, c, tolerance)) {
        return Complex.ONE.subtract(x).pow(b.negate());
    }
    if (F.isFuzzyEquals(b, c, tolerance)) {
        return Complex.ONE.subtract(x).pow(a.negate());
    }
    EvalEngine engine = EvalEngine.get();
    final int recursionLimit = engine.getRecursionLimit();
    try {
        if (recursionLimit > 0) {
            int counter = engine.incRecursionCounter();
            if (counter > recursionLimit) {
                // 
                RecursionLimitExceeded.throwIt(// 
                counter, F.Hypergeometric2F1(F.complexNum(a), F.complexNum(b), F.complexNum(c), F.complexNum(x)));
            }
        }
        // choose smallest absolute value of transformed argument
        // transformations from Abramowitz & Stegun p.559
        // fewer operations compared to dlmf.nist.gov/15.8
        double[] absArray = new double[] { // 
        x.norm(), // 
        x.divide(x.subtract(1)).norm(), // 
        new Complex(1).subtract(x).norm(), // 
        x.reciprocal().norm(), // 
        new Complex(1).subtract(x).reciprocal().norm(), new Complex(1).subtract(x.reciprocal()).norm() };
        double min = Double.POSITIVE_INFINITY;
        double newMin = Double.POSITIVE_INFINITY;
        int index = -1;
        for (int i = 0; i < absArray.length; i++) {
            newMin = Math.min(min, absArray[i]);
            if (newMin != min) {
                min = newMin;
                index = i;
            }
        }
        final Complex subtractCA = c.subtract(a);
        final Complex subtractCB = c.subtract(b);
        final Complex af = a;
        final Complex bf = b;
        final Complex cf = c;
        final Complex xf = x;
        switch(index) {
            case 0:
                break;
            case 1:
                return new Complex(1.0).subtract(x).pow(a.negate()).multiply(hypergeometric2F1(a, c.subtract(b), c, x.divide(x.subtract(1))));
            case 2:
                {
                    if (c.subtract(a.add(b)).isMathematicalInteger() || (subtractCA.isMathematicalInteger() && subtractCA.getReal() <= 0)) {
                        return complexAverage(v -> hypergeometric2F1(v, bf, cf, xf), af);
                    }
                    if (subtractCB.isMathematicalInteger() && subtractCB.getReal() <= 0) {
                        return complexAverage(v -> hypergeometric2F1(af, v, cf, xf), bf);
                    }
                    Complex t1 = Arithmetic.lanczosApproxGamma(c).multiply(Arithmetic.lanczosApproxGamma(c.subtract(a.add(b)))).multiply(Arithmetic.lanczosApproxGamma(subtractCA).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(c.subtract(b)).reciprocal()).multiply(hypergeometric2F1(a, b, a.add(b).add(c.negate()).add(1), new Complex(1).subtract(x)));
                    Complex t2 = new Complex(1).subtract(x).pow(c.subtract(a.add(b))).multiply(Arithmetic.lanczosApproxGamma(c)).multiply(Arithmetic.lanczosApproxGamma(a.add(b).subtract(c))).multiply(Arithmetic.lanczosApproxGamma(a).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(b).reciprocal()).multiply(hypergeometric2F1(subtractCA, c.subtract(b), a.add(a.negate()).add(b.negate()).add(1), new Complex(1).subtract(x)));
                    return t1.add(t2);
                }
            case 3:
                {
                    if (a.subtract(b).isMathematicalInteger() || (subtractCA.isMathematicalInteger() && subtractCA.getReal() <= 0)) {
                        return complexAverage(v -> hypergeometric2F1(v, bf, cf, xf), af);
                    }
                    if (subtractCB.isMathematicalInteger() && subtractCB.getReal() <= 0) {
                        return complexAverage(v -> hypergeometric2F1(af, v, cf, xf), bf);
                    }
                    Complex t1 = Arithmetic.lanczosApproxGamma(c).multiply(Arithmetic.lanczosApproxGamma(b.subtract(a))).multiply(Arithmetic.lanczosApproxGamma(b).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(subtractCA).reciprocal()).multiply(x.negate().pow(a.negate())).multiply(hypergeometric2F1(a, a.add(1).add(c.negate()), a.add(1).add(b.negate()), x.reciprocal()));
                    Complex t2 = Arithmetic.lanczosApproxGamma(c).multiply(Arithmetic.lanczosApproxGamma(a.subtract(b))).multiply(Arithmetic.lanczosApproxGamma(a).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(c.subtract(b)).reciprocal()).multiply(x.negate().pow(b.negate())).multiply(hypergeometric2F1(b, b.add(1).add(c.negate()), b.add(1).add(a.negate()), x.reciprocal()));
                    return t1.add(t2);
                }
            case 4:
                {
                    if (a.subtract(b).isMathematicalInteger() || (subtractCA.isMathematicalInteger() && subtractCA.getReal() <= 0)) {
                        return complexAverage(v -> hypergeometric2F1(v, bf, cf, xf), af);
                    }
                    if (subtractCB.isMathematicalInteger() && subtractCB.getReal() <= 0) {
                        return complexAverage(v -> hypergeometric2F1(af, v, cf, xf), bf);
                    }
                    Complex t1 = new Complex(1.0).subtract(x).pow(a.negate()).multiply(Arithmetic.lanczosApproxGamma(c)).multiply(Arithmetic.lanczosApproxGamma(b.subtract(a))).multiply(Arithmetic.lanczosApproxGamma(b).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(subtractCA).reciprocal()).multiply(hypergeometric2F1(a, c.subtract(b), a.add(b.negate()).add(1), new Complex(1).subtract(x).reciprocal()));
                    Complex t2 = new Complex(1).subtract(x).pow(b.negate()).multiply(Arithmetic.lanczosApproxGamma(c)).multiply(Arithmetic.lanczosApproxGamma(a.subtract(b))).multiply(Arithmetic.lanczosApproxGamma(a).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(c.subtract(b)).reciprocal()).multiply(hypergeometric2F1(b, subtractCA, b.add(a.negate()).add(1), new Complex(1).subtract(x).reciprocal()));
                    return t1.add(t2);
                }
            case 5:
                {
                    if (c.subtract(a.add(b)).isMathematicalInteger() || (subtractCA.isMathematicalInteger() && subtractCA.getReal() <= 0)) {
                        return complexAverage(v -> hypergeometric2F1(v, bf, cf, xf), af);
                    }
                    if (subtractCB.isMathematicalInteger() && subtractCB.getReal() <= 0) {
                        return complexAverage(v -> hypergeometric2F1(af, v, cf, xf), bf);
                    }
                    Complex t1 = Arithmetic.lanczosApproxGamma(c).multiply(Arithmetic.lanczosApproxGamma(c.subtract(a.add(b)))).multiply(Arithmetic.lanczosApproxGamma(subtractCA).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(c.subtract(b)).reciprocal()).multiply(x.pow(a.negate())).multiply(hypergeometric2F1(a, a.add(c.negate()).add(1), a.add(b).add(c.negate()).add(1), new Complex(1).subtract(x.reciprocal())));
                    Complex t2 = Arithmetic.lanczosApproxGamma(c).multiply(Arithmetic.lanczosApproxGamma(a.add(b).subtract(c))).multiply(Arithmetic.lanczosApproxGamma(a).reciprocal()).multiply(Arithmetic.lanczosApproxGamma(b).reciprocal()).multiply(new Complex(1).subtract(x).pow(c.subtract(a.add(b)))).multiply(x.pow(a.subtract(c))).multiply(hypergeometric2F1(subtractCA, new Complex(1).subtract(a), c.add(a.negate()).add(b.negate()).add(1), new Complex(1).subtract(x.reciprocal())));
                    return t1.add(t2);
                }
        }
        if (c.isMathematicalInteger() && c.getReal() <= 0) {
            throw new ResultException(F.CComplexInfinity);
        // throw new ArgumentTypeException("hypergeometric function pole");
        }
        Complex s = Complex.ONE;
        Complex p = Complex.ONE;
        int i = 1;
        long iterationLimit = engine.getIterationLimit();
        while (Math.abs(p.getReal()) > tolerance || Math.abs(p.getImaginary()) > tolerance) {
            p = p.multiply(x).multiply(a).multiply(b).multiply(c.reciprocal()).divide(i);
            s = s.add(p);
            a = a.add(1);
            b = b.add(1);
            c = c.add(1);
            if (i++ > iterationLimit && iterationLimit > 0) {
                IterationLimitExceeded.throwIt(i, S.Hypergeometric2F1);
            }
        }
        return s;
    } finally {
        if (recursionLimit > 0) {
            engine.decRecursionCounter();
        }
    }
}
Also used : Gamma(org.hipparchus.special.Gamma) EvalEngine(org.matheclipse.core.eval.EvalEngine) ResultException(org.matheclipse.core.eval.exception.ResultException) F(org.matheclipse.core.expression.F) Complex(org.hipparchus.complex.Complex) Config(org.matheclipse.core.basic.Config) Math.abs(java.lang.Math.abs) Function(java.util.function.Function) ArrayList(java.util.ArrayList) S(org.matheclipse.core.expression.S) Arithmetic(org.matheclipse.core.builtin.Arithmetic) RecursionLimitExceeded(org.matheclipse.core.eval.exception.RecursionLimitExceeded) IterationLimitExceeded(org.matheclipse.core.eval.exception.IterationLimitExceeded) IntFunction(java.util.function.IntFunction) ArgumentTypeException(org.matheclipse.core.eval.exception.ArgumentTypeException) EvalEngine(org.matheclipse.core.eval.EvalEngine) ResultException(org.matheclipse.core.eval.exception.ResultException) Complex(org.hipparchus.complex.Complex)

Aggregations

Math.abs (java.lang.Math.abs)1 ArrayList (java.util.ArrayList)1 Function (java.util.function.Function)1 IntFunction (java.util.function.IntFunction)1 Complex (org.hipparchus.complex.Complex)1 Gamma (org.hipparchus.special.Gamma)1 Config (org.matheclipse.core.basic.Config)1 Arithmetic (org.matheclipse.core.builtin.Arithmetic)1 EvalEngine (org.matheclipse.core.eval.EvalEngine)1 ArgumentTypeException (org.matheclipse.core.eval.exception.ArgumentTypeException)1 IterationLimitExceeded (org.matheclipse.core.eval.exception.IterationLimitExceeded)1 RecursionLimitExceeded (org.matheclipse.core.eval.exception.RecursionLimitExceeded)1 ResultException (org.matheclipse.core.eval.exception.ResultException)1 F (org.matheclipse.core.expression.F)1 S (org.matheclipse.core.expression.S)1