Skip to content

Commit

Permalink
Improve BellB, BellY, BernoulliB / Javadoc
Browse files Browse the repository at this point in the history
  • Loading branch information
@axkr
axkr committed Nov 9, 2023
1 parent ff20144 commit 922098d
Show file tree
Hide file tree
Showing 5 changed files with 122 additions and 82 deletions.
105 changes: 45 additions & 60 deletions ...oid_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/NumberTheory.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,20 +133,13 @@ private static class BellB extends AbstractFunctionEvaluator {
* @param index an int number >= 0
* @return
*/
private static IInteger bellNumber(int index) {
private static IExpr bellNumber(int index) {
if (index < BELLB_25.length) {
return AbstractIntegerSym.valueOf(BELLB_25[index]);
}

// Sum[StirlingS2[n, k], {k, 0, n}]
IInteger sum = F.C1;
for (int ki = 0; ki < index; ki++) {
sum = sum.add(stirlingS2(index, F.ZZ(ki), ki));
if (sum.bitLength() > Config.MAX_BIT_LENGTH / 100) {
BigIntegerLimitExceeded.throwIt(Config.MAX_BIT_LENGTH / 100);
}
}
return sum;
return F.sum(k -> stirlingS2(index, k, k.toIntDefault()), 0, index, 1);
}

/**
Expand All @@ -169,20 +162,21 @@ private static IExpr bellBPolynomial(int n, IExpr z) {
return z;
}

IASTAppendable sum = F.PlusAlloc(n + 1);
for (int k = 0; k <= n; k++) {
sum.append(F.Times(F.StirlingS2(F.ZZ(n), F.ZZ(k)), F.Power(z, k)));
}
return sum;
return F.sum(k -> F.Times(F.StirlingS2(F.ZZ(n), k), F.Power(z, k)), 0, n + 1, 1);
}

@Override
public IExpr evaluate(final IAST ast, EvalEngine engine) {
try {
IExpr arg1 = ast.arg1();
if (arg1.isNegative()) {
// Non-negative machine-sized integer expected at position `2` in `1`
return Errors.printMessage(S.BellB, "intnm", F.list(ast, F.C1), engine);
}
int n = arg1.toIntDefault();
if (n < 0) {
if (arg1.isNumber()) {
// Non-negative machine-sized integer expected at position `2` in `1`
return Errors.printMessage(S.BellB, "intnm", F.list(ast, F.C1), engine);
}
}
Expand Down Expand Up @@ -265,53 +259,48 @@ private static class BernoulliB extends AbstractFunctionEvaluator {

@Override
public IExpr evaluate(IAST ast, EvalEngine engine) {
if (ast.arg1().isNegative()) {
// Non-negative machine-sized integer expected at position `2` in `1`.
return Errors.printMessage(S.BernoulliB, "intnm", F.List(ast, F.C1), engine);
}
if (ast.isAST1()) {
try {
int bn = ast.arg1().toIntDefault();
if (bn >= 0) {
return bernoulliNumber(bn);
}
IExpr temp = engine.evaluate(F.Subtract(ast.arg1(), F.C3));
if (temp.isIntegerResult() && temp.isPositiveResult() && temp.isEvenResult()) {
// http://fungrim.org/entry/a98234/
return F.C0;
}

} catch (RuntimeException rex) {
Errors.printMessage(S.BernoulliB, rex, engine);
int bn = ast.arg1().toIntDefault();
if (bn >= 0) {
return bernoulliNumber(bn);
}
IExpr temp = engine.evaluate(F.Subtract(ast.arg1(), F.C3));
if (temp.isIntegerResult() && temp.isPositiveResult() && temp.isEvenResult()) {
// http://fungrim.org/entry/a98234/
return F.C0;
}
return F.NIL;
}
if (ast.isAST2()) {
try {
IExpr n = ast.arg1();
IExpr x = ast.arg2();
int xInt = x.toIntDefault();
if (xInt != Integer.MIN_VALUE) {
if (xInt == 0) {
// http://fungrim.org/entry/a1d2d7/
return F.BernoulliB(ast.arg1());
}
if (xInt == 1 && n.isIntegerResult()) {
// http://fungrim.org/entry/829185/
return F.Times(F.Power(F.CN1, n), F.BernoulliB(n));
}
IExpr n = ast.arg1();
IExpr x = ast.arg2();
int xInt = x.toIntDefault();
if (xInt != Integer.MIN_VALUE) {
if (xInt == 0) {
// http://fungrim.org/entry/a1d2d7/
return F.BernoulliB(ast.arg1());
}
if (n.isInteger() && n.isNonNegativeResult()) {
if (x.isNumEqualRational(F.C1D2)) {
// http://fungrim.org/entry/03ee0b/
return F.Times(F.Subtract(F.Power(F.C2, F.Subtract(F.C1, n)), F.C1), F.BernoulliB(n));
}
int bn = n.toIntDefault();
if (bn >= 0) {
// http://fungrim.org/entry/555e10/
return F.sum(
k -> F.Times(F.Binomial(n, k), F.BernoulliB(F.Subtract(n, k)), F.Power(x, k)), 0,
bn);
}
if (xInt == 1 && n.isIntegerResult()) {
// http://fungrim.org/entry/829185/
return F.Times(F.Power(F.CN1, n), F.BernoulliB(n));
}
}
if (n.isInteger() && n.isNonNegativeResult()) {
if (x.isNumEqualRational(F.C1D2)) {
// http://fungrim.org/entry/03ee0b/
return F.Times(F.Subtract(F.Power(F.C2, F.Subtract(F.C1, n)), F.C1), F.BernoulliB(n));
}
int bn = n.toIntDefault();
if (bn >= 0) {
// http://fungrim.org/entry/555e10/
return F.sum(
k -> F.Times(F.Binomial(n, k), F.BernoulliB(F.Subtract(n, k)), F.Power(x, k)), 0,
bn);
}
} catch (RuntimeException rex) {
Errors.printMessage(S.BernoulliB, rex, engine);
}
}
return F.NIL;
Expand Down Expand Up @@ -1709,10 +1698,6 @@ private static IExpr divisorSigma(IExpr arg1, IInteger n) {
// general formula
IASTAppendable sum = F.PlusAlloc(size);
return sum.appendArgs(size, i -> F.Power(list.get(i), arg1));
// for (int i = 1; i < size; i++) {
// sum.append(F.Power(list.get(i), arg1));
// }
// return sum;
}
return F.NIL;
}
Expand Down Expand Up @@ -5089,9 +5074,9 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
}
// try to convert into a fractional number
return rationalize(arg1, epsilon, useConvergenceMethod).orElse(arg1);
} catch (Exception ex) {
} catch (RuntimeException rex) {
// ex.printStackTrace();
Errors.printMessage(S.Rationalize, ex, engine);
Errors.printMessage(S.Rationalize, rex, engine);
}

return F.NIL;
Expand Down
28 changes: 22 additions & 6 deletions ...rary/matheclipse-core/src/main/java/org/matheclipse/core/builtin/PolynomialFunctions.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1516,9 +1516,14 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
if (arg1.isInteger() && ast.arg2().isInteger()) {
int n = arg1.toIntDefault();
int k = ast.arg2().toIntDefault();
if (n < 0 || k < 0 || !ast.arg3().isList() || ast.arg3().isMatrix() != null) {
if (n < 0 || k < 0) {
return F.NIL;
}
IExpr list = ast.arg3().normal(false);
if (!list.isList()) {
return F.NIL;
}

if (n == 0 && k == 0) {
return F.C1;
}
Expand All @@ -1533,7 +1538,7 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) {
}
int max = n - k + 2;
if (max >= 0) {
return bellY(n, k, (IAST) ast.arg3(), ast, engine);
return bellY(n, k, (IAST) list, ast, engine);
}
}

Expand Down Expand Up @@ -2259,13 +2264,24 @@ private static IAST monomialListModulus(IExpr polynomial, List<IExpr> variablesL
}
}

public static IExpr bellY(int n, int k, IAST symbols, IAST ast, EvalEngine engine) {
/**
* Calculate the partial Bell polynomial recursively.
*
* @param n
* @param k
* @param symbols
* @param listOfVariables
* @param engine
* @return
*
*/
public static IExpr bellY(int n, int k, IAST symbols, IAST listOfVariables, EvalEngine engine) {
final int recursionLimit = engine.getRecursionLimit();
try {
if (recursionLimit > 0) {
int counter = engine.incRecursionCounter();
if (counter > recursionLimit) {
RecursionLimitExceeded.throwIt(counter, ast);
RecursionLimitExceeded.throwIt(counter, listOfVariables);
}
}
if (n == 0 && k == 0) {
Expand All @@ -2280,11 +2296,11 @@ public static IExpr bellY(int n, int k, IAST symbols, IAST ast, EvalEngine engin

int iterationLimit = engine.getIterationLimit();
if (iterationLimit >= 0 && iterationLimit <= max) {
IterationLimitExceeded.throwIt(max, ast);
IterationLimitExceeded.throwIt(max, listOfVariables);
}
for (int m = 1; m < max; m++) {
if ((m < symbols.size()) && !symbols.get(m).isZero()) {
IExpr bellY = bellY(n - m, k - 1, symbols, ast, engine);
IExpr bellY = bellY(n - m, k - 1, symbols, listOfVariables, engine);
if (bellY.isPlus()) {
bellY = bellY.mapThread(F.Times(a, F.Slot1, symbols.get(m)), 2);
} else {
Expand Down
38 changes: 31 additions & 7 deletions symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/F.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,7 @@
import org.matheclipse.core.eval.EvalAttributes;
import org.matheclipse.core.eval.EvalEngine;
import org.matheclipse.core.eval.exception.ASTElementLimitExceeded;
import org.matheclipse.core.eval.exception.BigIntegerLimitExceeded;
import org.matheclipse.core.eval.interfaces.AbstractCoreFunctionEvaluator;
import org.matheclipse.core.eval.interfaces.ICoreFunctionEvaluator;
import org.matheclipse.core.eval.util.BiIntFunction;
Expand Down Expand Up @@ -8155,7 +8156,8 @@ public static IAST Product(final IExpr expr, final IExpr iterationSpecification)
* @param to
* @return
*/
public static IExpr product(final Function<IExpr, IExpr> function, final int from, final int to) {
public static IExpr product(final Function<IInteger, IExpr> function, final int from,
final int to) {
return intProduct(function, from, to, 1);
}

Expand Down Expand Up @@ -9565,7 +9567,8 @@ public static IRational sumRational(final IntFunction<IRational> function, final
* @param iMax
* @return
*/
public static IExpr sum(final Function<IExpr, IExpr> function, final int iMin, final int iMax) {
public static IExpr sum(final Function<IInteger, IExpr> function, final int iMin,
final int iMax) {
return intSum(function, iMin, iMax, 1);
}

Expand All @@ -9580,8 +9583,11 @@ public static IExpr sum(final Function<IExpr, IExpr> function, final int iMin, f
* @param step
* @return
*/
public static IExpr intProduct(final Function<IExpr, IExpr> function, final int from,
public static IExpr intProduct(final Function<IInteger, IExpr> function, final int from,
final int to, final int step) {
if (from > to && step > 0) {
return F.C1;
}
IASTAppendable result = ast(S.Times, 15);
long numberOfLeaves = 0;
INumber number = F.C1;
Expand All @@ -9592,6 +9598,10 @@ public static IExpr intProduct(final Function<IExpr, IExpr> function, final int
IExpr temp = engine.evaluate(function.apply(ZZ(i)));
if (temp.isNumber()) {
number = number.times((INumber) temp);
if (number instanceof IInteger //
&& ((IInteger) number).bitLength() > Config.MAX_BIT_LENGTH / 100) {
BigIntegerLimitExceeded.throwIt(Config.MAX_BIT_LENGTH / 100);
}
} else {
numberOfLeaves += temp.leafCount() + 1;
if (numberOfLeaves >= Config.MAX_AST_SIZE / 2) {
Expand All @@ -9614,6 +9624,9 @@ public static IExpr intProduct(final Function<IExpr, IExpr> function, final int
* @return
*/
public static IExpr intSum(final IntFunction<IExpr> function, final int iMin, final int iMax) {
if (iMin > iMax) {
return F.C0;
}
IASTAppendable result = ast(S.Plus, 15);
int numberOfLeaves = 0;
EvalEngine engine = EvalEngine.get();
Expand All @@ -9624,6 +9637,10 @@ public static IExpr intSum(final IntFunction<IExpr> function, final int iMin, fi
IExpr temp = engine.evaluate(function.apply(i));
if (temp.isNumber()) {
number = number.plus((INumber) temp);
if (number instanceof IInteger //
&& ((IInteger) number).bitLength() > Config.MAX_BIT_LENGTH / 100) {
BigIntegerLimitExceeded.throwIt(Config.MAX_BIT_LENGTH / 100);
}
} else {
numberOfLeaves += temp.leafCount() + 1;
if (numberOfLeaves >= Config.MAX_AST_SIZE / 2) {
Expand All @@ -9647,7 +9664,7 @@ public static IExpr intSum(final IntFunction<IExpr> function, final int iMin, fi
* @param step
* @return
*/
public static IExpr intSum(final Function<IExpr, IExpr> function, final int from, final int to,
public static IExpr intSum(final Function<IInteger, IExpr> function, final int from, final int to,
final int step) {
return intSum(function, from, to, step, false);
}
Expand All @@ -9663,8 +9680,11 @@ public static IExpr intSum(final Function<IExpr, IExpr> function, final int from
* @param expand expand {@link S#Plus}, {@link S#Times}, {@link S#Power} subexpressions
* @return
*/
public static IExpr intSum(final Function<IExpr, IExpr> function, final int from, final int to,
public static IExpr intSum(final Function<IInteger, IExpr> function, final int from, final int to,
final int step, boolean expand) {
if (from > to && step > 0) {
return F.C0;
}
IASTAppendable result = ast(S.Plus, 15);
long numberOfLeaves = 0;
INumber number = F.C0;
Expand All @@ -9675,6 +9695,10 @@ public static IExpr intSum(final Function<IExpr, IExpr> function, final int from
IExpr temp = engine.evaluate(function.apply(ZZ(i)));
if (temp.isNumber()) {
number = number.plus((INumber) temp);
if (number instanceof IInteger //
&& ((IInteger) number).bitLength() > Config.MAX_BIT_LENGTH / 100) {
BigIntegerLimitExceeded.throwIt(Config.MAX_BIT_LENGTH / 100);
}
} else {
numberOfLeaves += temp.leafCount() + 1;
if (numberOfLeaves >= Config.MAX_AST_SIZE / 2) {
Expand All @@ -9701,7 +9725,7 @@ public static IExpr intSum(final Function<IExpr, IExpr> function, final int from
* @param iStep
* @return
*/
public static IExpr sum(final Function<IExpr, IExpr> function, final int iMin, final int iMax,
public static IExpr sum(final Function<IInteger, IExpr> function, final int iMin, final int iMax,
final int iStep) {
return intSum(function, iMin, iMax, iStep);
}
Expand All @@ -9716,7 +9740,7 @@ public static IExpr sum(final Function<IExpr, IExpr> function, final int iMin, f
* @param expand expand {@link S#Plus}, {@link S#Times}, {@link S#Power} subexpressions
* @return
*/
public static IExpr sum(final Function<IExpr, IExpr> function, final int iMin, final int iMax,
public static IExpr sum(final Function<IInteger, IExpr> function, final int iMin, final int iMax,
final int iStep, boolean expand) {
return intSum(function, iMin, iMax, iStep, expand);
}
Expand Down
17 changes: 9 additions & 8 deletions ..._android_library/matheclipse-core/src/main/java/org/matheclipse/core/interfaces/IAST.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1511,19 +1511,20 @@ default IAST most() {
}

/**
* Get the argument at index 1, if the <code>size() == 2</code> or the complete ast if the <code>
* size() > 2</code> (useful for ASTs with attribute <code>OneIdentity</code> for example for
* <code>Plus[]</code> you can call <code>getOneIdentity(F.C0)</code> or for <code>Times[]</code>)
* you can call <code>getOneIdentity(F.C1)</code>.
* Return the argument at index 1, if the <code>argSize() == 1</code>. Or return the complete ast
* if the <code>argSize() > 1</code> If the <code>argSize() == 0</code> return
* <code>defaultValue</code> (useful for ASTs with attribute <code>OneIdentity</code> for example
* for <code>Plus[]</code> you can call <code>getOneIdentity(F.C0)</code> or for
* <code>Times[]</code>) you can call <code>getOneIdentity(F.C1)</code>.
*
* @param defaultValue default value, if <code>size() < 2</code>.
* @return
*/
public IExpr oneIdentity(IExpr defaultValue);

/**
* Get the argument at index 1, if the <code>size() == 2</code> or the complete ast if the <code>
* size() > 2</code> If the <code>size() == 1</code> return <code>0</code>.
* Return the argument at index 1, if the <code>argSize() == 1</code>. Or return the complete ast
* if the <code>argSize() > 1</code> If the <code>argSize() == 0</code> return <code>0</code>.
*
* @return
*/
Expand All @@ -1532,8 +1533,8 @@ default IExpr oneIdentity0() {
}

/**
* Get the argument at index 1, if the <code>size() == 2</code> or the complete ast if the <code>
* size() > 2</code> If the <code>size() == 1</code> return <code>1</code>.
* Return the argument at index 1, if the <code>argSize() == 1</code>. Or return the complete ast
* if the <code>argSize() > 1</code> If the <code>argSize() == 0</code> return <code>1</code>.
*
* @return
*/
Expand Down
Loading

0 comments on commit 922098d

Please sign in to comment.