Examples with S org.matheclipse.core.expression.S used on opensource projects

Example 1 with S

use of org.matheclipse.core.expression.S 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 2 with S

use of org.matheclipse.core.expression.S in project symja_android_library by axkr.

the class ZetaJS method hurwitzZeta.

/**
 * @param x
 * @param a
 * @return
 */
public static Complex hurwitzZeta(final Complex x, final Complex a) {
    // TODO INVALID at the moment
    if (x.getReal() == 1.0 && x.getImaginary() == 0.0) {
        throw new ArgumentTypeException("Hurwitz zeta pole");
    }
    // direct summation more accurate than dlmf.nist.gov/25.11.4 for positive a
    int iterationLimit = EvalEngine.get().getIterationLimit();
    if (a.getReal() < 0.0) {
        double m = -Math.floor(a.getReal());
        return hurwitzZeta(x, a.add(m)).add(summation(i -> a.add(i).pow(x.negate()), 0, m - 1.0, iterationLimit));
    }
    // Johansson arxiv.org/abs/1309.2877
    // recommendation of Vepstas, Efficient Algorithm, p.12
    int n = 15;
    // Euler-Maclaurin has differences of large values in left-hand plane
    boolean useArbitrary = x.getReal() < 0;
    // if ( useArbitrary ) {
    // 
    // }
    double switchForms = -5.0;
    if (x.getReal() < switchForms) {
        throw new ArgumentTypeException("Currently unsuppported complex Hurwitz zeta");
    }
    Complex S = summation(i -> a.add(i).pow(x.negate()), 0, n - 1, iterationLimit);
    Complex I = a.add(n).pow(Complex.ONE.subtract(x)).divide(x.subtract(1.0));
    Complex p = x.multiply(0.5).multiply(a.add(n).reciprocal());
    Complex t = p.multiply(bernoulliInt(2));
    int i = 1;
    // converges rather quickly
    while (Math.abs(p.getReal()) > Config.SPECIAL_FUNCTIONS_TOLERANCE || Math.abs(p.getImaginary()) > Config.SPECIAL_FUNCTIONS_TOLERANCE) {
        if (i++ > iterationLimit && iterationLimit > 0) {
            IterationLimitExceeded.throwIt(i, org.matheclipse.core.expression.S.HurwitzZeta);
        }
        if (i > MAX_VALUE_HALF) {
            throw new ArgumentTypeException("Hurwitz zeta: i > MAX_VALUE_HALF");
        }
        int iPlusi = i + i;
        p = p.multiply(x.add(iPlusi - 2.0).multiply(x.add(iPlusi - 3.0)).multiply(a.add(n).pow(2.0).multiply(iPlusi * (iPlusi - 1)).reciprocal()));
        t = t.add(p.multiply(bernoulliInt(iPlusi)));
    }
    Complex T = t.add(0.5).divide(a.add(n).pow(x));
    return S.add(I).add(T);
}

Example 3 with S

use of org.matheclipse.core.expression.S in project symja_android_library by axkr.

the class Pods method parseInput.

/**
 * package private
 */
static IExpr parseInput(String inputStr, EvalEngine engine) {
    engine.setPackageMode(false);
    IExpr inExpr = F.NIL;
    // }
    if (!inExpr.isPresent()) {
        final FuzzyParser parser = new FuzzyParser(engine);
        try {
            inExpr = parser.parseFuzzyList(inputStr);
        } catch (SyntaxError serr) {
            // this includes syntax errors
            LOGGER.debug("Pods: FuzzyParser.parseFuzzyList() failed", serr);
            TeXParser texConverter = new TeXParser(engine);
            inExpr = texConverter.toExpression(inputStr);
        }
    }
    if (inExpr == S.$Aborted) {
        return F.NIL;
    }
    if (inExpr.isList() && inExpr.size() == 2) {
        inExpr = inExpr.first();
    }
    if (inExpr.isTimes() && !inExpr.isNumericFunction(true) && inExpr.argSize() <= 4) {
        if (((IAST) inExpr).isEvalFlagOn(IAST.TIMES_PARSED_IMPLICIT)) {
            inExpr = flattenTimes((IAST) inExpr).orElse(inExpr);
            IAST rest = ((IAST) inExpr).setAtClone(0, S.List);
            IASTAppendable specialFunction = F.NIL;
            String stemForm = getStemForm(rest.arg1().toString().toLowerCase());
            IExpr head = rest.head();
            if (stemForm != null) {
                head = STEM.getSymbol(stemForm);
                if (head != null) {
                    specialFunction = rest.setAtClone(0, head);
                    specialFunction.remove(1);
                }
            }
            if (!specialFunction.isPresent()) {
                stemForm = getStemForm(rest.last().toString().toLowerCase());
                if (stemForm != null) {
                    head = STEM.getSymbol(stemForm);
                    if (head != null) {
                        specialFunction = rest.setAtClone(0, head);
                        specialFunction.remove(rest.size() - 1);
                    }
                }
            }
            if (specialFunction.isPresent()) {
                if (head != null) {
                    if (head == S.UnitConvert) {
                        IExpr temp = unitConvert(engine, rest.rest());
                        if (temp.isPresent()) {
                            return temp;
                        }
                    } else {
                        int i = 1;
                        while (i < specialFunction.size()) {
                            String argStr = specialFunction.get(i).toString().toLowerCase();
                            if (argStr.equalsIgnoreCase("by") || argStr.equalsIgnoreCase("for")) {
                                specialFunction.remove(i);
                                continue;
                            }
                            i++;
                        }
                        return specialFunction;
                    }
                }
            }
            if (rest.arg1().toString().equalsIgnoreCase("convert")) {
                rest = inExpr.rest();
            }
            if (rest.argSize() > 2) {
                rest = rest.removeIf(x -> x.toString().equals("in"));
            }
            IExpr temp = unitConvert(engine, rest);
            if (temp.isPresent()) {
                return temp;
            }
        }
    }
    return inExpr;
}

