Skip to content

Commit

Permalink
use apfloat erf, erfc, fresnelC, fresnelS for apfloat nums
Browse files Browse the repository at this point in the history 
TODO:
- mtommila/apfloat#31 (comment)
  • Loading branch information
@axkr
axkr committed Jan 14, 2024
1 parent 4c4d262 commit c94a4d0
Show file tree
Hide file tree
Showing 8 changed files with 225 additions and 135 deletions.
16 changes: 13 additions & 3 deletions .../matheclipse-core/src/main/java/org/matheclipse/core/builtin/HypergeometricFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -999,8 +999,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
}
if (engine.isDoubleMode()) {
} else if (engine.isDoubleMode()) {

double aDouble = Double.NaN;
double bDouble = Double.NaN;
Expand Down Expand Up @@ -1407,7 +1406,18 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
0, nMinus1));
}
}
if (engine.isDoubleMode()) {
if (engine.isArbitraryMode()) {
try {
IExpr res = a.hypergeometric1F1(b, z);
if (res.isNumber()) {
return res;
}
} catch (ValidateException ve) {
return Errors.printMessage(ast.topHead(), ve, engine);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
} else if (engine.isDoubleMode()) {
double aDouble = Double.NaN;
double bDouble = Double.NaN;
double zDouble = Double.NaN;
Expand Down
9 changes: 5 additions & 4 deletions ...ibrary/matheclipse-core/src/main/java/org/matheclipse/core/builtin/functions/GammaJS.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@
import org.matheclipse.core.builtin.ConstantDefinitions;
import org.matheclipse.core.eval.EvalEngine;
import org.matheclipse.core.eval.exception.ArgumentTypeException;
import org.matheclipse.core.expression.ApcomplexNum;
import org.matheclipse.core.expression.F;
import org.matheclipse.core.expression.NumberUtil;
import org.matheclipse.core.generic.UnaryNumerical;
Expand Down Expand Up @@ -227,13 +226,15 @@ public static INumber incompleteBeta(double x, double y, double z) {

public static Complex fresnelS(Complex x) {
Apcomplex apcomplex = new Apcomplex(new Apfloat(x.getReal()), new Apfloat(x.getImaginary()));
Apcomplex fresnelC = ApcomplexNum.fresnelS(apcomplex, EvalEngine.getApfloatDouble());
return new Complex(fresnelC.real().doubleValue(), fresnelC.imag().doubleValue());
Apcomplex fresnelS = EvalEngine.getApfloat().fresnelS(apcomplex);
// Apcomplex fresnelS = ApcomplexNum.fresnelS(apcomplex, EvalEngine.getApfloatDouble());
return new Complex(fresnelS.real().doubleValue(), fresnelS.imag().doubleValue());
}

public static Complex fresnelC(Complex x) {
Apcomplex apcomplex = new Apcomplex(new Apfloat(x.getReal()), new Apfloat(x.getImaginary()));
Apcomplex fresnelC = ApcomplexNum.fresnelC(apcomplex, EvalEngine.getApfloatDouble());
Apcomplex fresnelC = EvalEngine.getApfloat().fresnelC(apcomplex);
// Apcomplex fresnelC = ApcomplexNum.fresnelC(apcomplex, EvalEngine.getApfloatDouble());
return new Complex(fresnelC.real().doubleValue(), fresnelC.imag().doubleValue());
}

Expand Down
174 changes: 93 additions & 81 deletions ..._library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ApcomplexNum.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,6 @@
import org.apfloat.Apfloat;
import org.apfloat.ApfloatMath;
import org.apfloat.ApfloatRuntimeException;
import org.apfloat.Apint;
import org.apfloat.Aprational;
import org.apfloat.FixedPrecisionApfloatHelper;
import org.hipparchus.complex.Complex;
import org.matheclipse.core.basic.Config;
Expand Down Expand Up @@ -218,63 +216,65 @@ public IExpr hypergeometric0F1(IExpr arg2) {

@Override
public IExpr fresnelC() {
return valueOf(fresnelC(fApcomplex, EvalEngine.getApfloat()));
}

public static Apcomplex fresnelC(Apcomplex z, FixedPrecisionApfloatHelper apfloat) {
// Complex m1 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// new Complex(0, Math.PI / 2).multiply(x.multiply(x)));
// Complex m2 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// new Complex(0, -Math.PI / 2).multiply(x.multiply(x)));
//
// return mul(0.5, x, m1.add(m2));
// z*z
Apcomplex sqr = apfloat.multiply(z, z);
// Math.PI / 2
Apfloat piHalf = apfloat.divide(apfloat.pi(), new Apint(2));
// 1/2
Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// 3/2
Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
Apcomplex m1 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
new Apcomplex(Apfloat.ZERO, piHalf).multiply(sqr));
Apcomplex m2 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
new Apcomplex(Apfloat.ZERO, piHalf.negate()).multiply(sqr));

return apfloat.multiply(oneHalf, z).multiply(apfloat.add(m1, m2));
}
return valueOf(EvalEngine.getApfloat().fresnelC(fApcomplex));
// return valueOf(fresnelC(fApcomplex, EvalEngine.getApfloat()));
}

// public static Apcomplex fresnelC(Apcomplex z, FixedPrecisionApfloatHelper apfloat) {
// // Complex m1 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// // new Complex(0, Math.PI / 2).multiply(x.multiply(x)));
// // Complex m2 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// // new Complex(0, -Math.PI / 2).multiply(x.multiply(x)));
// //
// // return mul(0.5, x, m1.add(m2));
// // z*z
// Apcomplex sqr = apfloat.multiply(z, z);
// // Math.PI / 2
// Apfloat piHalf = apfloat.divide(apfloat.pi(), new Apint(2));
// // 1/2
// Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// // 3/2
// Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
// Apcomplex m1 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
// new Apcomplex(Apfloat.ZERO, piHalf).multiply(sqr));
// Apcomplex m2 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
// new Apcomplex(Apfloat.ZERO, piHalf.negate()).multiply(sqr));
//
// return apfloat.multiply(oneHalf, z).multiply(apfloat.add(m1, m2));
// }

@Override
public IExpr fresnelS() {
return valueOf(fresnelS(fApcomplex, EvalEngine.getApfloat()));
}

public static Apcomplex fresnelS(Apcomplex z, FixedPrecisionApfloatHelper apfloat) {
// Complex m1 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// new Complex(0, Math.PI / 2).multiply(x.multiply(x)));
// Complex m2 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// new Complex(0, -Math.PI / 2).multiply(x.multiply(x)));
//
// return mul(new Complex(0, -0.5), x, sub(m1, m2));

// z*z
Apcomplex sqr = apfloat.multiply(z, z);
// Math.PI / 2
Apfloat piHalf = apfloat.divide(apfloat.pi(), new Apint(2));
// 1/2
Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// -1/2
Aprational minusOneHalf = new Aprational(new Apint(-1), new Apint(2));
// 3/2
Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
Apcomplex m1 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
new Apcomplex(Apfloat.ZERO, piHalf).multiply(sqr));
Apcomplex m2 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
new Apcomplex(Apfloat.ZERO, piHalf.negate()).multiply(sqr));

return apfloat.multiply(new Apcomplex(Apfloat.ZERO, minusOneHalf), z)
.multiply(apfloat.subtract(m1, m2));
}
return valueOf(EvalEngine.getApfloat().fresnelS(fApcomplex));
// return valueOf(fresnelS(fApcomplex, EvalEngine.getApfloat()));
}

