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);
}
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);
}
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);
}
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);
}
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();
}
}
}
Aggregations