Skip to content

Commit

Permalink
Make numeric evaluation more compatible with MMA
Browse files Browse the repository at this point in the history 
Examples:
- ContinuedFraction(9.0/22) -> {0,2,2,3,1} instead of
{0.0,2.0,2.0,3.0,1.0}
- "MantissaExponent(3.4*10^30) -> {0.34,31}
- improve evaluation for NHoldFirst, NHolRest, NHoldAll
- improve `Listable` attribute evaluation: RealAbs({1.2, 1.5, 0}) ->
{1.2,1.5,0}
  • Loading branch information
@axkr
axkr committed Sep 3, 2023
1 parent 69b86d2 commit 05c1f16
Show file tree
Hide file tree
Showing 22 changed files with 402 additions and 233 deletions.
24 changes: 21 additions & 3 deletions ...droid_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/Arithmetic.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2569,6 +2569,9 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
@Override
public IExpr numericEval(final IAST ast, EvalEngine engine) {
IExpr arg1 = ast.arg1();
if (arg1.isList()) {
return ((IAST) arg1).mapThread(ast, 1);
}
if (ast.isAST1()) {
return numericEvalAST1(arg1, engine);
}
Expand All @@ -2588,7 +2591,7 @@ public IExpr numericEval(final IAST ast, EvalEngine engine) {
return Errors.printMessage(S.N, "precgt",
F.list(arg2, F.ZZ(Config.MAX_PRECISION_APFLOAT)), engine);
}
return numericEvalAST2(arg1, nDigitPrecision, engine);
return numericEvalAST2(arg1, arg2, nDigitPrecision, engine);
} finally {
engine.setNumericMode(oldNumericMode);
engine.setNumericPrecision(oldDigitPrecision);
Expand All @@ -2605,26 +2608,41 @@ private static IExpr numericEvalAST1(IExpr expr, EvalEngine engine) {
// avoid infinite recursions in symbolic mode
if (expr.isNumericFunction(true)) {
engine.setNumericMode(true, numericPrecision, oldSignificantFigures);
return engine.evalWithoutNumericReset(expr);
IExpr temp = engine.evalWithoutNumericReset(expr);
if (temp.isList()) {
return ((IAST) temp).mapThread(F.N(F.Slot1), 1);
}
return temp;
}
expr = engine.evaluate(expr);
if (expr.isInexactNumber()) {
return expr;
}
if (expr.isList()) {
return ((IAST) expr).mapThread(F.N(F.Slot1), 1);
}
engine.setNumericMode(true, numericPrecision, oldSignificantFigures);
if (expr.isAST()) {
ISymbol topSymbol = expr.topHead();
expr = engine.evalArgs((IAST) expr, topSymbol.getAttributes(), true).orElse(expr);
}
return engine.evalWithoutNumericReset(expr);
} finally {
engine.setNumericMode(oldNumericMode);
engine.setNumericPrecision(oldDigitPrecision);
}
}