// public static Apcomplex fresnelS(Apcomplex z, FixedPrecisionApfloatHelper apfloat) {
// // Complex m1 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// // new Complex(0, Math.PI / 2).multiply(x.multiply(x)));
// // Complex m2 = HypergeometricJS.hypergeometric1F1(new Complex(0.5), new Complex(1.5),
// // new Complex(0, -Math.PI / 2).multiply(x.multiply(x)));
// //
// // return mul(new Complex(0, -0.5), x, sub(m1, m2));
//
// // z*z
// Apcomplex sqr = apfloat.multiply(z, z);
// // Math.PI / 2
// Apfloat piHalf = apfloat.divide(apfloat.pi(), new Apint(2));
// // 1/2
// Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// // -1/2
// Aprational minusOneHalf = new Aprational(new Apint(-1), new Apint(2));
// // 3/2
// Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
// Apcomplex m1 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
// new Apcomplex(Apfloat.ZERO, piHalf).multiply(sqr));
// Apcomplex m2 = apfloat.hypergeometric1F1(oneHalf, threeHalf,
// new Apcomplex(Apfloat.ZERO, piHalf.negate()).multiply(sqr));
//
// return apfloat.multiply(new Apcomplex(Apfloat.ZERO, minusOneHalf), z)
// .multiply(apfloat.subtract(m1, m2));
// }

