Skip to content

Commit

Permalink
WIP #942 Tan(Pi/2) // N returns wrong result
Browse files Browse the repository at this point in the history 
- new method Arithmetic#intPowerFractionNumeric()
  evaluate `Power(base,1/n)` in "double numeric mode" by "widen the
input domain of big integer numbers" to Apfloat values and calculating
the nth root.
  • Loading branch information
@axkr
axkr committed Mar 16, 2024
1 parent ac0a4a0 commit 90740c0
Show file tree
Hide file tree
Showing 6 changed files with 122 additions and 47 deletions.
70 changes: 61 additions & 9 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 @@ -99,6 +99,7 @@
import org.matheclipse.core.interfaces.IAST;
import org.matheclipse.core.interfaces.IASTAppendable;
import org.matheclipse.core.interfaces.IASTMutable;
import org.matheclipse.core.interfaces.IBigNumber;
import org.matheclipse.core.interfaces.IBuiltInSymbol;
import org.matheclipse.core.interfaces.IComplex;
import org.matheclipse.core.interfaces.IComplexNum;
Expand Down Expand Up @@ -2656,14 +2657,14 @@ private static IExpr numericEvalAST1(IExpr expr, EvalEngine engine) {
final int oldSignificantFigures = engine.getSignificantFigures();
try {
long numericPrecision = oldDigitPrecision; // Config.MACHINE_PRECISION;
if (expr.isNumericFunction(true)) {
engine.setNumericMode(true, numericPrecision, oldSignificantFigures);
IExpr temp = engine.evalWithoutNumericReset(expr);
if (temp.isListOrAssociation() || temp.isRuleAST()) {
return ((IAST) temp).mapThread(F.N(F.Slot1), 1);
}
return temp;
}
// if (expr.isNumericFunction(true)) {
// engine.setNumericMode(true, numericPrecision, oldSignificantFigures);
// IExpr temp = engine.evalWithoutNumericReset(expr);
// if (temp.isListOrAssociation() || temp.isRuleAST()) {
// return ((IAST) temp).mapThread(F.N(F.Slot1), 1);
// }
// return temp;
// }
expr = engine.evaluate(expr);
if (expr.isInexactNumber()) {
return expr;
Expand Down Expand Up @@ -3561,7 +3562,7 @@ public static IExpr sqrtDenest(IRational arg1, IExpr arg2) {
return F.NIL;
}

public IExpr binaryOperator(IAST ast, final IExpr base, final IExpr exponent,
public static IExpr binaryOperator(IAST ast, final IExpr base, final IExpr exponent,
EvalEngine engine) {
try {
if (base.isInexactNumber() && exponent.isInexactNumber()) {
Expand Down Expand Up @@ -7034,6 +7035,57 @@ protected ISymbol getArithmeticSymbol() {

}

/**
* Eval in double numeric mode by "widen the input domain" to Apfloat values.
*
* @param powerAST2 "binary {@link S#Power} function"
* @return
*/
public static IExpr intPowerFractionNumeric(IAST powerAST2) {
IExpr base = powerAST2.base();
IExpr exponent = powerAST2.exponent();
if ((base instanceof IBigNumber) && exponent.isFraction()) {
IFraction exp = (IFraction) exponent;
int denom = exp.denominator().toIntDefault();
if (denom > 0L) {
if (base.isRational()) {
IRational iBase = (IRational) base;
double fNum = base.evalf();
if (!Double.isFinite(fNum) || fNum <= Double.MIN_VALUE || fNum >= Double.MAX_VALUE) {
if (iBase.isPositive()) {
// special case root of "rational base"
ApfloatNum apfloat = iBase.apfloatNumValue();
if (exp.numerator().isOne()) {
return F.num(apfloat.rootN(denom).doubleValue());
} else if (exp.numerator().isMinusOne()) {
return F.num(apfloat.rootN(denom).inverse().doubleValue());
}
} else if (iBase.isNegative()) {
ApcomplexNum apcomplex = iBase.apcomplexNumValue();
if (exp.numerator().isOne()) {
return F.complexNum(apcomplex.rootN(denom).evalfc());
} else if (exp.numerator().isMinusOne()) {
return F.complexNum(apcomplex.rootN(denom).inverse().evalfc());
}
}
}
} else if (base.isComplex()) {
IComplex iBase = (IComplex) base;
org.hipparchus.complex.Complex fComplex = base.evalfc();
if (!fComplex.isFinite()) {
ApcomplexNum apcomplex = iBase.apcomplexNumValue();
if (exp.numerator().isOne()) {
return F.complexNum(apcomplex.rootN(denom).evalfc());
} else if (exp.numerator().isMinusOne()) {
return F.complexNum(apcomplex.rootN(denom).inverse().evalfc());
}
}
}
}
}
return F.NIL;
}

/**
* Try simpplifying <code>(power0Arg1 ^ power0Arg2) * (power1Arg1 ^ power1Arg2)</code>
*
Expand Down
37 changes: 30 additions & 7 deletions ..._android_library/matheclipse-core/src/main/java/org/matheclipse/core/eval/EvalEngine.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -849,6 +849,12 @@ public IASTMutable evalArgs(final IAST ast, final int attributes, boolean numeri
final boolean isNumericFunction;
if ((ISymbol.NUMERICFUNCTION & attributes) == ISymbol.NUMERICFUNCTION) {
isNumericFunction = true;
if (isDoubleMode() && ast.isPower()) {
IExpr temp = Arithmetic.intPowerFractionNumeric(ast);
if (temp.isPresent()) {
return F.unaryAST1(S.Power, temp);
}
}
} else {
isNumericFunction = numericFunction;
}
Expand Down Expand Up @@ -1633,7 +1639,7 @@ public final double evalDouble(IExpr expr, Function<IExpr, IExpr> function, doub
if (function != null) {
expr = expr.accept(new VisitorReplaceEvalf(function)).orElse(expr);
}
IExpr result = evalN(expr);
IExpr result = evalNumericFunction(expr);
if (result.isReal()) {
return ((IReal) result).doubleValue();
}
Expand All @@ -1648,12 +1654,29 @@ public final double evalDouble(IExpr expr, Function<IExpr, IExpr> function, doub
if (F.isZero(cc.getImaginaryPart())) {
return cc.getRealPart();
}
}
if (result.isQuantity()) {
return result.evalReal().doubleValue();
}
if (result.isAST(S.Labeled, 3, 4)) {
return result.first().evalReal().doubleValue();
} else {
result = evalN(expr);
if (result.isReal()) {
return ((IReal) result).doubleValue();
}
if (result.isInfinity()) {
return Double.POSITIVE_INFINITY;
}
if (result.isNegativeInfinity()) {
return Double.NEGATIVE_INFINITY;
}
if (result.isComplexNumeric()) {
IComplexNum cc = (IComplexNum) result;
if (F.isZero(cc.getImaginaryPart())) {
return cc.getRealPart();
}
}
if (result.isQuantity()) {
return result.evalReal().doubleValue();
}
if (result.isAST(S.Labeled, 3, 4)) {
return result.first().evalReal().doubleValue();
}
}
} finally {
fQuietMode = quietMode;
Expand Down
14 changes: 7 additions & 7 deletions ...library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6580,15 +6580,15 @@ public void testEllipticF() {
@Test
public void testEllipticK() {

check("N(EllipticK(8/10), 50)", //
"2.2572053268208536550832560045233873972354192817399");

check("EllipticK(0.999999999999999990000000000000000)", //
"20.9582676515692789828830607366566");
// reducing the accuracy gives ComplexInfinity in MMA:
check("EllipticK(0.99999999999999999)", //
"20.958267651569278");

check("N(EllipticK(8/10), 50)", //
"2.2572053268208536550832560045233873972354192817399");

checkNumeric("EllipticK(2.5+I)", //
"1.1551450606569331+I*0.9528453714670536");

Expand Down Expand Up @@ -15997,8 +15997,8 @@ public void testDeterminePrecision() {
@Test
public void testN() {
// issue #942
// check("Tan(Pi/2) // N", //
// "ComplexInfinity");
check("Tan(Pi/2) // N", //
"ComplexInfinity");
// issue #937
check("(x==-157079632679/100000000000) // N", //
"x==-1.5708");
Expand Down Expand Up @@ -24208,11 +24208,11 @@ public void testSurd() {
checkNumeric("Surd(-3,3)", //
"-3^(1/3)");
checkNumeric("N((-3)^(1/3))", //
"0.7211247851537043+I*1.2490247664834064");
"0.7211247851537042+I*1.2490247664834064");
checkNumeric("Surd(-3,3)-(-3)^(1/3)", //
"-(-3)^(1/3)-3^(1/3)");
checkNumeric("Surd(-3.,3)-(-3)^(1/3)", //
"-2.1633743554611127+I*(-1.2490247664834064)");
"-2.1633743554611122+I*(-1.2490247664834064)");
checkNumeric("Surd(-3,3)", //
"-3^(1/3)");
checkNumeric("Surd(-3.,3)", //
Expand Down
6 changes: 3 additions & 3 deletions ...roid_library/matheclipse-core/src/test/java/org/matheclipse/core/system/MainTestCase.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -918,10 +918,10 @@ public void testSystem040() {
"2.3715183290419594E-16+I*3.872983346207417");
checkNumeric(".5^.5", //
"0.7071067811865476");
// checkNumeric("N((-15)^(1/2))", //
// "I*3.872983346207417");
checkNumeric("N((-15)^(1/2))", //
"2.3715183290419594E-16+I*3.872983346207417");
"I*3.872983346207417");
// checkNumeric("N((-15)^(1/2))", //
// "2.3715183290419594E-16+I*3.872983346207417");
checkNumeric("N(Sin(1/2))", //
"0.479425538604203");
checkNumeric("N(1/6*(I*44^(1/2)+2))", //
Expand Down
6 changes: 3 additions & 3 deletions ...y/matheclipse-core/src/test/java/org/matheclipse/core/system/PolynomialFunctionsTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,10 +74,10 @@ public void testGegenbauerC() {
"{ComplexInfinity,ComplexInfinity,1.17445+I*2.30854,ComplexInfinity}");
check("GegenbauerC(2, 0.5)", //
"-0.5");
// checkNumeric("GegenbauerC(5,1/8,7) //N", //
// "16839.531372070312");
checkNumeric("GegenbauerC(5,1/8,7) //N", //
"16839.53137207032");
"16839.531372070312");
// checkNumeric("GegenbauerC(5,1/8,7) //N", //
// "16839.53137207032");
checkNumeric("Table(GegenbauerC(10, x), {x, 1, 5})", //
"{1/5,262087/5,22619537/5,457470751/5,4517251249/5}");
check("GegenbauerC({1/3, 1/2}, 1/6, 0)", //
Expand Down
36 changes: 18 additions & 18 deletions ...android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SolveTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -393,15 +393,15 @@ public void testSolve() {
"{{x->-a},{x->b}}");

// github #261 - JUnit test for Apfloat switching to complex Power calculation
// checkNumeric("Solve(0.00004244131815783 == x^5 , x)", //
// "{{x->-0.10802279680851234+I*0.07848315587546606},"//
// + "{x->-0.10802279680851212+I*(-0.07848315587546606)},"//
// + "{x->0.04126103682102799+I*(-0.1269884137508598)},"//
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
check("Solve(0.00004244131815783 == x^5 , x)", //
"{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},"//
+ "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");
checkNumeric("Solve(0.00004244131815783 == x^5 , x)", //
"{{x->-0.10802279680851234+I*0.07848315587546606},"//
+ "{x->-0.10802279680851212+I*(-0.07848315587546606)},"//
+ "{x->0.04126103682102799+I*(-0.1269884137508598)},"//
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
// check("Solve(0.00004244131815783 == x^5 , x)", //
// "{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},"//
// + "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");

// github #247
check("Solve((k+3)/(4)==(k)/2,{k})", //
Expand Down Expand Up @@ -801,15 +801,15 @@ public void testNSolve() {
"{{x->3.0,y->5.0},{x->6.34781+I*(-1.02885),y->1.75806+I*(-2.7734)},{x->6.34781+I*1.02885,y->1.75806+I*2.7734},{x->4.30438,y->1.48389}}");

// github #261 - JUnit test for Apfloat switching to complex Power calculation
// checkNumeric("NSolve(0.00004244131815783 == x^5 , x)", //
// "{{x->-0.10802279680851234+I*0.07848315587546606}," //
// + "{x->-0.10802279680851212+I*(-0.07848315587546606)}," //
// + "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
check("NSolve(0.00004244131815783 == x^5 , x)", //
"{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},{x->0.04126103682102799+I*(-0.1269884137508598)},"
//
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");
checkNumeric("NSolve(0.00004244131815783 == x^5 , x)", //
"{{x->-0.10802279680851234+I*0.07848315587546606}," //
+ "{x->-0.10802279680851212+I*(-0.07848315587546606)}," //
+ "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
// check("NSolve(0.00004244131815783 == x^5 , x)", //
// "{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},{x->0.04126103682102799+I*(-0.1269884137508598)},"
// //
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");


// github #247
Expand Down

0 comments on commit 90740c0

Please sign in to comment.