Examples with ArgumentTypeException org.matheclipse.core.eval.exception.ArgumentTypeException used on opensource projects

Example 6 with ArgumentTypeException

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

the class BesselJS method besselYZero.

/**
 * @param n a positive value
 * @param m
 * @return
 */
public static double besselYZero(double n, int m) {
    if (n < 0.0) {
        throw new ArgumentTypeException(n + " < 0 (negative order)");
    }
    // approximations from dlmf.nist.gov/10.21#vi
    double delta = Math.PI / 4.0;
    // if (false) {
    // // use derivative
    // double b = (m + n / 2.0 - 0.25) * Math.PI;
    // double e = b - (4.0 * (n * n) + 3.0) / (8.0 * b);
    // 
    // // return findRoot( x => diff( x => besselY(n,x), x ), [ e-delta, e+delta ] );
    // BisectionSolver solver = new BisectionSolver();
    // ISymbol x = F.Dummy("x");
    // IExpr function= F.Times(F.C1D2,F.Subtract(F.BesselY(F.num(-1.0 + n), x),F.BesselY(F.num(1.0 +
    // n), x)));
    // UnivariateDifferentiableFunction f = new UnaryNumerical(function, x, EvalEngine.get(), true);
    // return solver.solve(200, f, e - delta, e + delta);
    // } else {
    double a = (m + n / 2.0 - 0.75) * Math.PI;
    double e = a - (4.0 * (n * n) - 1.0) / (8.0 * a);
    BisectionSolver solver = new BisectionSolver();
    ISymbol x = F.Dummy("x");
    UnivariateDifferentiableFunction f = new UnaryNumerical(F.BesselY(F.num(n), x), x, EvalEngine.get(), true);
    return solver.solve(100, f, e - delta, e + delta);
// }
}

Example 7 with ArgumentTypeException

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

the class BesselJS method besselKDouble.

public static double besselKDouble(double n, double x) {
    final int useAsymptotic = 10;
    if (x > useAsymptotic) {
        return Math.sqrt(0.5 * Math.PI / x) * Math.exp(-x) * HypergeometricJS.hypergeometric2F0(n + 0.5, 0.5 - n, -0.5 / x);
    }
    if (x < 0) {
        throw new ArgumentTypeException(x + " < 0.0");
    }
    // based on dlmf.nist.gov/10.2.3
    if (F.isNumIntValue(n)) {
        FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(15, 0.01);
        ISymbol nSymbol = F.Dummy("n");
        UnivariateDifferentiableFunction f = differentiator.differentiate(new UnaryNumerical(F.BesselI(nSymbol, F.num(x)), nSymbol, EvalEngine.get()));
        DSFactory factory = new DSFactory(1, 1);
        double d1 = f.value(factory.variable(0, n)).getPartialDerivative(1);
        double d2 = f.value(factory.variable(0, -n)).getPartialDerivative(1);
        return ((d1 + d2) / 2.0) * Math.pow(-1.0, Math.rint(n + 1));
    }
    return (besselIDouble(-n, x) - besselIDouble(n, x)) * Math.PI / (2.0 * Math.sin(n * Math.PI));
}

Example 8 with ArgumentTypeException

use of org.matheclipse.core.eval.exception.ArgumentTypeException 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();
        }
    }
}

Example 9 with ArgumentTypeException

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

the class GammaJS method logGamma.

public static Complex logGamma(Complex x) {
    if (F.isNumIntValue(x.getReal()) && x.getReal() <= 0 && F.isZero(x.getImaginary())) {
        throw new ArgumentTypeException("Gamma function pole");
    }
    // reflection formula with modified Hare correction to imaginary part
    if (x.getReal() < 0.0) {
        Complex t = new Complex(Math.PI).divide(x.multiply(Math.PI).sin()).log().subtract(logGamma(x.negate().add(1.0)));
        double s = x.getImaginary() < 0.0 ? -1.0 : 1.0;
        double d = F.isZero(x.getImaginary()) ? 0.25 : 0;
        double k = Math.ceil(x.getReal() / 2.0 - 0.75 + d);
        return t.add(new Complex(0.0, 2.0 * s * k * Math.PI));
    }
    Complex t = x.add(5.24218750000000000);
    t = x.add(0.5).multiply(t.log()).subtract(t);
    Complex s = new Complex(0.999999999999997092);
    for (int j = 0; j < 14; j++) {
        s = s.add(x.add(j + 1).reciprocal().multiply(c[j]));
    }
    Complex u = t.add(s.divide(x).multiply(2.5066282746310005).log());
    // adjustment to keep imaginary part on same sheet
    if (s.getReal() < 0.0) {
        if (x.getImaginary() < 0.0 && s.divide(x).getImaginary() < 0) {
            u = u.add(new Complex(0.0, Math.PI + Math.PI));
        }
        if (x.getImaginary() > 0 && s.divide(x).getImaginary() > 0) {
            u = u.add(new Complex(0.0, -Math.PI + Math.PI));
        }
    }
    return u;
}

Example 10 with ArgumentTypeException

use of org.matheclipse.core.eval.exception.ArgumentTypeException 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));
}