Examples with Function de.lab4inf.math.Function used on opensource projects

Example 1 with Function

use of de.lab4inf.math.Function 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));
}

Example 2 with Function

use of de.lab4inf.math.Function in project symja_android_library by axkr.

the class Integrator method integrate.

/**
 * Public method to integrate the given function over the interval [a,b]. The special cases where
 * one or both borders are infinite are handled by this method and transformed to a proper
 * integration kernel.
 *
 * @param a double lower integration border
 * @param b double upper integration border
 * @param eps double the relative tolerated error
 * @param fct Function to integrate
 * @return double the result of the integration
 */
public static double integrate(final double a, final double b, final double eps, final Function fct) {
    Function kernel = fct;
    double left = a, right = b;
    // the [0,1] or [-1,1] or [0,1/z] interval with a new kernel...
    if (isInfinite(a) || isInfinite(b)) {
        if ((a == 0 || b == 0) || (isInfinite(a) && isInfinite(b))) {
            kernel = new KernelZeroToInfinity(fct);
            if (isInfinite(a)) {
                left = -1;
            }
            if (isInfinite(b)) {
                right = 1;
            }
        } else {
            kernel = new KernelInfinity(fct);
            if (isInfinite(a)) {
                right = 0;
                left = 1.0 / b;
            } else {
                right = 1.0 / a;
                left = 0;
            }
        }
    }
    return integrateSelected(left, right, eps, kernel);
}

Example 3 with Function

use of de.lab4inf.math.Function in project symja_android_library by axkr.

the class SecondOrderSolver method solve.

/* (non-Javadoc)
   * @see de.lab4inf.math.ode.SecondOrderOdeSolver#solve()
   */
public double solve(final double x0, final double y0, final double dy0, final double x1, final Function f, final double eps) {
    // solve y'' = f(x,y,y') by setting
    // y' = z
    // z' = f(x,y,y')
    Function z = new Ode2Wrapper();
    Function[] sys = { z, f };
    double[] ybeg = { y0, dy0 };
    double[] yend = solve(x0, ybeg, x1, sys, eps);
    return yend[0];
}

Example 4 with Function

use of de.lab4inf.math.Function in project symja_android_library by axkr.

the class NewtonRootFinder method modnewton.

/**
 * Modified Newton method for multiple roots.
 *
 * @param fct the function to find the root of
 * @param dF the first derivative of the function
 * @param x0 the initial guess as starting value
 * @param eps the accuracy to reach
 * @return the root
 */
public static double modnewton(final Function fct, final Function dF, final double x0, final double eps) {
    Function d2F;
    if (dF instanceof Differentiable) {
        d2F = ((Differentiable) dF).getDerivative();
    } else {
        getLogger().warn(NO_DERIVATIVE_NOT_RECOMMENDED);
        d2F = new Differentiator(dF);
    }
    return modnewton(fct, dF, d2F, x0, eps);
}

Example 5 with Function

use of de.lab4inf.math.Function in project symja_android_library by axkr.

the class ProbabilityDistribution method quantile.

/**
 * Internal helper function to find the quantile xp with cdf(xp)=p for a given p value, via a
 * Newton method.
 *
 * @param x0 the initial starting value
 * @param err the error(x) = cdf(x) - p
 * @return xp with err(xp)=0
 */
private static double quantile(final double x0, final Differentiable err) {
    final Function pdf = err.getDerivative();
    double g, e, t, f, d, x, y = x0;
    int n = 0;
    do {
        g = -1;
        x = y;
        e = err.f(x);
        if (e == 0) {
            return x;
        }
        d = e / max(2 * abs(e / x), pdf.f(x));
        t = y + g * d;
        try {
            f = err.f(t);
            while (abs(f) > abs(e) && g <= 1) {
                g += .33;
                t = y + g * d;
                f = err.f(t);
            // logger.warn("shrinking dx: " + g);
            }
        } catch (final Throwable error) {
            getLogger().warn("error for " + t);
            if (t > 1) {
                g = (y - 1) / (2 * d);
            }
            t = y + g * d;
            getLogger().warn("using " + t);
        }
        if (y + g * d < 0) {
            getLogger().warn("dx too large: " + g * d);
            d = -x / (2 * g);
        }
        y = x + g * d;
    } while (!hasReachedAccuracy(y, x, PRECISION) && ++n < MAX_ITE);
    if (n >= MAX_ITE) {
        final double p = -err.f(0);
        final String msg = String.format("quantile(%f)=%f no convergence", p, x);
        getLogger().warn(msg);
    }
    return y;
}