Example 4 with S

use of org.matheclipse.core.expression.S in project symja_android_library by axkr.

the class MatrixD method evaluate.

/**
 * For the referenced formula numbers (XX) see:
 * <a href="https://archive.org/details/imm3274/">Internet Archive - The Matrix Cookbook</a>
 */
@Override
public IExpr evaluate(final IAST ast, EvalEngine engine) {
    IOFunctions.printExperimental(S.MatrixD);
    if (ast.size() < 3) {
        return F.NIL;
    }
    final IExpr fx = ast.arg1();
    if (fx.isIndeterminate()) {
        return S.Indeterminate;
    }
    if (ast.size() > 3) {
        // reduce arguments by folding MatrixD[fxy, x, y] to MatrixD[ MatrixD[fxy, x], y] ...
        return ast.foldLeft((x, y) -> engine.evaluateNIL(F.MatrixD(x, y)), fx, 2);
    }
    IExpr x = ast.arg2();
    if (!(x.isVariable() || x.isList())) {
        // `1` is not a valid variable.
        return IOFunctions.printMessage(ast.topHead(), "ivar", F.list(x), engine);
    }
    if (fx.isList()) {
        IAST list = (IAST) fx;
        // thread over first list
        return list.mapThreadEvaled(engine, F.ListAlloc(list.size()), ast, 1);
    }
    if (x.isList()) {
        // MatrixD[fx_, {...}]
        IAST xList = (IAST) x;
        if (xList.isAST1() && xList.arg1().isListOfLists()) {
            IAST subList = (IAST) xList.arg1();
            IASTAppendable result = F.ListAlloc(subList.size());
            result.appendArgs(subList.size(), i -> F.MatrixD(fx, F.list(subList.get(i))));
            return result;
        } else if (xList.isAST1() && xList.arg1().isList()) {
            IAST subList = (IAST) xList.arg1();
            return subList.mapLeft(F.ListAlloc(), (a, b) -> engine.evaluateNIL(F.MatrixD(a, b)), fx);
        } else if (xList.isAST2()) {
            if (xList.arg1().isList()) {
                x = F.list(xList.arg1());
            } else {
                x = xList.arg1();
            }
            IExpr arg2 = xList.arg2();
            int n = arg2.toIntDefault();
            if (n >= 0) {
                IExpr temp = fx;
                for (int i = 0; i < n; i++) {
                    temp = S.MatrixD.ofNIL(engine, temp, x);
                    if (!temp.isPresent()) {
                        return F.NIL;
                    }
                }
                return temp;
            }
            if (arg2.isFree(num -> num.isNumber(), false)) {
                if (fx.equals(x)) {
                    return S.$SingleEntryMatrix;
                }
                if (fx.isAST()) {
                    final IAST function = (IAST) fx;
                    if (function.isPlus()) {
                        // (35)
                        return function.mapThread(F.MatrixD(F.Slot1, xList), 1);
                    }
                }
                return F.NIL;
            }
            if (!x.isVariable()) {
                // `1` is not a valid variable.
                return IOFunctions.printMessage(ast.topHead(), "ivar", F.list(x), engine);
            }
            if (arg2.isAST()) {
                return F.NIL;
            }
            // symbolic expression or a non-negative integer.
            return IOFunctions.printMessage(ast.topHead(), "dvar", F.list(xList), engine);
        }
        return F.NIL;
    }
    if (!x.isVariable()) {
        // `1` is not a valid variable.
        return IOFunctions.printMessage(ast.topHead(), "ivar", F.list(x), engine);
    }
    return binaryMatrixD(fx, x);
}

