Skip to content

Commit

Permalink
refactor ProductLog - avoid StackOverflow because of recursion
Browse files Browse the repository at this point in the history
- introduce new EvalEngine#evalNumericFunction() method
  • Loading branch information
axkr committed Mar 15, 2024
1 parent 6413070 commit 461ba62
Show file tree
Hide file tree
Showing 4 changed files with 70 additions and 18 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -2120,12 +2120,55 @@ public void setUp(final ISymbol newSymbol) {

@Override
public IExpr e1ObjArg(final IExpr o) {
IExpr temp = functionExpandLogArg(o);
if (temp.isPresent()) {
return temp;
}
if (o.equals(F.C0) || o.equals(F.CD0)) {
return F.C0;
}
return F.NIL;
}

private static IExpr functionExpandLogArg(final IExpr o) {
if (o.isTimes() && o.first().isFraction() && o.argSize() == 3) {
// ProductLog(Rational(k_,n_)*b_^Rational(c_,n_)*Log(b_)) :=
// Module( {a, v},
// a = N( (n*ProductLog((b^(c/n)*k*Log(b))/n))/Log(b) );
// v = Rationalize(a);
// v*Log(b)/n
// /; IntegerQ(v) && v >= 1 && PossibleZeroQ( (((-b^(c/n))*k + b^(v/n)*v)*Log(b))/n ))
IAST times = (IAST) o;
if (times.arg2().isPower() && times.arg3().isLog()) {
IAST power = (IAST) times.arg2();
if (power.base().isInteger() && power.base().equals(times.arg3().first())
&& power.exponent().isFraction()) {
IInteger b = (IInteger) times.arg3().first();
IFraction arg1 = (IFraction) times.arg1();
IInteger k = arg1.numerator();
IInteger n = arg1.denominator();
IFraction powExponent = (IFraction) power.exponent();
if (n.equals(powExponent.denominator())) {
EvalEngine engine = EvalEngine.get();
IExpr a = engine
.evalNumericFunction(F.Times(n, F.ProductLog(F.Times(arg1, F.Power(b, powExponent), F.Log(b))),
F.Power(F.Log(b), F.CN1)));
IExpr v = engine.evaluate(F.Rationalize(a));
if (v.isInteger() && ((IInteger) v).isGE(F.C1)) {
IFraction resultFactor = F.QQ((IInteger) v, n);
IExpr isZero = engine.evaluate(F.Plus(F.Times(k.negate(), F.Power(b, powExponent)),
F.Times(v, F.Power(b, resultFactor))));
if (isZero.isZero()) {
return engine.evaluate(F.Times(resultFactor, F.Log(b)));
}
}
}
}
}
}
return F.NIL;
}

@Override
public IExpr e2ObjArg(IExpr k, IExpr z) {
int ki = Integer.MIN_VALUE;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -787,8 +787,8 @@ private IAST endTrace() {
*
* @param result0 store the result of the evaluation in the i-th argument of the ast in <code>
* result0[0]</code>. <code>result0[0]</code> should be <code>F.NIL</code> if no evaluation
* occured.
* @param ast the original <code>ast</code> for whixh the arguments should be evaluated
* occurred.
* @param ast the original <code>ast</code> for which the arguments should be evaluated
* @param arg the i-th argument of <code>ast</code>
* @param i <code>arg</code> is the i-th argument of <code>ast</code>
* @param isNumericFunction if <code>true</code> the <code>NumericFunction</code> attribute is set
Expand Down Expand Up @@ -1671,14 +1671,14 @@ public final int evalInt(final IExpr expr) throws ArgumentTypeException {
result = expr.toIntDefault();
}
if (expr.isNumericFunction(true)) {
IExpr numericResult = evalN(expr);
IExpr numericResult = evalNumericFunction(expr);
if (numericResult.isReal()) {
result = numericResult.toIntDefault();
}
} else {
IExpr temp = evaluateNIL(expr);
if (temp.isNumericFunction(true)) {
IExpr numericResult = evalN(temp);
IExpr numericResult = evalNumericFunction(temp);
if (numericResult.isReal()) {
result = numericResult.toIntDefault();
}
Expand Down Expand Up @@ -2093,6 +2093,26 @@ private void printOnOffTrace(IExpr unevaledExpr, IExpr evaledExpr) {
}
}

public final IExpr evalNumericFunction(final IExpr expr) {
if (expr.isNumericFunction(true)) {
final boolean oldNumericMode = isNumericMode();
final long oldDigitPrecision = getNumericPrecision();
final int oldSignificantFigures = getSignificantFigures();
try {
setNumericMode(true, oldDigitPrecision, oldSignificantFigures);
IExpr temp = evalWithoutNumericReset(expr);
if (temp.isListOrAssociation() || temp.isRuleAST()) {
return ((IAST) temp).mapThread(arg -> evalNumericFunction(arg));
}
return temp;
} finally {
setNumericMode(oldNumericMode);
setNumericPrecision(oldDigitPrecision);
}
}
return expr;
}

/**
* Evaluates <code>expr</code> numerically.
*
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ public class ProductLogRules {
* <li>index 0 - number of equal rules in <code>RULES</code></li>
* </ul>
*/
final public static int[] SIZES = { 12, 5 };
final public static int[] SIZES = { 12, 4 };

final public static IAST RULES = List(
IInit(ProductLog, SIZES),
Expand Down Expand Up @@ -64,9 +64,6 @@ public class ProductLogRules {
Condition(x,GreaterEqual(x,CN1))),
// ProductLog(-1,E^x_*x_):=x/;x<=-1
ISetDelayed(ProductLog(CN1,Times(Exp(x_),x_)),
Condition(x,LessEqual(x,CN1))),
// ProductLog(Log(b_)*k_/n_*b_^(c_/n_)):=Module({a,v},a=N((n*ProductLog((b^(c/n)*k*Log(b))/n))/Log(b));v=Rationalize(a);v*Log(b)/n/;IntegerQ(v)&&v>=1&&PossibleZeroQ(((-b^(c/n)*k+b^(v/n)*v)*Log(b))/n))
ISetDelayed(ProductLog(Times(Log(b_),Rational(k_,n_),Power(b_,Rational(c_,n_)))),
Module(list(a,v),CompoundExpression(Set(a,N(Times(n,Power(Log(b),CN1),ProductLog(Times(Power(b,Times(c,Power(n,CN1))),k,Power(n,CN1),Log(b)))))),Set(v,Rationalize(a)),Condition(Times(v,Power(n,CN1),Log(b)),And(IntegerQ(v),GreaterEqual(v,C1),PossibleZeroQ(Times(Power(n,CN1),Plus(Times(CN1,Power(b,Times(c,Power(n,CN1))),k),Times(Power(b,Times(Power(n,CN1),v)),v)),Log(b))))))))
Condition(x,LessEqual(x,CN1)))
);
}
10 changes: 1 addition & 9 deletions symja_android_library/rules/ProductLogRules.m
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,5 @@
ProductLog(x_*E^x_) := x
/; x>=-1,
ProductLog(-1, x_*E^x_) := x
/; x<=-1,
ProductLog(Rational(k_,n_)*b_^Rational(c_,n_)*Log(b_)) :=
Module( {a, v},
a = N( (n*ProductLog((b^(c/n)*k*Log(b))/n))/Log(b) );
v = Rationalize(a);
v*Log(b)/n
/; IntegerQ(v) &&
v >= 1 &&
PossibleZeroQ( (((-b^(c/n))*k + b^(v/n)*v)*Log(b))/n ))
/; x<=-1
}

0 comments on commit 461ba62

Please sign in to comment.