Examples with Complex org.hipparchus.complex.Complex used on opensource projects

Example 56 with Complex

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

the class Fourier method fourier.

/**
 * Uses decimation-in-time or Cooley-Tukey FFT
 *
 * @param vector of length of power of 2
 * @param b is +1 for forward, and -1 for inverse transform
 * @return discrete Fourier transform of given vector
 */
private static IAST fourier(IAST vector, int b) {
    final int n = vector.argSize();
    Complex[] array = Convert.list2Complex(vector);
    if (array == null) {
        return F.NIL;
    }
    // FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.UNITARY);
    // org.hipparchus.complex.Complex[] result = fft.transform(array, TransformType.FORWARD);
    // return Object2Expr.convertComplex(true, result);
    int j = 0;
    for (int i = 0; i < n; ++i) {
        if (j > i) {
            Complex val = array[i];
            array[i] = array[j];
            array[j] = val;
        }
        int m = n >> 1;
        while (m > 0 && j >= m) {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    int mmax = 1;
    while (n > mmax) {
        int istep = mmax << 1;
        final double thalf = b * Math.PI / istep;
        final double wtemp = Math.sin(thalf);
        Complex wp = new Complex(-2 * wtemp * wtemp, Math.sin(thalf + thalf));
        Complex w = Complex.ONE;
        for (int m = 0; m < mmax; ++m) {
            for (int i = m; i < n; i += istep) {
                j = i + mmax;
                Complex temp = array[j].multiply(w);
                array[j] = array[i].subtract(temp);
                array[i] = array[i].add(temp);
            }
            w = w.add(w.multiply(wp));
        }
        mmax = istep;
    }
    return F.Divide(Convert.toVector(array), F.Sqrt(n));
}

Example 57 with Complex

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

the class ND method partialDerivative.

private IExpr partialDerivative(IExpr function, ISymbol variable, int order, IExpr value, EvalEngine engine) {
    double a3Double = Double.NaN;
    try {
        a3Double = value.evalDouble();
    } catch (ValidateException ve) {
    }
    if (Double.isNaN(a3Double)) {
        Complex a3Complex = Complex.NaN;
        a3Complex = value.evalComplex();
        if (a3Complex != null) {
        // FDSFactory<Complex> factory = new FDSFactory<Complex>(ComplexField.getInstance(),
        // 1, order);
        // // FieldDerivativeStructure<Complex> f = factory.variable(0, a3Complex);
        // FiniteDifferencesDifferentiator differentiator =
        // new FiniteDifferencesDifferentiator(15, 0.01);
        // UnivariateDifferentiableFunction f =
        // differentiator.differentiate(new UnaryNumerical(arg1, arg2,
        // EvalEngine.get()));
        // // return F.complexNum(f.getPartialDerivative(order));
        }
    } else {
        DSFactory factory = new DSFactory(1, order);
        FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(15, 0.01);
        UnivariateDifferentiableFunction f = differentiator.differentiate(new UnaryNumerical(function, variable, EvalEngine.get()));
        return F.num(f.value(factory.variable(0, a3Double)).getPartialDerivative(order));
    }
    return F.NIL;
}

Example 58 with Complex

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

the class ManipulateFunction method complexRange.

/**
 * Convert a range of 2 complex numbers into a <code>ComplexPlot3D</code> compatible range.
 *
 * @param graphicControl
 * @param plotRange example <code>{z, -2-2*I, 2+2*I}</code>
 * @param steps an additional step parameter. If less <code>0</code> the parameter will be ignored
 * @param toJS the expression to JavaScript transpiler
 * @return an array for the x-(real)-range and the y-(imaginary)-range of the 3D plot
 */
private static double[] complexRange(StringBuilder graphicControl, IAST plotRange, int steps, JavaScriptFormFactory toJS) {
    double[] result = new double[2];
    IExpr zMin = plotRange.arg2();
    IExpr zMax = plotRange.arg3();
    Complex cMin = zMin.evalComplex();
    Complex cMax = zMax.evalComplex();
    if (cMin != null && cMax != null) {
        double reMin = cMin.getReal();
        double imMin = cMin.getImaginary();
        double reMax = cMax.getReal();
        double imMax = cMax.getImaginary();
        if (reMin > reMax) {
            // swap range values
            double r = reMax;
            reMax = reMin;
            reMin = r;
        }
        if (imMin > imMax) {
            // swap range values
            double i = imMax;
            imMax = imMin;
            imMin = i;
        }
        result[0] = reMax - reMin;
        result[1] = imMax - imMin;
        graphicControl.append(", [");
        toJS.convert(graphicControl, F.num(reMin));
        graphicControl.append(", ");
        toJS.convert(graphicControl, F.num(reMax));
        if (steps > 0) {
            graphicControl.append(", ");
            graphicControl.append(steps);
        }
        graphicControl.append("], [");
        toJS.convert(graphicControl, F.num(imMin));
        graphicControl.append(", ");
        toJS.convert(graphicControl, F.num(imMax));
        if (steps > 0) {
            graphicControl.append(", ");
            graphicControl.append(steps);
        }
        graphicControl.append("]");
    }
    return result;
}

Example 59 with Complex

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

the class EllipticFunctionsJS method jacobiTheta.

public static Complex jacobiTheta(int n, Complex x, Complex q, double tolerance) {
    if (q.norm() >= 1) {
        throw new ArgumentTypeException("unsupported elliptic nome");
    }
    if (n < 1 || n > 4) {
        throw new ArgumentTypeException("undefined Jacobi theta index");
    }
    if (F.isZero(q)) {
        switch(n) {
            case 1:
            case 2:
                return Complex.ZERO;
            case 3:
            case 4:
                return Complex.ONE;
        }
    }
    Complex piTau = q.log().divide(Complex.I);
    // dlmf.nist.gov/20.2 to reduce overflow
    if (Math.abs(x.getImaginary()) > Math.abs(piTau.getImaginary()) || Math.abs(x.getReal()) > Math.PI) {
        double pt = trunc(x.getImaginary() / piTau.getImaginary());
        x = x.subtract(piTau.multiply(pt));
        double p = trunc(x.getReal() / Math.PI);
        x = x.subtract(p * Math.PI);
        Complex qFactor = q.pow(-pt * pt);
        Complex eFactor = x.multiply(Complex.I).multiply(-2 * pt).exp();
        // factors can become huge, so chop spurious parts first
        switch(n) {
            case 1:
                return qFactor.multiply(eFactor).multiply(F.chopComplex(jacobiTheta(n, x, q), tolerance)).multiply(Math.pow((-1), (p + pt)));
            case 2:
                return qFactor.multiply(eFactor).multiply(F.chopComplex(jacobiTheta(n, x, q), tolerance)).multiply(Math.pow((-1), p));
            case 3:
                return qFactor.multiply(eFactor).multiply(F.chopComplex(jacobiTheta(n, x, q), tolerance));
            case 4:
                return qFactor.multiply(eFactor).multiply(F.chopComplex(jacobiTheta(n, x, q), tolerance)).multiply(Math.pow((-1), pt));
        }
    }
    Complex s = Complex.ZERO;
    Complex p = Complex.ONE;
    long iterationLimit = EvalEngine.get().getIterationLimit();
    int i = 0;
    switch(n) {
        case 1:
            while (Math.abs(p.getReal()) > tolerance || Math.abs(p.getImaginary()) > tolerance) {
                p = q.pow(i * i + i).multiply(x.multiply(2 * i + 1).sin()).multiply(Math.pow(-1, i));
                s = s.add(p);
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return q.pow(0.25).multiply(s).multiply(2);
        case 2:
            while (Math.abs(p.getReal()) > tolerance || Math.abs(p.getImaginary()) > tolerance) {
                p = q.pow(i * i + i).multiply(x.multiply(2 * i + 1).cos());
                s = s.add(p);
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return q.pow(0.25).multiply(s).multiply(2);
        case 3:
            i = 1;
            while (Math.abs(p.getReal()) > tolerance || Math.abs(p.getImaginary()) > tolerance) {
                p = q.pow(i * i).multiply(x.multiply(2 * i).cos());
                s = s.add(p);
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return s.multiply(2.0).add(1.0);
        case 4:
            i = 1;
            while (Math.abs(p.getReal()) > tolerance || Math.abs(p.getImaginary()) > tolerance) {
                p = q.negate().pow(i * i).multiply(x.multiply(2 * i).cos());
                s = s.add(p);
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return s.multiply(2.0).add(1.0);
    }
    throw new ArgumentTypeException("undefined Jacobi theta index");
}

Example 60 with Complex

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

the class EllipticFunctionsJS method jacobiTheta.

public static Complex jacobiTheta(int n, double x, double q, double tolerance) {
    if (Math.abs(q) >= 1) {
        throw new ArgumentTypeException("unsupported elliptic nome");
    }
    if (n < 1 || n > 4) {
        throw new ArgumentTypeException("undefined Jacobi theta index");
    }
    if (F.isZero(q)) {
        switch(n) {
            case 1:
            case 2:
                return Complex.ZERO;
            case 3:
            case 4:
                return Complex.ONE;
        }
    }
    // dlmf.nist.gov/20.2 to reduce overflow
    if (Math.abs(x) > Math.PI) {
        double p = trunc(x / Math.PI);
        x = x - p * Math.PI;
        switch(n) {
            case 1:
            case 2:
                return new Complex(Math.pow(-1, p)).multiply(jacobiTheta(n, x, q));
            case 3:
            case 4:
                return jacobiTheta(n, x, q);
        }
    }
    long iterationLimit = EvalEngine.get().getIterationLimit();
    switch(n) {
        case 1:
            if (q < 0) {
                return jacobiTheta(n, new Complex(x), new Complex(q));
            }
            double s = 0;
            double p = 1;
            int i = 0;
            while (Math.abs(p) > tolerance) {
                p = Math.pow(-1, i) * Math.pow(q, (i * i + i)) * Math.sin((2 * i + 1) * x);
                s += p;
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return new Complex(2 * Math.pow(q, 0.25) * s);
        case 2:
            if (q < 0) {
                return jacobiTheta(n, new Complex(x), new Complex(q));
            }
            s = 0;
            p = 1;
            i = 0;
            while (Math.abs(p) > tolerance) {
                p = Math.pow(q, (i * i + i)) * Math.cos((2 * i + 1) * x);
                s += p;
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return new Complex(2 * Math.pow(q, 0.25) * s);
        case 3:
            s = 0;
            p = 1;
            i = 1;
            while (Math.abs(p) > tolerance) {
                p = Math.pow(q, (i * i)) * Math.cos(2 * i * x);
                s += p;
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return new Complex(1 + 2 * s);
        case 4:
            s = 0;
            p = 1;
            i = 1;
            while (Math.abs(p) > tolerance) {
                p = Math.pow(-q, (i * i)) * Math.cos(2 * i * x);
                s += p;
                if (i++ > iterationLimit && iterationLimit > 0) {
                    IterationLimitExceeded.throwIt(i, S.EllipticTheta);
                }
            }
            return new Complex(1 + 2 * s);
    }
    throw new ArgumentTypeException("undefined Jacobi theta index");
}