Skip to content

Commit 90740c0

Browse files
axkraxkr
committed
WIP #942 Tan(Pi/2) // N returns wrong result
- 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.
1 parent ac0a4a0 commit 90740c0

File tree

6 files changed

+122
-47
lines changed

6 files changed

+122
-47
lines changed

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

+61-9
Original file line numberDiff line numberDiff line change
@@ -99,6 +99,7 @@
9999
import org.matheclipse.core.interfaces.IAST;
100100
import org.matheclipse.core.interfaces.IASTAppendable;
101101
import org.matheclipse.core.interfaces.IASTMutable;
102+
import org.matheclipse.core.interfaces.IBigNumber;
102103
import org.matheclipse.core.interfaces.IBuiltInSymbol;
103104
import org.matheclipse.core.interfaces.IComplex;
104105
import org.matheclipse.core.interfaces.IComplexNum;
@@ -2656,14 +2657,14 @@ private static IExpr numericEvalAST1(IExpr expr, EvalEngine engine) {
26562657
final int oldSignificantFigures = engine.getSignificantFigures();
26572658
try {
26582659
long numericPrecision = oldDigitPrecision; // Config.MACHINE_PRECISION;
2659-
if (expr.isNumericFunction(true)) {
2660-
engine.setNumericMode(true, numericPrecision, oldSignificantFigures);
2661-
IExpr temp = engine.evalWithoutNumericReset(expr);
2662-
if (temp.isListOrAssociation() || temp.isRuleAST()) {
2663-
return ((IAST) temp).mapThread(F.N(F.Slot1), 1);
2664-
}
2665-
return temp;
2666-
}
2660+
// if (expr.isNumericFunction(true)) {
2661+
// engine.setNumericMode(true, numericPrecision, oldSignificantFigures);
2662+
// IExpr temp = engine.evalWithoutNumericReset(expr);
2663+
// if (temp.isListOrAssociation() || temp.isRuleAST()) {
2664+
// return ((IAST) temp).mapThread(F.N(F.Slot1), 1);
2665+
// }
2666+
// return temp;
2667+
// }
26672668
expr = engine.evaluate(expr);
26682669
if (expr.isInexactNumber()) {
26692670
return expr;
@@ -3561,7 +3562,7 @@ public static IExpr sqrtDenest(IRational arg1, IExpr arg2) {
35613562
return F.NIL;
35623563
}
35633564

3564-
public IExpr binaryOperator(IAST ast, final IExpr base, final IExpr exponent,
3565+
public static IExpr binaryOperator(IAST ast, final IExpr base, final IExpr exponent,
35653566
EvalEngine engine) {
35663567
try {
35673568
if (base.isInexactNumber() && exponent.isInexactNumber()) {
@@ -7034,6 +7035,57 @@ protected ISymbol getArithmeticSymbol() {
70347035

70357036
}
70367037

7038+
/**
7039+
* Eval in double numeric mode by "widen the input domain" to Apfloat values.
7040+
*
7041+
* @param powerAST2 "binary {@link S#Power} function"
7042+
* @return
7043+
*/
7044+
public static IExpr intPowerFractionNumeric(IAST powerAST2) {
7045+
IExpr base = powerAST2.base();
7046+
IExpr exponent = powerAST2.exponent();
7047+
if ((base instanceof IBigNumber) && exponent.isFraction()) {
7048+
IFraction exp = (IFraction) exponent;
7049+
int denom = exp.denominator().toIntDefault();
7050+
if (denom > 0L) {
7051+
if (base.isRational()) {
7052+
IRational iBase = (IRational) base;
7053+
double fNum = base.evalf();
7054+
if (!Double.isFinite(fNum) || fNum <= Double.MIN_VALUE || fNum >= Double.MAX_VALUE) {
7055+
if (iBase.isPositive()) {
7056+
// special case root of "rational base"
7057+
ApfloatNum apfloat = iBase.apfloatNumValue();
7058+
if (exp.numerator().isOne()) {
7059+
return F.num(apfloat.rootN(denom).doubleValue());
7060+
} else if (exp.numerator().isMinusOne()) {
7061+
return F.num(apfloat.rootN(denom).inverse().doubleValue());
7062+
}
7063+
} else if (iBase.isNegative()) {
7064+
ApcomplexNum apcomplex = iBase.apcomplexNumValue();
7065+
if (exp.numerator().isOne()) {
7066+
return F.complexNum(apcomplex.rootN(denom).evalfc());
7067+
} else if (exp.numerator().isMinusOne()) {
7068+
return F.complexNum(apcomplex.rootN(denom).inverse().evalfc());
7069+
}
7070+
}
7071+
}
7072+
} else if (base.isComplex()) {
7073+
IComplex iBase = (IComplex) base;
7074+
org.hipparchus.complex.Complex fComplex = base.evalfc();
7075+
if (!fComplex.isFinite()) {
7076+
ApcomplexNum apcomplex = iBase.apcomplexNumValue();
7077+
if (exp.numerator().isOne()) {
7078+
return F.complexNum(apcomplex.rootN(denom).evalfc());
7079+
} else if (exp.numerator().isMinusOne()) {
7080+
return F.complexNum(apcomplex.rootN(denom).inverse().evalfc());
7081+
}
7082+
}
7083+
}
7084+
}
7085+
}
7086+
return F.NIL;
7087+
}
7088+
70377089
/**
70387090
* Try simpplifying <code>(power0Arg1 ^ power0Arg2) * (power1Arg1 ^ power1Arg2)</code>
70397091
*

‎symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/eval/EvalEngine.java

+30-7
Original file line numberDiff line numberDiff line change
@@ -849,6 +849,12 @@ public IASTMutable evalArgs(final IAST ast, final int attributes, boolean numeri
849849
final boolean isNumericFunction;
850850
if ((ISymbol.NUMERICFUNCTION & attributes) == ISymbol.NUMERICFUNCTION) {
851851
isNumericFunction = true;
852+
if (isDoubleMode() && ast.isPower()) {
853+
IExpr temp = Arithmetic.intPowerFractionNumeric(ast);
854+
if (temp.isPresent()) {
855+
return F.unaryAST1(S.Power, temp);
856+
}
857+
}
852858
} else {
853859
isNumericFunction = numericFunction;
854860
}
@@ -1633,7 +1639,7 @@ public final double evalDouble(IExpr expr, Function<IExpr, IExpr> function, doub
16331639
if (function != null) {
16341640
expr = expr.accept(new VisitorReplaceEvalf(function)).orElse(expr);
16351641
}
1636-
IExpr result = evalN(expr);
1642+
IExpr result = evalNumericFunction(expr);
16371643
if (result.isReal()) {
16381644
return ((IReal) result).doubleValue();
16391645
}
@@ -1648,12 +1654,29 @@ public final double evalDouble(IExpr expr, Function<IExpr, IExpr> function, doub
16481654
if (F.isZero(cc.getImaginaryPart())) {
16491655
return cc.getRealPart();
16501656
}
1651-
}
1652-
if (result.isQuantity()) {
1653-
return result.evalReal().doubleValue();
1654-
}
1655-
if (result.isAST(S.Labeled, 3, 4)) {
1656-
return result.first().evalReal().doubleValue();
1657+
} else {
1658+
result = evalN(expr);
1659+
if (result.isReal()) {
1660+
return ((IReal) result).doubleValue();
1661+
}
1662+
if (result.isInfinity()) {
1663+
return Double.POSITIVE_INFINITY;
1664+
}
1665+
if (result.isNegativeInfinity()) {
1666+
return Double.NEGATIVE_INFINITY;
1667+
}
1668+
if (result.isComplexNumeric()) {
1669+
IComplexNum cc = (IComplexNum) result;
1670+
if (F.isZero(cc.getImaginaryPart())) {
1671+
return cc.getRealPart();
1672+
}
1673+
}
1674+
if (result.isQuantity()) {
1675+
return result.evalReal().doubleValue();
1676+
}
1677+
if (result.isAST(S.Labeled, 3, 4)) {
1678+
return result.first().evalReal().doubleValue();
1679+
}
16571680
}
16581681
} finally {
16591682
fQuietMode = quietMode;

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

+7-7
Original file line numberDiff line numberDiff line change
@@ -6580,15 +6580,15 @@ public void testEllipticF() {
65806580
@Test
65816581
public void testEllipticK() {
65826582

6583+
check("N(EllipticK(8/10), 50)", //
6584+
"2.2572053268208536550832560045233873972354192817399");
6585+
65836586
check("EllipticK(0.999999999999999990000000000000000)", //
65846587
"20.9582676515692789828830607366566");
65856588
// reducing the accuracy gives ComplexInfinity in MMA:
65866589
check("EllipticK(0.99999999999999999)", //
65876590
"20.958267651569278");
65886591

6589-
check("N(EllipticK(8/10), 50)", //
6590-
"2.2572053268208536550832560045233873972354192817399");
6591-
65926592
checkNumeric("EllipticK(2.5+I)", //
65936593
"1.1551450606569331+I*0.9528453714670536");
65946594

@@ -15997,8 +15997,8 @@ public void testDeterminePrecision() {
1599715997
@Test
1599815998
public void testN() {
1599915999
// issue #942
16000-
// check("Tan(Pi/2) // N", //
16001-
// "ComplexInfinity");
16000+
check("Tan(Pi/2) // N", //
16001+
"ComplexInfinity");
1600216002
// issue #937
1600316003
check("(x==-157079632679/100000000000) // N", //
1600416004
"x==-1.5708");
@@ -24208,11 +24208,11 @@ public void testSurd() {
2420824208
checkNumeric("Surd(-3,3)", //
2420924209
"-3^(1/3)");
2421024210
checkNumeric("N((-3)^(1/3))", //
24211-
"0.7211247851537043+I*1.2490247664834064");
24211+
"0.7211247851537042+I*1.2490247664834064");
2421224212
checkNumeric("Surd(-3,3)-(-3)^(1/3)", //
2421324213
"-(-3)^(1/3)-3^(1/3)");
2421424214
checkNumeric("Surd(-3.,3)-(-3)^(1/3)", //
24215-
"-2.1633743554611127+I*(-1.2490247664834064)");
24215+
"-2.1633743554611122+I*(-1.2490247664834064)");
2421624216
checkNumeric("Surd(-3,3)", //
2421724217
"-3^(1/3)");
2421824218
checkNumeric("Surd(-3.,3)", //

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

+3-3
Original file line numberDiff line numberDiff line change
@@ -918,10 +918,10 @@ public void testSystem040() {
918918
"2.3715183290419594E-16+I*3.872983346207417");
919919
checkNumeric(".5^.5", //
920920
"0.7071067811865476");
921-
// checkNumeric("N((-15)^(1/2))", //
922-
// "I*3.872983346207417");
923921
checkNumeric("N((-15)^(1/2))", //
924-
"2.3715183290419594E-16+I*3.872983346207417");
922+
"I*3.872983346207417");
923+
// checkNumeric("N((-15)^(1/2))", //
924+
// "2.3715183290419594E-16+I*3.872983346207417");
925925
checkNumeric("N(Sin(1/2))", //
926926
"0.479425538604203");
927927
checkNumeric("N(1/6*(I*44^(1/2)+2))", //

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

+3-3
Original file line numberDiff line numberDiff line change
@@ -74,10 +74,10 @@ public void testGegenbauerC() {
7474
"{ComplexInfinity,ComplexInfinity,1.17445+I*2.30854,ComplexInfinity}");
7575
check("GegenbauerC(2, 0.5)", //
7676
"-0.5");
77-
// checkNumeric("GegenbauerC(5,1/8,7) //N", //
78-
// "16839.531372070312");
7977
checkNumeric("GegenbauerC(5,1/8,7) //N", //
80-
"16839.53137207032");
78+
"16839.531372070312");
79+
// checkNumeric("GegenbauerC(5,1/8,7) //N", //
80+
// "16839.53137207032");
8181
checkNumeric("Table(GegenbauerC(10, x), {x, 1, 5})", //
8282
"{1/5,262087/5,22619537/5,457470751/5,4517251249/5}");
8383
check("GegenbauerC({1/3, 1/2}, 1/6, 0)", //

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

+18-18
Original file line numberDiff line numberDiff line change
@@ -393,15 +393,15 @@ public void testSolve() {
393393
"{{x->-a},{x->b}}");
394394

395395
// github #261 - JUnit test for Apfloat switching to complex Power calculation
396-
// checkNumeric("Solve(0.00004244131815783 == x^5 , x)", //
397-
// "{{x->-0.10802279680851234+I*0.07848315587546606},"//
398-
// + "{x->-0.10802279680851212+I*(-0.07848315587546606)},"//
399-
// + "{x->0.04126103682102799+I*(-0.1269884137508598)},"//
400-
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
401-
check("Solve(0.00004244131815783 == x^5 , x)", //
402-
"{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},"//
403-
+ "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
404-
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");
396+
checkNumeric("Solve(0.00004244131815783 == x^5 , x)", //
397+
"{{x->-0.10802279680851234+I*0.07848315587546606},"//
398+
+ "{x->-0.10802279680851212+I*(-0.07848315587546606)},"//
399+
+ "{x->0.04126103682102799+I*(-0.1269884137508598)},"//
400+
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
401+
// check("Solve(0.00004244131815783 == x^5 , x)", //
402+
// "{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},"//
403+
// + "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
404+
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");
405405

406406
// github #247
407407
check("Solve((k+3)/(4)==(k)/2,{k})", //
@@ -801,15 +801,15 @@ public void testNSolve() {
801801
"{{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}}");
802802

803803
// github #261 - JUnit test for Apfloat switching to complex Power calculation
804-
// checkNumeric("NSolve(0.00004244131815783 == x^5 , x)", //
805-
// "{{x->-0.10802279680851234+I*0.07848315587546606}," //
806-
// + "{x->-0.10802279680851212+I*(-0.07848315587546606)}," //
807-
// + "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
808-
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
809-
check("NSolve(0.00004244131815783 == x^5 , x)", //
810-
"{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},{x->0.04126103682102799+I*(-0.1269884137508598)},"
811-
//
812-
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");
804+
checkNumeric("NSolve(0.00004244131815783 == x^5 , x)", //
805+
"{{x->-0.10802279680851234+I*0.07848315587546606}," //
806+
+ "{x->-0.10802279680851212+I*(-0.07848315587546606)}," //
807+
+ "{x->0.04126103682102799+I*(-0.1269884137508598)}," //
808+
+ "{x->0.04126103682102799+I*0.1269884137508598},{x->0.13352351997496842}}");
809+
// check("NSolve(0.00004244131815783 == x^5 , x)", //
810+
// "{{x->-0.10802279680851234+I*0.07848315587546605},{x->-0.10802279680851212+I*(-0.07848315587546605)},{x->0.04126103682102799+I*(-0.1269884137508598)},"
811+
// //
812+
// + "{x->0.04126103682102799+I*0.1269884137508598},{x->0.1335235199749684}}");
813813

814814

815815
// github #247

0 commit comments

Comments
 (0)