Skip to content

Commit ddba370

Browse files
axkraxkr
committed
WIP #929: use Apfloat's high precision polylog and isZero methods
- mtommila/apfloat#7 - mtommila/apfloat#47
1 parent e8eefde commit ddba370

File tree

6 files changed

+93
-70
lines changed

6 files changed

+93
-70
lines changed

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SpecialFunctions.java

+11-25
Original file line numberDiff line numberDiff line change
@@ -1949,29 +1949,7 @@ public IExpr numericFunction(IAST ast, final EvalEngine engine) {
19491949
return temp;
19501950
}
19511951
// issue #929
1952-
// return v.polyLog(z);
1953-
1954-
if (engine.isDoubleMode()) {
1955-
try {
1956-
double nDouble = Double.NaN;
1957-
double xDouble = Double.NaN;
1958-
try {
1959-
nDouble = v.evalf();
1960-
xDouble = z.evalf();
1961-
} catch (ValidateException ve) {
1962-
}
1963-
1964-
if (Double.isNaN(nDouble) || Double.isNaN(xDouble)) {
1965-
Complex nComplex = v.evalfc();
1966-
Complex xComplex = z.evalfc();
1967-
return F.complexNum(ZetaJS.polyLog(nComplex, xComplex));
1968-
} else {
1969-
return F.complexNum(ZetaJS.polyLog(nDouble, xDouble));
1970-
}
1971-
} catch (RuntimeException rex) {
1972-
Errors.printMessage(S.PolyLog, rex, engine);
1973-
}
1974-
}
1952+
return v.polyLog(z);
19751953
}
19761954
return F.NIL;
19771955
}
@@ -2001,6 +1979,14 @@ private IExpr polyLogSymbolic(IExpr n, IExpr z) {
20011979
}
20021980

