Search in sources :

Example 16 with Complex

use of org.hipparchus.complex.Complex in project symja_android_library by axkr.

the class BesselJS method airyAi.

/**
 * See: <a href="http://dlmf.nist.gov/9.2.ii"></a> <a href="http://dlmf.nist.gov/9.6.i"></a>
 *
 * @param x
 * @return
 */
public static Complex airyAi(Complex x) {
    if (F.isZero(x)) {
        return new Complex(1.0 / Math.pow(3.0, 2.0 / 3.0) / GammaJS.gamma(2.0 / 3.0));
    }
    if (x.getReal() < 0) {
        Complex xMinus = x.negate();
        Complex z = xMinus.pow(1.5).multiply(2.0 / 3.0);
        return xMinus.sqrt().divide(3.0).multiply(besselJ(1.0 / 3.0, z).add(besselJ(-1.0 / 3.0, z)));
    }
    Complex z = x.pow(1.5).multiply(2.0 / 3.0);
    return besselK(1.0 / 3.0, z).multiply(x.divide(3.0).sqrt()).multiply(1 / Math.PI);
}
Also used : Complex(org.hipparchus.complex.Complex)

Example 17 with Complex

use of org.hipparchus.complex.Complex in project symja_android_library by axkr.

the class EllipticIntegralsJS method carlsonRJ.

private static Complex carlsonRJ(double x, double y, double z, double p, double tolerance) {
    // adapted from mpmath / elliptic.py
    double xm = x;
    double ym = y;
    double zm = z;
    double pm = p;
    double Am = (x + y + z + 2 * p) / 5.0;
    double A0 = Am;
    double delta = (p - x) * (p - y) * (p - z);
    double Q = Math.pow(.25 * tolerance, -1.0 / 6.0) * Math.max(Math.abs(A0 - x), Math.max(Math.abs(A0 - y), Math.max(Math.abs(A0 - z), Math.abs(A0 - p))));
    double m = 0;
    double g = 0.25;
    double pow4 = 1;
    Complex S = Complex.ZERO;
    while (true) {
        double sx = Math.sqrt(xm);
        double sy = Math.sqrt(ym);
        double sz = Math.sqrt(zm);
        double sp = Math.sqrt(pm);
        double lm = sx * sy + sx * sz + sy * sz;
        double Am1 = (Am + lm) * g;
        xm = (xm + lm) * g;
        ym = (ym + lm) * g;
        zm = (zm + lm) * g;
        pm = (pm + lm) * g;
        double dm = (sp + sx) * (sp + sy) * (sp + sz);
        double em = delta * Math.pow(4.0, -3.0 * m) / Math.pow(dm, 2);
        if (pow4 * Q < Math.abs(Am)) {
            break;
        }
        Complex T = Complex.valueOf(carlsonRC(1, 1 + em)).multiply(pow4 / dm);
        S = S.add(T);
        pow4 *= g;
        m += 1;
        Am = Am1;
    }
    double t = Math.pow(2, -2 * m) / Am;
    double X = (A0 - x) * t;
    double Y = (A0 - y) * t;
    double Z = (A0 - z) * t;
    double P = (-X - Y - Z) / 2.0;
    double E2 = X * Y + X * Z + Y * Z - 3 * Math.pow(P, 2);
    double E3 = X * Y * Z + 2 * E2 * P + 4 * Math.pow(P, 3);
    double E4 = (2 * X * Y * Z + E2 * P + 3 * Math.pow(P, 3)) * P;
    double E5 = X * Y * Z * Math.pow(P, 2);
    P = 24024 - 5148 * E2 + 2457 * Math.pow(E2, 2) + 4004 * E3 - 4158 * E2 * E3 - 3276 * E4 + 2772 * E5;
    double v1 = Math.pow(g, m) * Math.pow(Am, -1.5) * P / 24024.0;
    return S.multiply(6.0).add(v1);
}
Also used : Complex(org.hipparchus.complex.Complex)

Example 18 with Complex

use of org.hipparchus.complex.Complex in project symja_android_library by axkr.

the class EllipticIntegralsJS method carlsonRJ.