@Override
public IExpr hypergeometric1F1(IExpr arg2, IExpr arg3) {
Expand All @@ -295,6 +295,15 @@ public IExpr hypergeometric2F1(IExpr arg2, IExpr arg3, IExpr arg4) {
return IComplexNum.super.hypergeometric2F1(arg2, arg3, arg4);
}

@Override
public IExpr hypergeometricU(IExpr arg2, IExpr arg3) {
if (arg2 instanceof INumber && arg3 instanceof INumber) {
return valueOf(EvalEngine.getApfloat().hypergeometric1F1(fApcomplex,
((INumber) arg2).apcomplexValue(), ((INumber) arg3).apcomplexValue()));
}
return IComplexNum.super.hypergeometricU(arg2, arg3);
}

@Override
public IComplexNum add(final IComplexNum val) {
return valueOf(EvalEngine.getApfloat().add(fApcomplex, val.apcomplexValue()));
Expand Down Expand Up @@ -397,36 +406,39 @@ public IExpr evaluate(EvalEngine engine) {

@Override
public IExpr erf() {
FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
try {
Apcomplex erf = erf(fApcomplex, h);
return F.complexNum(erf);
} catch (Exception ce) {
//
}
return F.NIL;
return valueOf(EvalEngine.getApfloat().erf(fApcomplex));
// FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
// try {
// Apcomplex erf = erf(fApcomplex, h);
// return F.complexNum(erf);
// } catch (Exception ce) {
// //
// }
// return F.NIL;
}

private static Apcomplex erf(Apcomplex x, FixedPrecisionApfloatHelper h) {
Apint two = new Apint(2);
Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// 3/2
Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
Apcomplex erf = h.hypergeometric1F1(oneHalf, threeHalf, h.multiply(x, x).negate()).multiply(two)
.multiply(x).divide(h.sqrt(h.pi()));
return erf;
}
// private static Apcomplex erf(Apcomplex x, FixedPrecisionApfloatHelper h) {
// Apint two = new Apint(2);
// Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// // 3/2
// Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
// Apcomplex erf = h.hypergeometric1F1(oneHalf, threeHalf, h.multiply(x,
// x).negate()).multiply(two)
// .multiply(x).divide(h.sqrt(h.pi()));
// return erf;
// }

@Override
public IExpr erfc() {
FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
try {
Apcomplex c = erf(fApcomplex, h);
return F.complexNum(h.subtract(Apcomplex.ONE, c));
} catch (Exception ce) {
//
}
return F.NIL;
return valueOf(EvalEngine.getApfloat().erfc(fApcomplex));
// FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
// try {
// Apcomplex c = erf(fApcomplex, h);
// return F.complexNum(h.subtract(Apcomplex.ONE, c));
// } catch (Exception ce) {
// //
// }
// return F.NIL;
}

@Override
Expand Down
82 changes: 51 additions & 31 deletions ...id_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ApfloatNum.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
import org.apfloat.ApfloatMath;
import org.apfloat.ApfloatRuntimeException;
import org.apfloat.Apint;
import org.apfloat.Aprational;
import org.apfloat.FixedPrecisionApfloatHelper;
import org.matheclipse.core.basic.Config;
import org.matheclipse.core.eval.EvalEngine;
Expand Down Expand Up @@ -87,12 +86,14 @@ public int hierarchy() {

@Override
public IExpr fresnelC() {
return ApcomplexNum.valueOf(ApcomplexNum.fresnelC(fApfloat, EvalEngine.getApfloat()));
return valueOf(EvalEngine.getApfloat().fresnelC(fApfloat));
// return ApcomplexNum.valueOf(ApcomplexNum.fresnelC(fApfloat, EvalEngine.getApfloat()));
}

@Override
public IExpr fresnelS() {
return ApcomplexNum.valueOf(ApcomplexNum.fresnelS(fApfloat, EvalEngine.getApfloat()));
return valueOf(EvalEngine.getApfloat().fresnelS(fApfloat));
// return ApcomplexNum.valueOf(ApcomplexNum.fresnelS(fApfloat, EvalEngine.getApfloat()));
}

@Override
Expand Down Expand Up @@ -156,6 +157,23 @@ public IExpr hypergeometric2F1(IExpr arg2, IExpr arg3, IExpr arg4) {
return INum.super.hypergeometric2F1(arg2, arg3, arg4);
}

@Override
public IExpr hypergeometricU(IExpr arg2, IExpr arg3) {
if (arg2 instanceof IReal && arg3 instanceof IReal) {
try {
return valueOf(EvalEngine.getApfloat().hypergeometricU(fApfloat,
((IReal) arg2).apfloatValue(), ((IReal) arg3).apfloatValue()));
} catch (ArithmeticException | ApfloatRuntimeException e) {
// try as computation with complex numbers
}
}
if (arg2 instanceof INumber && arg3 instanceof INumber) {
return F.complexNum(EvalEngine.getApfloat().hypergeometricU(fApfloat,
((INumber) arg2).apcomplexValue(), ((INumber) arg3).apcomplexValue()));
}
return INum.super.hypergeometricU(arg2, arg3);
}

/** {@inheritDoc} */
@Override
public boolean isNumEqualInteger(IInteger ii) throws ArithmeticException {
Expand Down Expand Up @@ -210,37 +228,39 @@ public IExpr evaluate(EvalEngine engine) {

@Override
public IExpr erf() {
FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
try {
Apfloat erf = erf(fApfloat, h);
return F.num(erf);
} catch (Exception ce) {
//
}
return F.NIL;
}

private static Apfloat erf(Apfloat x, FixedPrecisionApfloatHelper h) {
Apint two = new Apint(2);
// 1/2
Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// 3/2
Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
Apfloat erf = h.hypergeometric1F1(oneHalf, threeHalf, h.multiply(x, x).negate()).multiply(two)
.multiply(x).divide(h.sqrt(h.pi()));
return erf;
}
return valueOf(EvalEngine.getApfloat().erf(fApfloat));
// FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
// try {
// Apfloat erf = erf(fApfloat, h);
// return F.num(erf);
// } catch (Exception ce) {
// //
// }
// return F.NIL;
}

// private static Apfloat erf(Apfloat x, FixedPrecisionApfloatHelper h) {
// Apint two = new Apint(2);
// // 1/2
// Aprational oneHalf = new Aprational(Apint.ONE, new Apint(2));
// // 3/2
// Aprational threeHalf = new Aprational(new Apint(3), new Apint(2));
// Apfloat erf = h.hypergeometric1F1(oneHalf, threeHalf, h.multiply(x, x).negate()).multiply(two)
// .multiply(x).divide(h.sqrt(h.pi()));
// return erf;
// }

@Override
public IExpr erfc() {
FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
try {
Apfloat c = erf(fApfloat, h);
return F.num(h.subtract(Apcomplex.ONE, c));
} catch (Exception ce) {
//
}
return F.NIL;
return valueOf(EvalEngine.getApfloat().erfc(fApfloat));
// FixedPrecisionApfloatHelper h = EvalEngine.getApfloat();
// try {
// Apfloat c = erf(fApfloat, h);
// return F.num(h.subtract(Apcomplex.ONE, c));
// } catch (Exception ce) {
// //
// }
// return F.NIL;
}

@Override
Expand Down
Loading

0 comments on commit c94a4d0

Please sign in to comment.