Example 5 with S

use of org.matheclipse.core.expression.S in project symja_android_library by axkr.

the class HypergeometricJS method hypergeometric1F2.

public static Complex hypergeometric1F2(Complex a, Complex b, Complex c, Complex x) {
    final int useAsymptotic = 200;
    if (x.norm() > useAsymptotic) {
        Complex p = a.add(b.negate()).add(c.negate()).add(0.5).divide(2.0);
        ArrayList<Complex> ck = new ArrayList<Complex>();
        // 
        ck.add(Complex.ONE);
        ck.add(((a.multiply(3.0).add(b).add(c).add(-2.0)).multiply(a.subtract(b.add(c))).multiply(0.5)).add(b.multiply(c).multiply(2)).add(// 
        -3.0 / 8.0));
        ck.add((a.multiply(3.0).add(b).add(c).add(-2.0)).multiply(a.subtract(b.add(c))).multiply(0.25).add(b.multiply(c).add(-3.0 / 16.0)).pow(2).multiply(// 
        2));
        // 
        ck.add(new Complex(-1.0).multiply(a.multiply(2.0).subtract(3.0)).multiply(b).multiply(c));
        ck.add(a.pow(2.0).multiply(-8.0).add(a.multiply(11.0)).add(b).add(c).add(-2.0).multiply(a.subtract(b.add(c))).multiply(// 
        0.25));
        ck.add(new Complex(-3.0 / 16.0));
        IntFunction<Complex> w = k -> ck.get(k).multiply(x.negate().pow(-k / 2.0)).divide(Math.pow(2.0, k));
        Complex u1 = Complex.I.multiply(p.multiply(Math.PI).add(x.negate().sqrt().multiply(2.0))).exp();
        Complex u2 = new Complex(0.0, -1.0).multiply(p.multiply(Math.PI).add(x.negate().sqrt().multiply(2.0))).exp();
        Complex wLast = w.apply(2);
        Complex w2Negate = wLast.negate();
        Complex s = // 
        u1.multiply(new Complex(0.0, -1.0).multiply(w.apply(1)).add(w2Negate).add(1.0)).add(u2.multiply(Complex.I.multiply(w.apply(1)).add(w2Negate).add(1.0)));
        int k = 3;
        while (wLast.norm() > w.apply(k).norm()) {
            // 
            ck.add(a.multiply(-6.0).add(b.multiply(2)).add(c.multiply(2.0)).add(-4.0).multiply(k).add(a.pow(a).multiply(3.0)).add(b.subtract(c).pow(2.0).negate()).add(a.multiply(b.add(c).add(-2)).multiply(2.0).negate()).add(0.25).add(3.0 * k * k).multiply(1.0 / (2.0 * k)).multiply(// 
            ck.get(k - 1)).subtract(a.negate().add(b).add(c.negate()).add(-0.5).add(k).multiply(a.negate().add(b.negate()).add(c).add(-0.5).add(k)).multiply(a.negate().add(b).add(c).add(-2.5).add(k)).multiply(// 
            ck.get(k - 2))));
            wLast = w.apply(k);
            s = s.add(// 
            u1.multiply(new Complex(0.0, -1.0).pow(k)).multiply(wLast).add(u2.multiply(Complex.I.pow(k)).multiply(wLast)));
            k++;
        }
        Complex t1 = Arithmetic.lanczosApproxGamma(a).reciprocal().multiply(x.negate().pow(p)).multiply(s).divide(2.0 * Math.sqrt(Math.PI));
        Complex t2 = Arithmetic.lanczosApproxGamma(b.subtract(a)).reciprocal().multiply(Arithmetic.lanczosApproxGamma(c.subtract(a)).reciprocal()).multiply(x.negate().pow(a.negate())).multiply(hypergeometricSeries(new Complex[] { a, a.add(b.negate()).add(1), a.add(c.negate().add(1.0)) }, new Complex[] {}, // , true ) );
        x.reciprocal()));
        // hypergeometricSeries( [ a, add(a,neg(b),1), add(a,neg(c),1) ], [], inv(x), true ) );
        return Arithmetic.lanczosApproxGamma(b).multiply(Arithmetic.lanczosApproxGamma(c)).multiply(t1.add(t2));
    }
    return hypergeometricSeries(new Complex[] { a }, new Complex[] { b, c }, x);
}