private static Complex carlsonRJ(Complex x, Complex y, Complex z, Complex p, double tolerance) {
    // if ( isComplex(x) || isComplex(y) || isComplex(z) || isComplex(p) ) {
    Complex xm = x;
    Complex ym = y;
    Complex zm = z;
    Complex pm = p;
    Complex Am = x.add(y).add(z).add(p.multiply(2)).divide(5.0);
    Complex A0 = Am;
    Complex delta = p.subtract(x).multiply(p.subtract(y)).multiply(p.subtract(z));
    double Q = Math.pow(0.25 * tolerance, -1.0 / 6.0) * Math.max(A0.subtract(x).norm(), Math.max(A0.subtract(y).norm(), Math.max(A0.subtract(z).norm(), A0.subtract(p).norm())));
    double m = 0.0;
    double g = 0.25;
    double pow4 = 1.0;
    Complex S = Complex.ZERO;
    while (true) {
        Complex sx = xm.sqrt();
        Complex sy = ym.sqrt();
        Complex sz = zm.sqrt();
        Complex sp = pm.sqrt();
        Complex lm = sx.multiply(sy).add(sx.multiply(sz)).add(sy.multiply(sz));
        Complex Am1 = Am.add(lm).multiply(g);
        xm = xm.add(lm).multiply(g);
        ym = ym.add(lm).multiply(g);
        zm = zm.add(lm).multiply(g);
        pm = pm.add(lm).multiply(g);
        Complex dm = sp.add(sx).multiply(sp.add(sy)).multiply(sp.add(sz));
        Complex em = dm.reciprocal().multiply(dm.reciprocal()).multiply(delta).multiply(Math.pow(4.0, -3.0 * m));
        if (pow4 * Q < Am.norm()) {
            break;
        }
        Complex T = carlsonRC(Complex.ONE, em.add(1)).multiply(pow4).multiply(dm.reciprocal());
        S = S.add(T);
        pow4 *= g;
        m += 1;
        Am = Am1;
    }
    Complex t = Am.reciprocal().multiply(Math.pow(2, -2 * m));
    Complex X = A0.subtract(x).multiply(t);
    Complex Y = A0.subtract(y).multiply(t);
    Complex Z = A0.subtract(z).multiply(t);
    Complex P = X.add(Y.add(Z)).divide(-2);
    Complex E2 = X.multiply(Y).add(X.multiply(Z)).add(Y.multiply(Z)).add(P.multiply(P).multiply(-3));
    Complex E3 = X.multiply(Y).multiply(Z).add(E2.multiply(P).multiply(2)).add(P.multiply(P).multiply(P).multiply(4));
    Complex E4 = X.multiply(Y).multiply(Z).multiply(2).add(E2.multiply(P)).add(P.multiply(P).multiply(P).multiply(3)).multiply(P);
    Complex E5 = X.multiply(Y).multiply(Z).multiply(P).multiply(P);
    P = E2.multiply(-5148).add(E2.multiply(E2).multiply(2457)).add(E3.multiply(4004)).add(E2.multiply(E3).multiply(-4158)).add(E4.multiply(-3276)).add(E5.multiply(2772)).add(24024);
    Complex v1 = Am.pow(-1.5).multiply(Math.pow(g, m)).multiply(P).multiply(1.0 / 24024.0);
    return S.multiply(6.0).add(v1);
}
Also used : Complex(org.hipparchus.complex.Complex)

Example 19 with Complex

use of org.hipparchus.complex.Complex in project symja_android_library by axkr.

the class EllipticIntegralsJS method ellipticPi.

public static Complex ellipticPi(double n, double x, double m) {
    if (n > 1 && Math.abs(x) > Math.asin(1 / Math.sqrt(n))) {
        return ellipticPi(new Complex(n), new Complex(x), new Complex(m));
    }
    if (m > 1 && Math.abs(x) > Math.asin(1 / Math.sqrt(m))) {
        return ellipticPi(new Complex(n), new Complex(x), new Complex(m));
    }
    // return Complex.valueOf(LegendreEllipticIntegral.bigPi(n, x, m));
    Complex period = Complex.ZERO;
    if (Math.abs(x) > Math.PI / 2.0) {
        long p = Math.round(x / Math.PI);
        x = x - p * Math.PI;
        period = ellipticPi(n, Math.PI / 2.0, m).multiply(p + p);
    }
    double sinX = Math.sin(x);
    double cosX = Math.cos(x);
    double sqrSinX = sinX * sinX;
    double sqrCosX = cosX * cosX;
    double p3SqrSinX = sqrSinX * sinX;
    double mSqrSinX = 1.0 - m * sqrSinX;
    double nSqrSinX = 1.0 - n * sqrSinX;
    if (mSqrSinX < 0) {
        return carlsonRF(new Complex(sqrCosX), new Complex(mSqrSinX), Complex.ONE).multiply(sinX).add(carlsonRJ(new Complex(sqrCosX), new Complex(mSqrSinX), Complex.ONE, new Complex(nSqrSinX)).multiply(n / 3.0 * p3SqrSinX)).add(period);
    }
    return Complex.valueOf(carlsonRF(sqrCosX, mSqrSinX, 1)).multiply(sinX).add(Complex.valueOf(carlsonRJ(sqrCosX, mSqrSinX, 1, nSqrSinX)).multiply(n / 3.0 * p3SqrSinX)).add(period);
}
Also used : Complex(org.hipparchus.complex.Complex)

Example 20 with Complex

use of org.hipparchus.complex.Complex 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

Complex (org.hipparchus.complex.Complex)74 ArgumentTypeException (org.matheclipse.core.eval.exception.ArgumentTypeException)12 EvalEngine (org.matheclipse.core.eval.EvalEngine)8 IExpr (org.matheclipse.core.interfaces.IExpr)6 Config (org.matheclipse.core.basic.Config)4 F (org.matheclipse.core.expression.F)4 IAST (org.matheclipse.core.interfaces.IAST)4 Gamma (org.hipparchus.special.Gamma)3 Arithmetic (org.matheclipse.core.builtin.Arithmetic)3 IterationLimitExceeded (org.matheclipse.core.eval.exception.IterationLimitExceeded)3 ValidateException (org.matheclipse.core.eval.exception.ValidateException)3 S (org.matheclipse.core.expression.S)3 IComplex (org.matheclipse.core.interfaces.IComplex)3 ISymbol (org.matheclipse.core.interfaces.ISymbol)3 Math.abs (java.lang.Math.abs)2 ArrayList (java.util.ArrayList)2 Function (java.util.function.Function)2 IntFunction (java.util.function.IntFunction)2 DSFactory (org.hipparchus.analysis.differentiation.DSFactory)2 FiniteDifferencesDifferentiator (org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator)2