private static IExpr numericEvalAST2(IExpr expr, long nDigitPrecision, EvalEngine engine) {
private static IExpr numericEvalAST2(IExpr expr, IExpr arg2, long nDigitPrecision,
EvalEngine engine) {
// first try symbolic evaluation
expr = engine.evaluate(expr);
if (expr.isInexactNumber()) {
return expr;
}
if (expr.isList()) {
return ((IAST) expr).mapThread(F.N(F.Slot1, arg2), 1);
}
final int maxSize =
(Config.MAX_OUTPUT_SIZE > Short.MAX_VALUE) ? Short.MAX_VALUE : Config.MAX_OUTPUT_SIZE;
int significantFigures = (nDigitPrecision > maxSize) ? maxSize : (int) nDigitPrecision;
Expand Down
4 changes: 2 additions & 2 deletions ...library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/BooleanFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4092,7 +4092,7 @@ private static final class SameQ extends AbstractCoreFunctionEvaluator
@Override
public IExpr evaluate(final IAST ast, EvalEngine engine) {
if (ast.size() > 2) {
IAST temp = engine.evalArgs(ast, ISymbol.NOATTRIBUTE).orElse(ast);
IAST temp = engine.evalArgs(ast, ISymbol.NOATTRIBUTE, false).orElse(ast);
if (temp.isAST2()) {
return temp.arg1().isSame(temp.arg2()) ? S.True : S.False;
}
Expand Down Expand Up @@ -4720,7 +4720,7 @@ private static final class UnsameQ extends AbstractCoreFunctionEvaluator
@Override
public IExpr evaluate(final IAST ast, EvalEngine engine) {
if (ast.size() > 2) {
IAST temp = engine.evalArgs(ast, ISymbol.NOATTRIBUTE).orElse(ast);
IAST temp = engine.evalArgs(ast, ISymbol.NOATTRIBUTE, false).orElse(ast);
if (temp.isAST2()) {
return temp.arg1().isSame(temp.arg2()) ? S.False : S.True;
}
Expand Down
4 changes: 3 additions & 1 deletion ...rary/matheclipse-core/src/main/java/org/matheclipse/core/builtin/ClusteringFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,9 @@ public double compute(double[] a, double[] b) throws MathIllegalArgumentExceptio

@Override
public IExpr distance(IExpr a, IExpr b, EvalEngine engine) {
return F.Divide(F.Total(F.Abs(F.Subtract(a, b))), F.Total(F.Abs(F.Plus(a, b))));
IExpr divide = F.Divide(F.Total(F.Abs(F.Subtract(a, b))), F.Total(F.Abs(F.Plus(a, b))));
divide = engine.evaluate(divide);
return divide;
}
}

Expand Down
86 changes: 40 additions & 46 deletions ...ibrary/matheclipse-core/src/main/java/org/matheclipse/core/builtin/EllipticIntegrals.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -978,15 +978,13 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
if (engine.isDoubleMode() && x.isNumber() && m.isNumber()) {
if (x.isReal() && m.isReal()) {
try {
return F
.complexNum(EllipticFunctionsJS.jacobiTheta(a, x.evalf(), m.evalf()));
return F.complexNum(EllipticFunctionsJS.jacobiTheta(a, x.evalf(), m.evalf()));
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
} else if (x.isInexactNumber() && m.isInexactNumber()) {
try {
return F.complexNum(
EllipticFunctionsJS.jacobiTheta(a, x.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.jacobiTheta(a, x.evalfc(), m.evalfc()));
} catch (ValidateException ve) {
return Errors.printMessage(ast.topHead(), ve, engine);
} catch (RuntimeException rex) {
Expand Down Expand Up @@ -1105,8 +1103,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
// } catch (ValidateException ve) {
// }
// if (Double.isNaN(zDouble) || Double.isNaN(mDouble)) {
return F
.complexNum(EllipticFunctionsJS.inverseJacobiCD(z.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.inverseJacobiCD(z.evalfc(), m.evalfc()));
// } else {
// return F.num(EllipticFunctionsJS.inverseJacobiCD(zDouble, mDouble));
// }
Expand Down Expand Up @@ -1165,8 +1162,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
// } catch (ValidateException ve) {
// }
// if (Double.isNaN(zDouble) || Double.isNaN(mDouble)) {
return F
.complexNum(EllipticFunctionsJS.inverseJacobiCN(z.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.inverseJacobiCN(z.evalfc(), m.evalfc()));
// } else {
// return F.num(EllipticFunctionsJS.inverseJacobiCN(zDouble, mDouble));
// }
Expand Down Expand Up @@ -1216,8 +1212,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
// } catch (ValidateException ve) {
// }
// if (Double.isNaN(zDouble) || Double.isNaN(mDouble)) {
return F
.complexNum(EllipticFunctionsJS.inverseJacobiDN(z.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.inverseJacobiDN(z.evalfc(), m.evalfc()));
// } else {
// return F.num(EllipticFunctionsJS.inverseJacobiDN(zDouble, mDouble));
// }
Expand Down Expand Up @@ -1281,8 +1276,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
// } catch (ValidateException ve) {
// }
// if (Double.isNaN(zDouble) || Double.isNaN(mDouble)) {
return F
.complexNum(EllipticFunctionsJS.inverseJacobiSC(z.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.inverseJacobiSC(z.evalfc(), m.evalfc()));
// } else {
// return F.num(EllipticFunctionsJS.inverseJacobiSC(zDouble, mDouble));
// }
Expand Down Expand Up @@ -1336,8 +1330,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
// } catch (ValidateException ve) {
// }
// if (Double.isNaN(zDouble) || Double.isNaN(mDouble)) {
return F
.complexNum(EllipticFunctionsJS.inverseJacobiSD(z.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.inverseJacobiSD(z.evalfc(), m.evalfc()));
// } else {
// return F.num(EllipticFunctionsJS.inverseJacobiSD(zDouble, mDouble));
// }
Expand Down Expand Up @@ -1397,8 +1390,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
// } catch (ValidateException ve) {
// }
// if (Double.isNaN(zDouble) || Double.isNaN(mDouble)) {
return F
.complexNum(EllipticFunctionsJS.inverseJacobiSN(z.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.inverseJacobiSN(z.evalfc(), m.evalfc()));
// } else {
// return F.num(EllipticFunctionsJS.inverseJacobiSN(zDouble, mDouble));
// }
Expand Down Expand Up @@ -1466,11 +1458,9 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
if (engine.isDoubleMode() && z.isNumber() && m.isNumber()) {
try {
if (z.isReal() && m.isReal()) {
return F
.complexNum(EllipticFunctionsJS.jacobiAmplitude(z.evalf(), m.evalf()));
return F.complexNum(EllipticFunctionsJS.jacobiAmplitude(z.evalf(), m.evalf()));
}
return F
.complexNum(EllipticFunctionsJS.jacobiAmplitude(z.evalfc(), m.evalfc()));
return F.complexNum(EllipticFunctionsJS.jacobiAmplitude(z.evalfc(), m.evalfc()));
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
return F.NIL;
Expand Down Expand Up @@ -2075,14 +2065,16 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
IAST list = (IAST) ast.arg1();
IExpr g2 = list.arg1();
IExpr g3 = list.arg2();
// numric mode isn't set here
if (g2.isInexactNumber() && g3.isInexactNumber()) {
try {
org.hipparchus.complex.Complex[] invariants =
EllipticFunctionsJS.weierstrassHalfPeriods(g2.evalfc(), g3.evalfc());
return Object2Expr.convertComplex(false, invariants);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
// numeric mode isn't set here
if (g2.isInexactNumber() || g3.isInexactNumber()) {
if (g2.isNumber() && g3.isNumber()) {
try {
org.hipparchus.complex.Complex[] invariants =
EllipticFunctionsJS.weierstrassHalfPeriods(g2.evalfc(), g3.evalfc());
return Object2Expr.convertComplex(false, invariants);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
}
}
}
Expand All @@ -2109,16 +2101,18 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
IAST list = (IAST) ast.arg1();
IExpr g2 = list.arg1();
IExpr g3 = list.arg2();
// numric mode isn't set here
if (g2.isInexactNumber() && g3.isInexactNumber()) {
try {
org.hipparchus.complex.Complex[] invariants =
EllipticFunctionsJS.weierstrassInvariants(g2.evalfc(), g3.evalfc());
return Object2Expr.convertComplex(false, invariants);
} catch (ValidateException ve) {
return Errors.printMessage(ast.topHead(), ve, engine);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
// numeric mode isn't set here
if (g2.isInexactNumber() || g3.isInexactNumber()) {
if (g2.isNumber() && g3.isNumber()) {
try {
org.hipparchus.complex.Complex[] invariants =
EllipticFunctionsJS.weierstrassInvariants(g2.evalfc(), g3.evalfc());
return Object2Expr.convertComplex(false, invariants);
} catch (ValidateException ve) {
return Errors.printMessage(ast.topHead(), ve, engine);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
}
}
}
Expand Down Expand Up @@ -2156,11 +2150,11 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
// 1 + (3/2) Cot(Sqrt(3/2)*u)^2
return F.Plus(F.C1, F.Times(F.C3D2, F.Sqr(F.Cot(F.Times(F.Sqrt(F.C3D2), u)))));
}
// numric mode isn't set here
if (u.isInexactNumber() && g2.isInexactNumber() && g3.isInexactNumber()) {
// numeric mode isn't set here
if (u.isInexactNumber() && g2.isNumber() && g3.isNumber()) {
try {
return F.complexNum(EllipticFunctionsJS.weierstrassP(u.evalfc(), g2.evalfc(),
g3.evalfc()));
return F
.complexNum(EllipticFunctionsJS.weierstrassP(u.evalfc(), g2.evalfc(), g3.evalfc()));
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
Expand Down Expand Up @@ -2201,11 +2195,11 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
return F.Times(F.CN3, F.Sqrt(F.C3D2), F.Cot(F.Times(F.Sqrt(F.C3D2), u)),
F.Sqr(F.Csc(F.Times(F.Sqrt(F.C3D2), u))));
}
// numric mode isn't set here
if (u.isInexactNumber() && g2.isInexactNumber() && g3.isInexactNumber()) {
// numeric mode isn't set here
if (u.isInexactNumber() && g2.isNumber() && g3.isNumber()) {
try {
return F.complexNum(EllipticFunctionsJS.weierstrassPPrime(u.evalfc(),
g2.evalfc(), g3.evalfc()));
return F.complexNum(
EllipticFunctionsJS.weierstrassPPrime(u.evalfc(), g2.evalfc(), g3.evalfc()));
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
Expand Down
63 changes: 31 additions & 32 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 @@ -1181,44 +1181,43 @@ public IExpr evaluate(IAST ast, EvalEngine engine) {
b = bVector;
}

// if (engine.isArbitraryMode() && a.isVector() > 0 && b.isVector() > 0) {
// try {
// HypergeometricHelper.hypergeometricPFQ();
// } catch (ValidateException ve) {
// return IOFunctions.printMessage(ast.topHead(), ve, engine);
// } catch (RuntimeException rex) {
// LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
// }
// }
// numeric mode isn't set here

if (engine.isDoubleMode() && a.isVector() > 0 && b.isVector() > 0) {
try {
double A[] = a.toDoubleVector();
double B[] = b.toDoubleVector();
double cDouble = Double.NaN;
try {
cDouble = c.evalf();
} catch (ValidateException ve) {
}
if (A == null || B == null || Double.isNaN(cDouble)) {
Complex AC[] = a.toComplexVector();
Complex BC[] = b.toComplexVector();
if (AC != null && BC != null) {
return F.complexNum(
HypergeometricJS.hypergeometricPFQ(AC, BC, c.evalfc(), Config.DOUBLE_TOLERANCE));
}
} else {
INum result = F.num(HypergeometricJS.hypergeometricPFQ(A, B, cDouble));
if (c.isInexactNumber() || a.isInexactVector() > 0 || b.isInexactVector() > 0) {
return numericHypergeometricPFQ(a, b, c, ast, engine);
}

return result;
}
return F.NIL;
}

private IExpr numericHypergeometricPFQ(IExpr a, IExpr b, IExpr c, IAST ast, EvalEngine engine) {
try {
double A[] = a.toDoubleVector();
double B[] = b.toDoubleVector();
double cDouble = Double.NaN;
try {
cDouble = c.evalf();
} catch (ValidateException ve) {
return Errors.printMessage(ast.topHead(), ve, engine);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}
if (A == null || B == null || Double.isNaN(cDouble)) {
Complex AC[] = a.toComplexVector();
Complex BC[] = b.toComplexVector();
if (AC != null && BC != null) {
return F.complexNum(
HypergeometricJS.hypergeometricPFQ(AC, BC, c.evalfc(), Config.DOUBLE_TOLERANCE));
}
} else {
INum result = F.num(HypergeometricJS.hypergeometricPFQ(A, B, cDouble));

return result;
}

} catch (ValidateException ve) {
return Errors.printMessage(ast.topHead(), ve, engine);
} catch (RuntimeException rex) {
LOGGER.log(engine.getLogLevel(), ast.topHead(), rex);
}

return F.NIL;
}

Expand Down
13 changes: 8 additions & 5 deletions ...id_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/LinearAlgebra.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2034,14 +2034,14 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
Times(C4, matrix.getEntry(0, 1), matrix.getEntry(1, 0)),
Times(CN2, matrix.getEntry(0, 0), matrix.getEntry(1, 1)),
Sqr(matrix.getEntry(1, 1))));
IExpr eigenValues = List(
IAST eigenValues = List(
Times(C1D2,
Plus(Negate(sqrtExpr), matrix.getEntry(0, 0), matrix.getEntry(1, 1))),
Times(C1D2, Plus(sqrtExpr, matrix.getEntry(0, 0), matrix.getEntry(1, 1))));
return sortValuesIfNumeric((IAST) eigenValues, arg2);
return sortValuesIfNumeric(eigenValues, arg2);
}
} else {
boolean hasNumericArgument = arg1.isEvalFlagOn(IAST.CONTAINS_NUMERIC_ARG);
boolean hasNumericArgument = arg1.isNumericArgument(true);// (IAST.CONTAINS_NUMERIC_ARG);
if (!hasNumericArgument) {
ISymbol x = F.Dummy("x");
IExpr m = engine.evaluate(F.CharacteristicPolynomial(arg1, x));
Expand All @@ -2065,7 +2065,8 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
// switch to numeric calculation
IExpr eigenValues = numericEval(ast, engine);
if (eigenValues.isList()) {
return sortValuesIfNumeric((IAST) eigenValues, arg2);
IAST sortFunction = sortValuesIfNumeric((IASTMutable) eigenValues, arg2);
return engine.evaluate(sortFunction);
}
return F.NIL;
}
Expand All @@ -2079,8 +2080,10 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
* elements are sorted in order of decreasing absolute value.
* @return
*/
private IExpr sortValuesIfNumeric(IAST eigenValuesList, final IExpr arg2) {
private IAST sortValuesIfNumeric(IAST eigenValuesList, final IExpr arg2) {
if (eigenValuesList.forAll(v -> v.isNumericFunction())) {
eigenValuesList = eigenValuesList.copy();
((IASTMutable) eigenValuesList).sortInplace();
if (arg2 != null && arg2.isPresent()) {
int n = arg2.toIntDefault();
if (n < 0) {
Expand Down
Loading

0 comments on commit 05c1f16

Please sign in to comment.