20031981
if (n.isReal()) {
1982+
if (z.isInfinity() || z.isNegativeInfinity()) {
1983+
if (n.isZero() || n.isMathematicalIntegerNegative()) {
1984+
return S.Indeterminate;
1985+
}
1986+
if (n.isMathematicalIntegerNonNegative()) {
1987+
return F.CNInfinity;
1988+
}
1989+
}
20041990
if (n.isZero()) {
20051991
// arg2/(1 - arg2)
20061992
return Times(z, Power(Plus(C1, Negate(z)), -1));
@@ -2150,8 +2136,8 @@ private static IExpr functionExpandLogArg(final IExpr o) {
21502136
IFraction powExponent = (IFraction) power.exponent();
21512137
if (n.equals(powExponent.denominator())) {
21522138
EvalEngine engine = EvalEngine.get();
2153-
IExpr a = engine
2154-
.evalNumericFunction(F.Times(n, F.ProductLog(F.Times(arg1, F.Power(b, powExponent), F.Log(b))),
2139+
IExpr a = engine.evalNumericFunction(
2140+
F.Times(n, F.ProductLog(F.Times(arg1, F.Power(b, powExponent), F.Log(b))),
21552141
F.Power(F.Log(b), F.CN1)));
21562142
IExpr v = engine.evaluate(F.Rationalize(a));
21572143
if (v.isInteger() && ((IInteger) v).isGE(F.C1)) {

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ApcomplexNum.java

+10-9
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
import org.apfloat.ApfloatMath;
88
import org.apfloat.ApfloatRuntimeException;
99
import org.apfloat.FixedPrecisionApfloatHelper;
10+
import org.apfloat.LossOfPrecisionException;
1011
import org.apfloat.OverflowException;
1112
import org.hipparchus.complex.Complex;
1213
import org.matheclipse.core.basic.Config;
@@ -1001,14 +1002,12 @@ public INumber inverse() {
10011002

10021003
@Override
10031004
public boolean isMinusOne() {
1004-
return fApcomplex.real().compareTo(ApfloatNum.MINUS_ONE) == 0 //
1005-
&& fApcomplex.imag().signum() == 0;
1005+
return fApcomplex.imag().isZero() && fApcomplex.real().compareTo(ApfloatNum.MINUS_ONE) == 0;
10061006
}
10071007

10081008
@Override
10091009
public boolean isOne() {
1010-
return fApcomplex.real().compareTo(Apfloat.ONE) == 0 //
1011-
&& fApcomplex.imag().signum() == 0;
1010+
return fApcomplex.imag().isZero() && fApcomplex.real().compareTo(Apfloat.ONE) == 0;
10121011
}
10131012

10141013
@Override
@@ -1021,9 +1020,7 @@ public boolean isSame(IExpr expression, double epsilon) {
10211020

10221021
@Override
10231022
public boolean isZero() {
1024-
// return fApcomplex.isZero();
1025-
return fApcomplex.real().signum() == 0 //
1026-
&& fApcomplex.imag().signum() == 0;
1023+
return fApcomplex.isZero();
10271024
}
10281025

10291026
@Override
@@ -1165,8 +1162,12 @@ public IExpr polyGamma(long n) {
11651162
@Override
11661163
public IExpr polyLog(IExpr arg2) {
11671164
if (arg2 instanceof INumber) {
1168-
return valueOf(
1169-
EvalEngine.getApfloat().polylog(fApcomplex, ((INumber) arg2).apcomplexValue()));
1165+
try {
1166+
return valueOf(
1167+
EvalEngine.getApfloat().polylog(fApcomplex, ((INumber) arg2).apcomplexValue()));
1168+
} catch (LossOfPrecisionException lope) {
1169+
// Complete loss of precision
1170+
}
11701171
}
11711172
return IComplexNum.super.polyLog(arg2);
11721173
}

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ApfloatNum.java

+15-14
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
import org.apfloat.ApfloatRuntimeException;
99
import org.apfloat.Apint;
1010
import org.apfloat.FixedPrecisionApfloatHelper;
11+
import org.apfloat.LossOfPrecisionException;
1112
import org.apfloat.OverflowException;
1213
import org.matheclipse.core.basic.Config;
1314
import org.matheclipse.core.eval.Errors;
@@ -1085,10 +1086,6 @@ public boolean isOne() {
10851086
return fApfloat.compareTo(Apfloat.ONE) == 0;
10861087
}
10871088

1088-
// public ApfloatNum apfloatNumValue() {
1089-
// return this;
1090-
// }
1091-
10921089
/** {@inheritDoc} */
10931090
@Override
10941091
public boolean isPi() {
@@ -1118,7 +1115,7 @@ public boolean isSame(IExpr expression, double epsilon) {
11181115
/** {@inheritDoc} */
11191116
@Override
11201117
public boolean isZero() {
1121-
return fApfloat.signum() == 0;
1118+
return fApfloat.isZero();
11221119
}
11231120

11241121
@Override
@@ -1287,18 +1284,22 @@ public IExpr polyGamma(long n) {
12871284

12881285
@Override
12891286
public IExpr polyLog(IExpr arg2) {
1290-
if (arg2 instanceof IReal) {
1287+
if (arg2 instanceof INumber) {
1288+
if (arg2 instanceof IReal) {
1289+
try {
1290+
return valueOf(EvalEngine.getApfloat().polylog(fApfloat, ((IReal) arg2).apfloatValue()));
1291+
} catch (ArithmeticException | ApfloatRuntimeException e) {
1292+
// java.lang.ArithmeticException: Result would be complex
1293+
}
1294+
}
1295+
12911296
try {
1292-
return valueOf(
1293-
EvalEngine.getApfloat().polylog(fApfloat, ((IReal) arg2).apfloatValue()));
1294-
} catch (ArithmeticException | ApfloatRuntimeException e) {
1295-
// try as computation with complex numbers
1297+
return F.complexNum(
1298+
EvalEngine.getApfloat().polylog(fApfloat, ((INumber) arg2).apcomplexValue()));
1299+
} catch (LossOfPrecisionException lope) {
1300+
// Complete loss of precision
12961301
}
12971302
}
1298-
if (arg2 instanceof INumber) {
1299-
return F.complexNum(
1300-
EvalEngine.getApfloat().polylog(fApfloat, ((INumber) arg2).apcomplexValue()));
1301-
}
13021303
return INum.super.polyLog(arg2);
13031304
}
13041305

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/ComplexNum.java

+8-3
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
import org.apfloat.Apcomplex;
55
import org.apfloat.Apfloat;
66
import org.apfloat.ApfloatRuntimeException;
7+
import org.apfloat.LossOfPrecisionException;
78
import org.apfloat.OverflowException;
89
import org.hipparchus.complex.Complex;
910
import org.hipparchus.exception.NullArgumentException;
@@ -1326,9 +1327,13 @@ public IExpr polyGamma(long n) {
13261327
@Override
13271328
public IExpr polyLog(IExpr arg2) {
13281329
if (arg2 instanceof INumber) {
1329-
Apcomplex polylog = EvalEngine.getApfloatDouble().polylog(apcomplexValue(),
1330-
((INumber) arg2).apcomplexValue());
1331-
return F.complexNum(polylog.real().doubleValue(), polylog.imag().doubleValue());
1330+
try {
1331+
Apcomplex polylog = EvalEngine.getApfloatDouble().polylog(apcomplexValue(),
1332+
((INumber) arg2).apcomplexValue());
1333+
return F.complexNum(polylog.real().doubleValue(), polylog.imag().doubleValue());
1334+
} catch (LossOfPrecisionException lope) {
1335+
// Complete loss of precision
1336+
}
13321337
}
13331338
return IComplexNum.super.polyLog(arg2);
13341339
}

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/Num.java

+15-6
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
import org.apfloat.Apcomplex;
55
import org.apfloat.Apfloat;
66
import org.apfloat.ApfloatRuntimeException;
7+
import org.apfloat.LossOfPrecisionException;
78
import org.apfloat.OverflowException;
89
import org.hipparchus.complex.Complex;
910
import org.hipparchus.util.MathUtils;
@@ -1345,13 +1346,21 @@ public IExpr polyGamma(long n) {
13451346
public IExpr polyLog(IExpr arg2) {
13461347
if (arg2 instanceof INumber) {
13471348
if (arg2 instanceof IReal) {
1348-
Apfloat polylog =
1349-
EvalEngine.getApfloatDouble().polylog(apfloatValue(), ((IReal) arg2).apfloatValue());
1350-
return F.num(polylog.doubleValue());
1349+
try {
1350+
Apfloat polylog =
1351+
EvalEngine.getApfloatDouble().polylog(apfloatValue(), ((IReal) arg2).apfloatValue());
1352+
return F.num(polylog.doubleValue());
1353+
} catch (ArithmeticException | ApfloatRuntimeException e) {
1354+
// java.lang.ArithmeticException: Result would be complex
1355+
}
1356+
}
1357+
try {
1358+
Apcomplex polylog = EvalEngine.getApfloatDouble().polylog(apfloatValue(),
1359+
((INumber) arg2).apcomplexValue());
1360+
return F.complexNum(polylog.real().doubleValue(), polylog.imag().doubleValue());
1361+
} catch (LossOfPrecisionException lope) {
1362+
// Complete loss of precision
13511363
}
1352-
Apcomplex polylog = EvalEngine.getApfloatDouble().polylog(apfloatValue(),
1353-
((INumber) arg2).apcomplexValue());
1354-
return F.complexNum(polylog.real().doubleValue(), polylog.imag().doubleValue());
13551364
}
13561365
return INum.super.polyLog(arg2);
13571366
}

‎symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java

+34-13
Original file line numberDiff line numberDiff line change
@@ -18649,12 +18649,37 @@ public void testPolyGamma() {
1864918649

1865018650
@Test
1865118651
public void testPolyLog() {
18652+
check("PolyLog(-42,Infinity)", //
18653+
"Indeterminate");
18654+
check("PolyLog(0,Infinity)", //
18655+
"Indeterminate");
18656+
check("PolyLog(2,Infinity)", //
18657+
"-Infinity");
18658+
check("PolyLog(42,-Infinity)", //
18659+
"-Infinity");
18660+
18661+
checkNumeric("PolyLog(1, 0.333333)", //
18662+
"0.4054646081082894");
18663+
18664+
18665+
checkNumeric("PolyLog(2, 0.9)", //
18666+
"1.2997147230049588");
18667+
checkNumeric("PolyLog(0, 5.0)", //
18668+
"-1.25");
18669+
18670+
checkNumeric("N(PolyLog(1, 1/3), 50)", //
18671+
"0.40546510810816438197801311546434913657199042346249");
18672+
checkNumeric("PolyLog(2, 0.300000000000000000)", //
18673+
"0.326129510075476069");
18674+
checkNumeric("PolyLog(0.2 + I, 0.5 - I)", //
18675+
"-0.08985258966284129+I*(-0.5958648241210646)");
18676+
1865218677
// issue #929
1865318678
// message Infinite or NaN number in z1 calculation.
18654-
check("(PolyLog(2,E^(I*1/270*Pi^2)))/Pi // N", //
18655-
"0.31831*PolyLog(2.0,0.999332+I*0.036546)");
18656-
check("PolyLog(2,0.9993319736282411 + 0.03654595031305655*I)", //
18657-
"PolyLog(2,0.9993319736282411+I*0.03654595031305655)");
18679+
checkNumeric("(PolyLog(2,E^(I*1/270*Pi^2)))/Pi // N", //
18680+
"0.5054280619497805+I*0.0501372676125785");
18681+
checkNumeric("PolyLog(2,0.9993319736282411 + 0.03654595031305655*I)", //
18682+
"1.5878490863395573+I*0.1575108716027421");
1865818683

1865918684
check("PolyLog(2,z) + PolyLog(2,1-z)", //
1866018685
"Pi^2/6-Log(1-z)*Log(z)");
@@ -18667,6 +18692,11 @@ public void testPolyLog() {
1866718692
"Pi^2/6");
1866818693
check("PolyLog(2,E^(41*I*Pi))", //
1866918694
"-Pi^2/12");
18695+
18696+
18697+
18698+
18699+
1867018700
// TODO https://github.com/mtommila/apfloat/issues/34
1867118701
// check("PolyLog(2147483647,-3.1415)", //
1867218702
// " ");
@@ -18683,15 +18713,6 @@ public void testPolyLog() {
1868318713
check("PolyLog(-2147483648,2.718281828459045)", //
1868418714
"PolyLog(-2.14748*10^9,2.71828)");
1868518715

18686-
check("PolyLog(0.2 + I, 0.5 - I)", //
18687-
"-0.0898526+I*(-0.595865)");
18688-
check("PolyLog(2, 0.9)", //
18689-
"1.29971");
18690-
check("PolyLog(0, 5.0)", //
18691-
"-1.25");
18692-
check("PolyLog(1, 0.333333)", //
18693-
"0.405465");
18694-
1869518716
check("PolyLog(2,0)", //
1869618717
"0");
1869718718
check("PolyLog(2,-1)", //

0 commit comments

Comments
 (0)