From d56f0842a0d796eae65eb58b3144ef838e3a658a Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 12:30:49 +0200 Subject: [PATCH 1/8] SeriesCoefficient: closed Binomial form for (a+b*x)^p with symbolic index SeriesCoefficient[(a+b*x)^p, {x, 0, n}] with a symbolic index n emitted an opaque DifferenceRoot for finite binomial powers, and returned unevaluated for fractional exponents such as Sqrt(1+x). Add a binomialPowerSeriesCoefficient fast-path (run after the rational path and before the holonomic/DifferenceRoot path) that emits the closed form Piecewise({{b^n*a^(p-n)*Binomial(p,n), cond}}, 0), with cond = 0<=n<=p for a non-negative integer p and n>=0 otherwise. Matches Mathematica/WOXI. Guarded to symbolic index, x0==0, nonzero constant term, and non-negative/non-integer exponents (negative integers stay on the rational path for a cleaner form). Add testSeriesCoefficientBinomialWoxi to SeriesTest. Co-Authored-By: Claude Opus 4.8 --- .../core/builtin/SeriesFunctions.java | 53 +++++++++++++++++++ .../matheclipse/core/system/SeriesTest.java | 33 ++++++++++++ 2 files changed, 86 insertions(+) diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java index 4ab9f81a98..86cee03342 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java @@ -1000,6 +1000,13 @@ private static IExpr functionCoefficient(final IAST ast, IExpr function, IExpr x return rationalCoeff; } + // Binomial Power Fast-Path: (a + b*x)^p with a symbolic index yields a closed Binomial form + // (avoids the opaque DifferenceRoot the holonomic engine would otherwise emit). + IExpr binomialCoeff = binomialPowerSeriesCoefficient(function, x, x0, n, engine); + if (binomialCoeff.isPresent()) { + return binomialCoeff; + } + IExpr temp = polynomialSeriesCoefficient(function, x, x0, n, ast, engine); if (temp.isPresent()) { return temp; @@ -1179,6 +1186,52 @@ private static IExpr functionCoefficient(final IAST ast, IExpr function, IExpr x return taylorCoefficient(function, x, x0, n, engine); } + /** + * Series coefficient of a binomial power (a + b*x)^p with a symbolic index + * n, expanded around x==0. + * + *

+ * The coefficient of x^n in (a + b*x)^p is + * b^n * a^(p-n) * Binomial(p, n). When p is a non-negative integer the + * expansion is finite, so the index is bounded by 0 <= n <= p; otherwise the + * series runs for all n >= 0. Returning this closed Binomial form is a + * readability improvement over the equivalent (but opaque) DifferenceRoot the + * holonomic engine would otherwise emit for a finite binomial power. + * + * @return a {@link S#Piecewise} expression, or {@link F#NIL} if function is not a + * binomial power with a symbolic index expanded around zero + */ + private static IExpr binomialPowerSeriesCoefficient(IExpr function, IExpr x, IExpr x0, IExpr n, + EvalEngine engine) { + if (n.isNumber() || !x0.isZero() || !function.isPower()) { + return F.NIL; + } + IExpr p = function.exponent(); + if (!p.isFree(x)) { + return F.NIL; + } + // Negative integer exponents are rational functions handled by rationalSeriesCoefficient. + if (p.isInteger() && p.isNegative()) { + return F.NIL; + } + IExpr[] linear = function.base().linear(x); + if (linear == null) { + return F.NIL; + } + IExpr a = linear[0]; + IExpr b = linear[1]; + // a==0 would leave the messy 0^(p-n) factor; b==0 means the base is free of x. + if (a.isZero() || b.isZero()) { + return F.NIL; + } + IExpr coeff = engine + .evaluate(F.Times(F.Power(b, n), F.Power(a, F.Subtract(p, n)), F.Binomial(p, n))); + IExpr condition = p.isInteger() // + ? F.LessEqual(F.C0, n, p) // + : F.GreaterEqual(n, F.C0); + return F.Piecewise(F.list(F.list(coeff, condition)), F.C0); + } + /** * Constructs a DifferenceRoot (Holonomic sequence) representation for the series coefficients * of rational functions P(x)/Q(x) when 'n' is symbolic, and matches known identities (like diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java index 49b9740431..55d888fd9c 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java @@ -654,6 +654,39 @@ public void testSeriesCoefficient() { + "0}},0)"); } + @Test + public void testSeriesCoefficientBinomialWoxi() { + // Finite binomial powers (a+b*x)^p with a symbolic index n now return a closed Binomial + // Piecewise instead of an opaque DifferenceRoot (matches Mathematica / WOXI). + check("SeriesCoefficient((1+x)^5, {x, 0, n})", // + "Piecewise({{Binomial(5,n),0<=n<=5}},0)"); + check("SeriesCoefficient((1+2*x)^3, {x, 0, n})", // + "Piecewise({{2^n*Binomial(3,n),0<=n<=3}},0)"); + check("SeriesCoefficient((x+2)^3, {x, 0, n})", // + "Piecewise({{2^(3-n)*Binomial(3,n),0<=n<=3}},0)"); + + // Fractional exponent: infinite series bounded by n>=0 (previously returned unevaluated). + check("SeriesCoefficient(Sqrt(1+x), {x, 0, n})", // + "Piecewise({{Binomial(1/2,n),n>=0}},0)"); + + // Symbolic exponent p (not known to be a non-negative integer) => n>=0. + check("SeriesCoefficient((1+x)^p, {x, 0, n})", // + "Piecewise({{Binomial(p,n),n>=0}},0)"); + // Symbolic constant term in the base. + check("SeriesCoefficient((c+x)^3, {x, 0, n})", // + "Piecewise({{c^(3-n)*Binomial(3,n),0<=n<=3}},0)"); + + // A concrete index still returns the plain numeric coefficient (no Piecewise). + check("SeriesCoefficient((1+x)^5, {x, 0, 3})", // + "10"); + check("SeriesCoefficient((1+x)^5, {x, 0, 7})", // + "0"); + + // Negative-integer exponents stay on the rational-function path (cleaner (-1)^n form). + check("SeriesCoefficient((1+x)^(-1), {x, 0, n})", // + "Piecewise({{(-1)^n,n>=0}},0)"); + } + @Test public void testSeriesCoefficientDLMF() { check("SeriesCoefficient(BesselK(4,x), {x, 0, 3})", // From 12c5ef7e5caf17d0f2c73f9479c92841e64dea45 Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 14:56:50 +0200 Subject: [PATCH 2/8] Fix and test ZTransform, Laplace, Pade and discrete-difference functions Mine the Apache-licensed WOXI test suite for Symja calculus gaps, fix the stragglers found, and add WOXI-derived JUnit tests throughout. ZTransform/InverseZTransform: fix the n^k infinite recursion (guard the ExpandAll branch, route pure powers through the (-z d/dz)^k operator), combine polynomial results via Together, add a robust single-linear-pole inverse handler covering every sign/scaling variant and symbolic poles, and support UnitStep. PadeApproximant: when the Series is flagged Puiseux but expands into integer powers of x (e.g. Sqrt(1+x)), compute the approximant from the Normal Taylor polynomial instead of returning unevaluated. DifferenceDelta: canonicalize the generic forward difference (Factor for a numeric step, Expand for a symbolic step), leaving symbolic-function differences such as f(x+1)-f(x) untouched. DiscreteShift: for integer shifts, combine rational summands over a common denominator and expand top-level Plus results. DiscreteRatio: support the {k, order, step} spec and reduce the ratio via Together(FunctionExpand(...)). Co-Authored-By: Claude Opus 4.8 --- .../core/builtin/SeriesFunctions.java | 12 +- .../reflection/system/DifferenceDelta.java | 24 ++- .../core/reflection/system/DiscreteRatio.java | 55 ++--- .../core/reflection/system/DiscreteShift.java | 27 ++- .../reflection/system/InverseZTransform.java | 118 +++++++++++ .../core/reflection/system/ZTransform.java | 77 +++++-- .../reflection/system/ZTransformTest.java | 188 +++++++++++++++++- .../core/system/DifferenceDeltaTest.java | 36 ++++ .../core/system/LaplaceTransformTest.java | 99 +++++++++ .../core/system/LowercaseTestCase.java | 67 +++++++ .../matheclipse/core/system/SeriesTest.java | 149 ++++++++++++++ 11 files changed, 796 insertions(+), 56 deletions(-) diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java index 86cee03342..00c3b38657 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/builtin/SeriesFunctions.java @@ -142,7 +142,17 @@ private static IExpr seriesPade(IExpr function, IExpr x, IExpr x0, int m, int n, } ASTSeriesData series = (ASTSeriesData) seriesExpr; if (series.puiseuxDenominator() != 1) { - // a Puiseux series has fractional exponents which a rational approximant can't represent + // The series is flagged Puiseux (e.g. because of a Sqrt), but functions such as + // Sqrt(1+x) or Cos(Sqrt(x)) still expand into integer powers of x. When the normal + // form is an ordinary polynomial in x, the rational approximant is well defined, so + // compute it from that truncated Taylor polynomial (whose coefficients are exactly the + // ones a Pade approximant of order {m, n} consumes). + IExpr normal = engine.evaluate(F.Normal(seriesExpr)); + if (!normal.equals(function) && engine.evaluate(F.PolynomialQ(normal, x)).isTrue()) { + return engine.evaluate( + F.PadeApproximant(normal, F.List(x, x0, F.List(F.ZZ(m), F.ZZ(n))))); + } + // a genuine Puiseux series has fractional exponents a rational approximant can't represent return F.NIL; } int poleOrder = series.minExponent() < 0 ? -series.minExponent() : 0; diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DifferenceDelta.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DifferenceDelta.java index 8ed8ac9401..44e2fbd911 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DifferenceDelta.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DifferenceDelta.java @@ -98,7 +98,7 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { F.xreplace(arg1, symbolX, F.Plus(F.Times(F.ZZ(i), stepH), symbolX)); result.append(F.Times(F.Power(F.CN1, n - i), F.Binomial(n, i), diffTerm)); } - return result; + return normalizeDifference(result, symbolX, stepH, engine); } } } else { @@ -118,13 +118,33 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { } IExpr f2 = F.xreplace(arg1, arg2, F.Plus(F.C1, arg2)); - return F.Subtract(f2, arg1); + return normalizeDifference(F.Subtract(f2, arg1), arg2, F.C1, engine); } // All quantities that do not explicitly depend on the variables given are taken to have zero // partial difference. return F.C0; } + /** + * Bring a raw forward-difference result into the canonical output form. A genuine polynomial or + * rational function of x is canonicalized: with a numeric step it is factored + * (matching wolframscript, e.g. x^2 + x -> 2*(1 + x), and rational summands combine + * over a common denominator), with a symbolic step it is only expanded. A difference of unknown + * functions such as f(x + 1) - f(x) is already in simplest form and is returned + * unchanged (Factor must not reorder it or fold it into an opaque product). + */ + private static IExpr normalizeDifference(IExpr result, IExpr x, IExpr stepH, EvalEngine engine) { + IExpr together = engine.evaluate(F.Together(result)); + boolean rational = engine.evaluate(F.PolynomialQ(F.Numerator(together), x)).isTrue() + && engine.evaluate(F.PolynomialQ(F.Denominator(together), x)).isTrue(); + if (!rational) { + return result; + } + return stepH.isNumber() // + ? engine.evaluate(F.Factor(result)) + : engine.evaluate(F.ExpandAll(result)); + } + /** * Attempts to generate a closed-form solution for the n-th order forward difference of known * functions. diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java index 4786da7e6f..0b83f8801d 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java @@ -4,6 +4,7 @@ import org.matheclipse.core.eval.interfaces.AbstractFunctionEvaluator; import org.matheclipse.core.expression.F; import org.matheclipse.core.expression.ImplementationStatus; +import org.matheclipse.core.expression.S; import org.matheclipse.core.interfaces.IAST; import org.matheclipse.core.interfaces.IExpr; import org.matheclipse.core.interfaces.ISymbol; @@ -27,15 +28,21 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { ISymbol variable = null; int order = 1; + IExpr step = F.C1; - // Parse argument: either "k" or "{k, order}" + // Parse argument: "k", "{k, order}" or "{k, order, step}" if (arg.isSymbol()) { variable = (ISymbol) arg; } else if (arg.isList()) { IAST list = (IAST) arg; - if (list.isAST2() && list.arg1().isSymbol()) { + if (list.size() >= 2 && list.arg1().isSymbol()) { variable = (ISymbol) list.arg1(); - order = list.arg2().toIntDefault(); + if (list.size() >= 3) { + order = list.arg2().toIntDefault(); + } + if (list.size() >= 4) { + step = list.arg3(); + } } } @@ -46,38 +53,32 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { continue; } - // Apply ratio logic 'order' times + // Apply the ratio operator R[g] = g(k + step) / g(k) 'order' times for (int k = 0; k < order; k++) { - result = computeSingleRatio(result, variable, engine); - // If result becomes exactly 0 or 1, we can often stop or optimize, - // but for safety we continue (e.g. 0/0 cases typically don't occur here structurally). + result = computeSingleRatio(result, variable, step, engine); } } - return result; + // Reduce to wolframscript's canonical single-fraction ratio form. FunctionExpand collapses + // factorial/Gamma/Pochhammer/Binomial ratios (e.g. DiscreteRatio(n!, n) -> 1 + n), and + // Together cancels common polynomial factors (e.g. DiscreteRatio(n^2 + n, n) -> (2 + n)/n) + // while keeping a single grouped fraction (Cancel would rewrite (1 + n)^2/n^2 as (1 + 1/n)^2). + // Only adopt the FunctionExpand result when it fully reduces (no Gamma/Factorial residue left, + // e.g. a scaled argument like (2 n)! that woxi/Symja cannot reduce further). + IExpr expanded = engine.evaluate(F.Together(F.FunctionExpand(result))); + if (expanded.isFree(S.Gamma, true) && expanded.isFree(S.Factorial, true)) { + return expanded; + } + return engine.evaluate(F.Together(result)); } /** - * Calculates Simplify( f(k) / f(k-1) ) + * Compute the single-step ratio f(k + step) / f(k). */ - private IExpr computeSingleRatio(IExpr expr, ISymbol variable, EvalEngine engine) { - // Calculate denominator: f(k-1) - // We use subst to substitute k -> 1+k - IExpr numerator = F.subst(expr, variable, F.Plus(F.C1, variable)); - - IExpr fraction = F.Divide(numerator, expr); - - // Standard evaluate() is not enough to cancel terms like x^n / x^(n-1). - // F.Cancel divides out common factors from numerator and denominator. - // F.Simplify is more powerful but slower; usually Cancel is sufficient for Ratios. - - // We try Cancel first, as it handles the polynomial/multiplicative cancellation best. - IExpr cancelled = engine.evaluate(F.FunctionExpand(fraction)); - - // Fallback: If Cancel didn't reduce it effectively (e.g. Gamma functions), - // Simplify or FunctionExpand might be needed. - // For now, Cancel is the standard approach for rational simplification. - return cancelled; + private static IExpr computeSingleRatio(IExpr expr, ISymbol variable, IExpr step, + EvalEngine engine) { + IExpr shifted = F.subst(expr, variable, F.Plus(variable, step)); + return engine.evaluate(F.Divide(shifted, expr)); } @Override diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteShift.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteShift.java index 68cd6a0f5f..f17e0c4276 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteShift.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteShift.java @@ -31,6 +31,7 @@ public IExpr evaluate(IAST ast, EvalEngine engine) { // Use IASTAppendable to collect the transformation rules. // We allocate size based on ast.size() - 1 args. IASTAppendable rulesList = F.ListAlloc(ast.argSize()); + boolean allIntegerShifts = true; for (int i = 2; i < ast.size(); i++) { IExpr arg = ast.get(i); @@ -67,14 +68,34 @@ public IExpr evaluate(IAST ast, EvalEngine engine) { return F.NIL; } + if (!engine.evaluate(shift).isInteger()) { + allIntegerShifts = false; + } // Create the rule: variable -> variable + shift // Example: k -> k + 1 or k -> k + m rulesList.append(F.Rule(variable, F.Plus(variable, shift))); } - // Apply all rules to the expression using F.subst - // This performs the structural shift without invoking Simplification logic - return F.subst(expression, rulesList); + // Apply all rules to the expression structurally. + IExpr result = engine.evaluate(F.subst(expression, rulesList)); + // With an integer shift wolframscript combines rational summands over a common denominator, so + // DiscreteShift(1/(2 n + 1), n) becomes (3 + 2 n)^-1 rather than (1 + 2 (1 + n))^-1. A symbolic + // shift stays unfolded (e.g. (1 + 2 (k + n))^-1). A pure polynomial product such as + // (1 + m) (1 + n) has no denominator to combine and must not be expanded here. + boolean hasDenominator = !result.isFree(x -> x.isPower() && x.exponent().isNegative(), false); + if (allIntegerShifts && hasDenominator) { + IExpr together = engine.evaluate(F.Together(result)); + IExpr denominator = engine.evaluate(F.Denominator(together)); + // Together leaves a lone denominator such as 1 + 2 (1 + n) unexpanded, so expand it. + result = denominator.isOne() // + ? together + : engine.evaluate(F.Divide(F.Numerator(together), F.Expand(denominator))); + } + // wolframscript expands the result only when its top-level head is Plus. + if (result.isPlus()) { + result = engine.evaluate(F.Expand(result)); + } + return result; } @Override diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/InverseZTransform.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/InverseZTransform.java index 27445bb6f0..e97147d684 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/InverseZTransform.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/InverseZTransform.java @@ -109,6 +109,18 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { } } + // ======================================================================== + // Single linear-pole base case (geometric and polynomial-geometric): + // any F(z) whose F(z)/z has the shape numerator(z) / (b*z + a)^p is inverted + // in closed form. Handles every sign/scaling variant of the pole (including a + // symbolic pole such as (a - z)), which neither the pattern rules in + // InverseZTransformRules.m nor Apart (for symbolic poles) can resolve. + // ======================================================================== + IExpr geometric = inverseLinearPole(fx, z, n, engine); + if (geometric.isPresent()) { + return geometric; + } + if (fx.isTimes()) { IAST function = (IAST) fx; IASTAppendable constantArgs = F.TimesAlloc(); @@ -167,6 +179,112 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { return F.NIL; } + /** + * Inverse Z transform of any F(z) whose F(z)/z is a proper rational + * function with a single linear pole, i.e. numerator(z) / (b*z + a)^p where + * b, a are free of z and p is a positive + * integer. + * + *

+ * The numerator is re-expanded in powers of the pole L = b*z + a via the + * substitution z -> (w - a)/b, splitting F(z)/z into simple partial + * fractions d_j * L^(j-p). Multiplying back by z maps each term through + * the elementary pair + * + * 

+   * Z^-1{ z / (b*z + a)^k } = b^(-k) * Binomial(n, k-1) * r^(n-k+1),   r = -a/b
+   *
+ * + * This covers every sign and scaling variant of the pole (e.g. z/(z-1)^2, + * z/(1-z)^2, (3*z)/(-1+3*z), -(z/(a-z))) and non-trivial + * numerators over a symbolic pole (e.g. -(a*z*(a+z))/(a-z)^3), which neither the + * fixed pattern rules nor {@link org.matheclipse.core.expression.S#Apart} (for symbolic poles) + * resolve. + * + * @return the inverse transform, or {@link F#NIL} if fx is not of this shape + */ + private static IExpr inverseLinearPole(IExpr fx, IExpr z, IExpr n, EvalEngine engine) { + if (fx.isFree(z)) { + return F.NIL; + } + // Read numerator/denominator from the raw input so a factored pole such as (z - 1)^3 + // stays factored (Together would expand it into -1 + 3*z - 3*z^2 + z^3). + IExpr num = engine.evaluate(F.Numerator(fx)); + IExpr den = engine.evaluate(F.Denominator(fx)); + if (den.isOne()) { + return F.NIL; + } + // numG = numerator(F(z)/z): F(z)/z is the standard quotient used in the inverse transform, + // and a genuine Z-transform carries a factor z in the numerator that must cancel cleanly. + IExpr numG = engine.evaluate(F.Cancel(F.Divide(num, z))); + if (!engine.evaluate(F.PolynomialQ(numG, z)).isTrue()) { + return F.NIL; + } + // Recover a factored form if the denominator arrived expanded. + if (den.isPlus()) { + den = engine.evaluate(F.Factor(den)); + } + // Peel off any factor free of z (e.g. Factor may yield 9*(z - 1/3)^2), folding it into numG. + if (den.isTimes()) { + IASTAppendable constPart = F.TimesAlloc(); + IASTAppendable rest = F.TimesAlloc(); + for (IExpr factor : (IAST) den) { + (factor.isFree(z) ? constPart : rest).append(factor); + } + if (constPart.argSize() > 0 && rest.argSize() > 0) { + numG = engine.evaluate(F.Divide(numG, constPart.oneIdentity1())); + den = engine.evaluate(rest.oneIdentity1()); + } + } + // denominator must be linear(z)^p with p a positive integer + IExpr base; + int p; + if (den.isPower() && den.exponent().isInteger() && den.exponent().isPositive()) { + base = den.base(); + p = den.exponent().toIntDefault(); + } else { + base = den; + p = 1; + } + if (p < 1) { + return F.NIL; + } + // base must be linear in z: b*z + a (b, a free of z, a != 0 so the pole r = -a/b != 0) + IExpr b = engine.evaluate(F.Coefficient(base, z, F.C1)); + IExpr a = engine.evaluate(F.Coefficient(base, z, F.C0)); + if (b.isZero() || a.isZero() || !b.isFree(z) || !a.isFree(z)) { + return F.NIL; + } + if (!engine.evaluate(F.Subtract(base, F.Plus(F.Times(b, z), a))).isZero()) { + return F.NIL; + } + IExpr r = engine.evaluate(F.Divide(F.Negate(a), b)); + + // Re-expand numG as a polynomial in w = base, via z -> (w - a)/b: + // numG(z) = Sum_j d_j * base^j, so F(z)/z = Sum_j d_j * base^(j - p). + IExpr w = F.Dummy("w"); + IExpr q = engine.evaluate(F.ExpandAll(F.subst(numG, z, F.Divide(F.Subtract(w, a), b)))); + int degQ = engine.evaluate(F.Exponent(q, w)).toIntDefault(); + if (degQ < 0 || degQ >= p) { + // improper: a polynomial part would contribute DiscreteDelta terms - not handled here + return F.NIL; + } + IASTAppendable sum = F.PlusAlloc(degQ + 1); + for (int j = 0; j <= degQ; j++) { + IExpr dj = engine.evaluate(F.Coefficient(q, w, F.ZZ(j))); + if (dj.isZero()) { + continue; + } + long k = p - j; // k >= 1 because j <= degQ < p + // Z^-1{ z / base^k } = b^(-k) * Binomial(n, k-1) * r^(n-k+1) + sum.append(F.Times(dj, // + F.Power(b, F.ZZ(-k)), // + F.Binomial(n, F.ZZ(k - 1L)), // + F.Power(r, F.Plus(n, F.ZZ(1L - k))))); + } + return engine.evaluate(F.Simplify(F.FunctionExpand(sum))); + } + @Override public int[] expectedArgSize(IAST ast) { return IFunctionEvaluator.ARGS_3_3; diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/ZTransform.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/ZTransform.java index d2de2a939e..7f4bf6b843 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/ZTransform.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/ZTransform.java @@ -42,8 +42,18 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { if (expr.equals(n)) { return F.Times(F.Power(F.Plus(F.CN1, z), F.CN2), z); } + if (expr.isAST(S.UnitStep, 2) && expr.first().equals(n)) { + // UnitStep(n) == 1 for n >= 0, so ZTransform(UnitStep(n), n, z) == ZTransform(1, n, z) + return engine.evaluate(F.Divide(z, F.Subtract(z, F.C1))); + } if (expr.isPlus()) { - return ((IAST) expr).mapThread(F.ZTransform(F.Slot1, n, z), 1); + IExpr mapped = ((IAST) expr).mapThread(F.ZTransform(F.Slot1, n, z), 1); + // For polynomial inputs in n, combine the partial results into a single + // rational function (matches Mathematica's ZTransform output). + if (engine.evaluate(F.PolynomialQ(expr, n)).isTrue()) { + return engine.evaluate(F.Together(mapped)); + } + return mapped; } if (expr.isTimes()) { @@ -52,9 +62,21 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { return temp; } } else if (expr.isPower()) { - // RESTORED: Expand polynomial powers + IExpr base = expr.base(); + IExpr exponent = expr.exponent(); + if (base.equals(n) && exponent.isInteger() && exponent.isPositive()) { + // ZTransform(n^k, n, z) = (-z * d/dz)^k ZTransform(1, n, z) for integer k > 0 + int k = exponent.toIntDefault(); + if (k > 0) { + return polynomialTransform(engine, k, F.C1, n, z); + } + } + // Expand polynomial powers with a compound base, e.g. (1 + n)^2 if (!expr.isExpanded()) { - return engine.evaluate(F.ZTransform(engine.evaluate(F.ExpandAll(expr)), n, z)); + IExpr expanded = engine.evaluate(F.ExpandAll(expr)); + if (!expanded.equals(expr)) { + return engine.evaluate(F.ZTransform(expanded, n, z)); + } } } @@ -192,22 +214,7 @@ private IExpr times(EvalEngine engine, IAST timesAST, IExpr n, IExpr z) { } else { f_n = restNArgs; } - - // Use a secure dummy variable for differentiation. - // If 'z' is a compound fractional expression like (1/s), D(..., 1/s) is invalid! - IExpr zSym = z.isSymbol() ? z : F.Dummy("zDummy"); - IExpr result = engine.evaluate(F.ZTransform(f_n, n, zSym)); - - // Iteratively apply the (-z * D(..., z)) operator k times - for (int i = 0; i < kExponent; i++) { - result = engine.evaluate(F.Times(F.CN1, zSym, F.D(result, zSym))); - } - - // Safely project the solved derivative result back onto the original compound variable - if (!zSym.equals(z)) { - result = engine.evaluate(F.subst(result, zSym, z)); - } - return engine.evaluate(F.Simplify(result)); + return polynomialTransform(engine, kExponent, f_n, n, z); } // Expand polynomials multiplied by functions (e.g., (1+n)^2 * f(n)) @@ -223,6 +230,38 @@ private IExpr times(EvalEngine engine, IAST timesAST, IExpr n, IExpr z) { return F.NIL; } + /** + * Apply the polynomial-multiplication property + * ZTransform(n^k * f(n), n, z) == (-z * d/dz)^k ZTransform(f(n), n, z). + * + * @param k a positive integer exponent + * @param f_n the remaining factor f(n) (or {@link F#C1}) + */ + private IExpr polynomialTransform(EvalEngine engine, int k, IExpr f_n, IExpr n, IExpr z) { + // Use a secure dummy variable for differentiation. + // If 'z' is a compound fractional expression like (1/s), D(..., 1/s) is invalid! + IExpr zSym = z.isSymbol() ? z : F.Dummy("zDummy"); + IExpr result = engine.evaluate(F.ZTransform(f_n, n, zSym)); + + // Iteratively apply the (-z * D(..., z)) operator k times + for (int i = 0; i < k; i++) { + result = engine.evaluate(F.Times(F.CN1, zSym, F.D(result, zSym))); + } + + // Safely project the solved derivative result back onto the original compound variable + if (!zSym.equals(z)) { + result = engine.evaluate(F.subst(result, zSym, z)); + } + if (f_n.isOne()) { + // Pure power ZTransform(n^k, n, z) is a rational function of z. Together first so that + // higher-order derivatives (n^4, ...) collapse into a single fraction that Simplify alone + // would otherwise leave expanded. (Only for the pure-power case, to avoid reshuffling the + // symbolic D(ZTransform(f(n),...)) forms produced when f(n) is an unknown sequence.) + result = engine.evaluate(F.Together(result)); + } + return engine.evaluate(F.Simplify(result)); + } + @Override public int[] expectedArgSize(IAST ast) { return IFunctionEvaluator.ARGS_3_3; diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/ZTransformTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/ZTransformTest.java index 82a5ae753b..2a3a2afd8c 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/ZTransformTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/ZTransformTest.java @@ -108,9 +108,189 @@ public void testZTransform003() { "((a+b+c)*z)/(-1+z)"); } - // ========================================================== - // Transform Boundary Limit Tests - // ========================================================== + @Test + public void testZTransform004() { + check("ZTransform(n^2 + n + 1, n, z)", // + "(z+z^3)/(-1+z)^3"); + } + + @Test + public void testZTransformMonomials() { + check("ZTransform(1, n, z)", // + "z/(-1+z)"); + check("ZTransform(n, n, z)", // + "z/(1-z)^2"); + check("ZTransform(n^2, n, z)", // + "(z*(1+z))/(-1+z)^3"); + check("ZTransform(n^3, n, z)", // + "(z*(1+4*z+z^2))/(1-z)^4"); + check("ZTransform(n^4, n, z)", // + "(z*(1+11*z+11*z^2+z^3))/(-1+z)^5"); + } + + @Test + public void testZTransformConstantMultiples() { + check("ZTransform(2, n, z)", // + "(2*z)/(-1+z)"); + check("ZTransform(2*n, n, z)", // + "(2*z)/(1-z)^2"); + check("ZTransform(2*n^2, n, z)", // + "(2*z*(1+z))/(-1+z)^3"); + // Symbols free of n act as constants + check("ZTransform(x, n, z)", // + "(x*z)/(-1+z)"); + } + + @Test + public void testZTransformGeometric() { + check("ZTransform(a^n, n, z)", // + "z/(-a+z)"); + check("ZTransform(2^n, n, z)", // + "z/(-2+z)"); + // Rational base: z/(-1/3+z) == (3*z)/(-1+3*z) + check("ZTransform((1/3)^n, n, z)", // + "z/(-1/3+z)"); + } + + @Test + public void testZTransformPolynomialTimesGeometric() { + check("ZTransform(n*a^n, n, z)", // + "(a*z)/(-a+z)^2"); + check("ZTransform(n*2^n, n, z)", // + "(2*z)/(2-z)^2"); + check("ZTransform(n/3^n, n, z)", // + "(3*z)/(1-3*z)^2"); + check("ZTransform(n^2*a^n, n, z)", // + "(a*z*(a+z))/(-a+z)^3"); + check("ZTransform(n^2*2^n, n, z)", // + "(2*z*(2+z))/(-2+z)^3"); + // The documentation example: (z*(2+4*z))/(-1+2*z)^3 == (2*z*(1+2*z))/(-1+2*z)^3 + check("ZTransform(n^2/2^n, n, z)", // + "(z*(2+4*z))/(-1+2*z)^3"); + } + + @Test + public void testZTransformInverseFactorial() { + check("ZTransform(1/n!, n, z)", // + "E^(1/z)"); + check("ZTransform(3^n/n!, n, z)", // + "E^(3/z)"); + } + + @Test + public void testZTransformTrig() { + check("ZTransform(Sin(n), n, z)", // + "(z*Sin(1))/(1+z^2-2*z*Cos(1))"); + check("ZTransform(Cos(n), n, z)", // + "(z*(z-Cos(1)))/(1+z^2-2*z*Cos(1))"); + check("ZTransform(Sin(a*n), n, z)", // + "(z*Sin(a))/(1+z^2-2*z*Cos(a))"); + check("ZTransform(Cos(a*n), n, z)", // + "(z*(z-Cos(a)))/(1+z^2-2*z*Cos(a))"); + check("ZTransform(Sin(2*n), n, z)", // + "(z*Sin(2))/(1+z^2-2*z*Cos(2))"); + check("ZTransform(Sin(n/2), n, z)", // + "(z*Sin(1/2))/(1+z^2-2*z*Cos(1/2))"); + } + + @Test + public void testZTransformUnitStep() { + // UnitStep(n) == 1 for n >= 0, so its Z-transform is z/(z-1) + check("ZTransform(UnitStep(n), n, z)", // + "z/(-1+z)"); + } + + @Test + public void testZTransformUnsupported() { + // Sin(n^2) is not a linear a*n argument, so there is no closed form + check("ZTransform(Sin(n^2), n, z)", // + "ZTransform(Sin(n^2),n,z)"); + } + + @Test + public void testInverseZTransformGeometric() { + check("InverseZTransform(z/(z - a), z, n)", // + "a^n"); + check("InverseZTransform(z/(z - 1), z, n)", // + "1"); + check("InverseZTransform(z/(z - 2), z, n)", // + "2^n"); + // (1/3)^n == 3^(-n) + check("InverseZTransform((3*z)/(-1 + 3*z), z, n)", // + "(1/3)^n"); + // Sign-wrapped spelling of the same transform + check("InverseZTransform(-(z/(a - z)), z, n)", // + "a^n"); + } + + @Test + public void testInverseZTransformPolynomialSequences() { + check("InverseZTransform(z/(z - 1)^2, z, n)", // + "n"); + check("InverseZTransform((z*(1 + z))/(-1 + z)^3, z, n)", // + "n^2"); + } + + @Test + public void testInverseZTransformBinomialSequences() { + // 1/2*(-1+n)*n == ((-1 + n)*n)/2 + check("InverseZTransform(z/(z - 1)^3, z, n)", // + "1/2*(-1+n)*n"); + check("InverseZTransform(z/(z - 1)^4, z, n)", // + "1/6*(-2+n)*(-1+n)*n"); + } + + @Test + public void testInverseZTransformPolynomialTimesGeometric() { + check("InverseZTransform((a*z)/(a - z)^2, z, n)", // + "a^n*n"); + check("InverseZTransform((2*z)/(-2 + z)^2, z, n)", // + "2^n*n"); + check("InverseZTransform((3*z)/(1 - 3*z)^2, z, n)", // + "n/3^n"); + check("InverseZTransform((2*z*(2 + z))/(-2 + z)^3, z, n)", // + "2^n*n^2"); + check("InverseZTransform(-((a*z*(a + z))/(a - z)^3), z, n)", // + "a^n*n^2"); + check("InverseZTransform((2*z*(1 + 2*z))/(-1 + 2*z)^3, z, n)", // + "n^2/2^n"); + } + + @Test + public void testInverseZTransformExponentialForms() { + // Gamma(1+n) == n! + check("InverseZTransform(E^(1/z), z, n)", // + "1/Gamma(1+n)"); + check("InverseZTransform(E^(3/z), z, n)", // + "3^n/Gamma(1+n)"); + } + + @Test + public void testInverseZTransformConstants() { + check("InverseZTransform((2*z)/(-1 + z), z, n)", // + "2"); + check("InverseZTransform((x*z)/(-1 + z), z, n)", // + "x"); + // z-free input is a DiscreteDelta impulse + check("InverseZTransform(1, z, n)", // + "DiscreteDelta(n)"); + check("InverseZTransform(5, z, n)", // + "5*DiscreteDelta(n)"); + } + + @Test + public void testZTransformRoundTrip() { + check("InverseZTransform(ZTransform(n^2/2^n, n, z), z, n)", // + "n^2/2^n"); + check("ZTransform(InverseZTransform(z/(z - 2), z, n), n, z)", // + "z/(-2+z)"); + } + + @Test + public void testInverseZTransformUnsupported() { + check("InverseZTransform(Sin(z), z, n)", // + "InverseZTransform(Sin(z),z,n)"); + } @Test public void testZTransformFractionalShift() { @@ -129,7 +309,7 @@ public void testEGFTransformFractionalShift() { } @Test - public void testExponentialFractions() { + public void testInverseZTransformExponentialFractions() { // Z^-1 { E^(a/z) } check("InverseZTransform(E^(a/z), z, n)", // "a^n/Gamma(1+n)"); diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java index 52aef73566..fdc3839de0 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java @@ -71,4 +71,40 @@ public void testDifferenceDeltaIndefiniteSum() { check("DifferenceDelta(Sum(k*k!,k),k)", // "k*k!"); } + + // ========================================================== + // Ported from the Woxi project (tests/interpreter_tests/calculus.rs, mod difference_delta). + // The raw forward difference is now brought into canonical form: a numeric step is factored + // (x^2 + x -> 2*(1 + x), and rational summands combine over a common denominator), a symbolic + // step is expanded. Symja's form is sometimes cosmetically different but mathematically equal, + // e.g. 1/((-1-n)*n) == -1/(n*(1+n)). + // ========================================================== + @Test + public void testDifferenceDeltaWoxi() { + check("DifferenceDelta(5, x)", "0"); + check("DifferenceDelta(x, x)", "1"); + check("DifferenceDelta(a*x + b, x)", "a"); + check("DifferenceDelta(x^2, x)", "1+2*x"); + check("DifferenceDelta(x^3, x)", "1+3*x+3*x^2"); + check("DifferenceDelta(f(x), x)", "-f(x)+f(1+x)"); + // Wolfram simplifies this trig difference to 2*Sin[1/2]*Sin[(1+Pi)/2+x] + check("DifferenceDelta(Sin(x), x)", "2*Sin(1/2)*Sin(1/2*(1+Pi)+x)"); + // second-order difference of x^2 is 2 + check("DifferenceDelta(x^2, {x, 2})", "2"); + // zeroth-order difference returns the expression itself + check("DifferenceDelta(x^2, {x, 0})", "x^2"); + check("DifferenceDelta(x^2, {x, 1, h})", "h^2+2*h*x"); + check("DifferenceDelta(y, x)", "0"); + check("DifferenceDelta(2^x, x)", "2^x"); + // rational summands combine over a common denominator + check("DifferenceDelta(1/(2*n + 1), n)", "-2/((1+2*n)*(3+2*n))"); + // 1/((-1-n)*n) == -1/(n*(1+n)) + check("DifferenceDelta(1/n, n)", "1/((-1-n)*n)"); + // a numeric step factors the polynomial result + check("DifferenceDelta(x^2 + x, x)", "2*(1+x)"); + check("DifferenceDelta(a*x^2, x)", "a*(1+2*x)"); + check("DifferenceDelta(x^2, {x, 1, 2})", "4*(1+x)"); + // a symbolic step stays expanded (not factored) + check("DifferenceDelta(x^2 + x, {x, 1, h})", "h+h^2+2*h*x"); + } } diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java index 16aef9c278..0ac6469840 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java @@ -324,4 +324,103 @@ public void testLaplaceTransformEdgeCases() { "(a*b)/s"); } + // ========================================================== + // Ported from the Woxi project (laplace_transform module). Symja's canonical output + // form sometimes differs cosmetically from wolframscript / Woxi but is mathematically + // equal, e.g. 1/s == s^(-1), 1/(1+s^2) == (1 + s^2)^(-1), E^(-a*t) == E^(-(a*t)). + // ========================================================== + + @Test + public void testLaplaceTransformWoxiMonomials() { + // 1/s == s^(-1) + check("LaplaceTransform(1, t, s)", // + "1/s"); + // 1/s^2 == s^(-2) + check("LaplaceTransform(t, t, s)", // + "1/s^2"); + check("LaplaceTransform(t^2, t, s)", // + "2/s^3"); + check("LaplaceTransform(t^3, t, s)", // + "6/s^4"); + check("LaplaceTransform(5*t, t, s)", // + "5/s^2"); + } + + @Test + public void testLaplaceTransformWoxiTrig() { + // 1/(1+s^2) == (1 + s^2)^(-1) + check("LaplaceTransform(Sin(t), t, s)", // + "1/(1+s^2)"); + check("LaplaceTransform(Cos(t), t, s)", // + "s/(1+s^2)"); + check("LaplaceTransform(Sin(3*t), t, s)", // + "3/(9+s^2)"); + } + + @Test + public void testLaplaceTransformWoxiExponential() { + // 1/(a+s) == (a + s)^(-1) + check("LaplaceTransform(Exp(-a*t), t, s)", // + "1/(a+s)"); + // 1/(-a+s) == (-a + s)^(-1) + check("LaplaceTransform(Exp(a*t), t, s)", // + "1/(-a+s)"); + } + + @Test + public void testLaplaceTransformWoxiLinearity() { + check("LaplaceTransform(3*t^2 + 2*Sin(t), t, s)", // + "6/s^3+2/(1+s^2)"); + } + + @Test + public void testLaplaceTransformWoxiUnevaluated() { + // Unknown functions should return unevaluated + check("LaplaceTransform(BesselJ(0, t), t, s)", // + "LaplaceTransform(BesselJ(0,t),t,s)"); + } + + // ========================================================== + // InverseLaplaceTransform, ported from the Woxi project. + // ========================================================== + + @Test + public void testInverseLaplaceTransformWoxiMonomials() { + check("InverseLaplaceTransform(1/s, s, t)", // + "1"); + check("InverseLaplaceTransform(1/s^2, s, t)", // + "t"); + check("InverseLaplaceTransform(2/s^3, s, t)", // + "t^2"); + check("InverseLaplaceTransform(6/s^4, s, t)", // + "t^3"); + } + + @Test + public void testInverseLaplaceTransformWoxiTrig() { + check("InverseLaplaceTransform(1/(s^2 + 1), s, t)", // + "Sin(t)"); + check("InverseLaplaceTransform(s/(s^2 + 1), s, t)", // + "Cos(t)"); + check("InverseLaplaceTransform(a/(s^2 + a^2), s, t)", // + "Sin(a*t)"); + check("InverseLaplaceTransform(s/(s^2 + a^2), s, t)", // + "Cos(a*t)"); + } + + @Test + public void testInverseLaplaceTransformWoxiExponential() { + check("InverseLaplaceTransform(1/(s - a), s, t)", // + "E^(a*t)"); + // E^(-a*t) == E^(-(a*t)) + check("InverseLaplaceTransform(1/(s + a), s, t)", // + "E^(-a*t)"); + } + + @Test + public void testInverseLaplaceTransformWoxiUnevaluated() { + check("InverseLaplaceTransform(Log(s), s, t)", // + "InverseLaplaceTransform(Log(s),s,t)"); + } + } diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java index cf9040e77a..5e8c99237c 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java @@ -5517,6 +5517,16 @@ public void testDifferenceQuotient() { check("dq = DifferenceQuotient(f(x), {x, h})", // "(-x*Sin(x)+(h+x)*Sin(h+x))/h"); + // ported from the Woxi project (mod difference_quotient). The bare {x} form (no step) + // stays unevaluated; only the {x, h} step form evaluates. + check("DifferenceQuotient(x, x)", // + "DifferenceQuotient(x,x)"); + check("DifferenceQuotient(x^2, x)", // + "DifferenceQuotient(x^2,x)"); + check("DifferenceQuotient(x^3, {x, h})", // + "h^2+3*h*x+3*x^2"); + check("DifferenceQuotient(5, {x, h})", // + "0"); } @Test @@ -5658,6 +5668,40 @@ public void testDiscreteRatio() { "4"); check("DiscreteRatio(x^n * y^m, n, m)", // "1"); + + // ported from the Woxi project (mod discrete_ratio). The ratio is reduced to a single + // grouped fraction; (1+2*n+n^2)/n^2 == Wolfram's (1+n)^2/n^2. + check("DiscreteRatio(n^2, n)", // + "(1+2*n+n^2)/n^2"); + // common polynomial factors cancel + check("DiscreteRatio(n^2 + n, n)", // + "(2+n)/n"); + check("DiscreteRatio(f(n), n)", // + "f(1+n)/f(n)"); + check("DiscreteRatio(a^n, n)", // + "a"); + check("DiscreteRatio(2^n*3^n, n)", // + "6"); + // higher order applies the operator repeatedly + check("DiscreteRatio(f(n), {n, 3})", // + "(f(1+n)^3*f(3+n))/(f(n)*f(2+n)^3)"); + check("DiscreteRatio(f(n), {n, 0})", // + "f(n)"); + // step in a three-element spec {var, order, step} + check("DiscreteRatio(f(n), {n, 1, 2})", // + "f(2+n)/f(n)"); + check("DiscreteRatio(f(n), {n, 2, 3})", // + "(f(n)*f(6+n))/f(3+n)^2"); + check("DiscreteRatio(f(n, m), n, m)", // + "(f(n,m)*f(1+n,1+m))/(f(n,1+m)*f(1+n,m))"); + // factorial/Gamma/Pochhammer/Binomial ratios reduce + check("DiscreteRatio(Pochhammer(n, 3), n)", // + "(3+n)/n"); + check("DiscreteRatio(Binomial(n, 2), n)", // + "(1+n)/(-1+n)"); + // a symbolic order stays unevaluated + check("DiscreteRatio(f(n), {n, h})", // + "DiscreteRatio(f(n),{n,h})"); } @Test @@ -5687,6 +5731,29 @@ public void testDiscreteShift() { "f(h+x)"); check("DiscreteShift(f(x), {x,2,h})", // "f(2*h+x)"); + + // ported from the Woxi project (mod discrete_shift). A top-level Plus result is expanded; a + // single power/product is kept; an integer shift combines rational summands over a common + // denominator; a symbolic shift stays unfolded. + check("DiscreteShift(n^2 + 3*n + 1, n)", // + "5+5*n+n^2"); + check("DiscreteShift(a*n^2 + b*n, n)", // + "a+b+2*a*n+b*n+a*n^2"); + check("DiscreteShift(2*n^2, n)", // + "2*(1+n)^2"); + check("DiscreteShift({n, n^2}, n)", // + "{1+n,(1+n)^2}"); + check("DiscreteShift(1/(2*n + 1), n)", // + "1/(3+2*n)"); + check("DiscreteShift(1/(3*n - 2), n)", // + "1/(1+3*n)"); + check("DiscreteShift(1/(2*n + 1), {n, 2})", // + "1/(5+2*n)"); + check("DiscreteShift(1/(2*n + 1) + n, n)", // + "(4+5*n+2*n^2)/(3+2*n)"); + // a symbolic shift is left unfolded + check("DiscreteShift(1/(2*n + 1), {n, k})", // + "1/(1+2*(k+n))"); } @Test diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java index 55d888fd9c..b6891ca522 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java @@ -1810,4 +1810,153 @@ public void testHolonomicShiftedExpansion() { check("SeriesCoefficient(1/(2-x), {x, 0, n})", // "Piecewise({{(1/2)^(1+n),n>=0}},0)"); } + + // ========================================================== + // PadeApproximant: functions whose Series is flagged Puiseux (because of a Sqrt) + // but which still expand into integer powers of x. Previously these returned + // unevaluated; now the approximant is computed from the normal Taylor polynomial. + // (gap found via the Woxi test suite) + // ========================================================== + @Test + public void testPadeApproximantPuiseuxNormal() { + check("PadeApproximant(Sqrt(1+x), {x, 0, {1, 1}})", // + "(1+3/4*x)/(1+x/4)"); + check("PadeApproximant(Sqrt(1+x), {x, 0, {2, 2}})", // + "(1+5/4*x+5/16*x^2)/(1+3/4*x+x^2/16)"); + // single-order spec is the diagonal approximant + check("PadeApproximant(Sqrt(1+x), {x, 0, 1})", // + "(1+3/4*x)/(1+x/4)"); + // Cos(Sqrt(x)) is likewise integer-power in x + check("PadeApproximant(Cos(Sqrt(x)), {x, 0, {1, 1}})", // + "(1-5/12*x)/(1+x/12)"); + } + + // ========================================================== + // Ported from the Woxi project (tests/interpreter_tests/calculus.rs). Symja's + // canonical output form sometimes differs cosmetically but is mathematically equal, + // e.g. -I*1/2 == -I/2, 1/n! == n!^(-1), (-1)^(1+n)/n == -((-1)^n/n). + // Series results print in Symja's normal series form rather than SeriesData[...]. + // ========================================================== + @Test + public void testResidueWoxi() { + check("Residue(1/z, {z, 0})", "1"); + check("Residue(1/(z - 2), {z, 2})", "1"); + check("Residue(1/(z^2 - 1), {z, 1})", "1/2"); + check("Residue(Sin(z)/z^2, {z, 0})", "1"); + check("Residue(1/z^2, {z, 0})", "0"); + check("Residue(Exp(z)/z^3, {z, 0})", "1/2"); + check("Residue((z + 1)/(z - 1)^2, {z, 1})", "1"); + check("Residue(z^2, {z, 0})", "0"); + check("Residue(Cot(z), {z, 0})", "1"); + // -I*1/2 == -I/2 + check("Residue(1/(z^2 + 1), {z, I})", "-I*1/2"); + check("Residue(1/(z - a), {z, a})", "1"); + check("Residue(b/(z - a)^2, {z, a})", "0"); + check("Residue(1/(z^4 - 1), {z, 1})", "1/4"); + check("Residue(z/(z^2 + 1), {z, I})", "1/2"); + check("Residue(Gamma(z), {z, 0})", "1"); + check("Residue(Gamma(z), {z, -1})", "-1"); + check("Residue(Zeta(z), {z, 1})", "1"); + check("Residue(1/Sin(z), {z, Pi})", "-1"); + check("Residue(Tan(z), {z, Pi/2})", "-1"); + check("Residue(Sec(z), {z, Pi/2})", "-1"); + check("Residue(1/Sin(z)^2, {z, 0})", "0"); + check("Residue(f(z)/z, {z, 0})", "f(0)"); + } + + @Test + public void testInverseSeriesWoxi() { + // reversion of Sin gives the ArcSin series + check("InverseSeries(Series(Sin(x), {x, 0, 5}))", // + "x+x^3/6+3/40*x^5+O(x)^6"); + // reversion of Exp[x]-1 gives Log[1+x] + check("InverseSeries(Series(Exp(x) - 1, {x, 0, 4}))", // + "x-x^2/2+x^3/3-x^4/4+O(x)^5"); + // reversion of Tan gives ArcTan + check("InverseSeries(Series(Tan(x), {x, 0, 6}))", // + "x-x^3/3+x^5/5+O(x)^7"); + // polynomial reversion (Catalan numbers) + check("InverseSeries(Series(x + x^2, {x, 0, 5}))", // + "x-x^2+2*x^3-5*x^4+14*x^5+O(x)^6"); + check("InverseSeries(Series(2*x + 3*x^2, {x, 0, 4}))", // + "x/2-3/8*x^2+9/16*x^3-135/128*x^4+O(x)^5"); + // two-argument form renames the variable + check("InverseSeries(Series(Sin(x), {x, 0, 5}), y)", // + "y+y^3/6+3/40*y^5+O(y)^6"); + check("Normal(InverseSeries(Series(Sin(x), {x, 0, 5})))", // + "x+x^3/6+3/40*x^5"); + } + + @Test + public void testComposeSeriesWoxi() { + check("ComposeSeries(Series(Exp(x), {x, 0, 3}), Series(Sin(x), {x, 0, 3}))", // + "1+x+x^2/2+O(x)^4"); + check("ComposeSeries(Series(Exp(y), {y, 0, 4}), Series(x + x^2, {x, 0, 4}))", // + "1+x+3/2*x^2+7/6*x^3+25/24*x^4+O(x)^5"); + // inner nmin = 2 caps the result order + check("ComposeSeries(Series(1/(1-y), {y, 0, 2}), Series(x^2, {x, 0, 5}))", // + "1+x^2+O(x)^4"); + // truncation capped by inner accuracy + check("ComposeSeries(Series(1/(1-y), {y, 0, 5}), Series(Sin(x), {x, 0, 2}))", // + "1+x+x^2+O(x)^3"); + check("ComposeSeries(Series(Log(1+y), {y, 0, 4}), Series(x + x^3, {x, 0, 4}))", // + "x-x^2/2+4/3*x^3-5/4*x^4+O(x)^5"); + // symbolic outer coefficients + check("ComposeSeries(Series(f(y), {y, 0, 3}), Series(x + x^2, {x, 0, 3}))", // + "f(0)+f'(0)*x+(f'(0)+f''(0)/2)*x^2+(f''(0)+Derivative(3)[f][0]/6)*x^3+O(x)^4"); + } + + @Test + public void testPadeApproximantWoxi() { + check("PadeApproximant(Exp(x), {x, 0, {2, 2}})", // + "(1+x/2+x^2/12)/(1-x/2+x^2/12)"); + check("PadeApproximant(Exp(x), {x, 0, {1, 1}})", // + "(1+x/2)/(1-x/2)"); + // single-order spec is the diagonal approximant + check("PadeApproximant(Exp(x), {x, 0, 2})", // + "(1+x/2+x^2/12)/(1-x/2+x^2/12)"); + check("PadeApproximant(Sin(x), {x, 0, 4})", // + "(x-31/294*x^3)/(1+3/49*x^2+11/5880*x^4)"); + check("PadeApproximant(Exp(x), {x, 0, {3, 2}})", // + "(1+3/5*x+3/20*x^2+x^3/60)/(1-2/5*x+x^2/20)"); + check("PadeApproximant(Cos(x), {x, 0, {2, 2}})", // + "(1-5/12*x^2)/(1+x^2/12)"); + check("PadeApproximant(Log(1 + x), {x, 0, {2, 2}})", // + "(x+x^2/2)/(1+x+x^2/6)"); + check("PadeApproximant(ArcTan(x), {x, 0, {3, 2}})", // + "(x+4/15*x^3)/(1+3/5*x^2)"); + // zero denominator degree is the Taylor polynomial + check("PadeApproximant(Exp(x), {x, 0, {2, 0}})", // + "1+x+x^2/2"); + } + + @Test + public void testSeriesCoefficientWoxi() { + check("SeriesCoefficient(1/(1-x), {x, 0, 5})", "1"); + check("SeriesCoefficient(Exp(x), {x, 0, 3})", "1/6"); + check("SeriesCoefficient(Sin(x), {x, 0, 5})", "1/120"); + check("SeriesCoefficient(Log(1+x), {x, 0, 4})", "-1/4"); + // Sin has no even-order terms + check("SeriesCoefficient(Sin(x), {x, 0, 4})", "0"); + check("SeriesCoefficient(2*x, {x, 0, 2})", "0"); + check("SeriesCoefficient(Exp(Sin(x)), {x, 0, 4})", "-1/8"); + // query a computed SeriesData with the {x, x0, n} spec + check("SeriesCoefficient(Series(Exp(x), {x, 0, 10}), {x, 0, 5})", "1/120"); + check("SeriesCoefficient(Series(Sin(x), {x, 0, 10}), {x, 0, 3})", "-1/6"); + // symbolic index returns the general term as a Piecewise; 1/n! == n!^(-1) + check("SeriesCoefficient(Exp(x), {x, 0, n})", // + "Piecewise({{1/n!,n>=0}},0)"); + check("SeriesCoefficient(1/(1 - x), {x, 0, n})", // + "Piecewise({{1,n>=0}},0)"); + check("SeriesCoefficient(1/(1 - x)^2, {x, 0, n})", // + "Piecewise({{1+n,n>=0}},0)"); + check("SeriesCoefficient(Cosh(x), {x, 0, n})", // + "Piecewise({{1/n!,Mod(n,2)==0&&n>=0}},0)"); + // (-1)^(1+n)/n == -((-1)^n/n) + check("SeriesCoefficient(Log(1 + x), {x, 0, n})", // + "Piecewise({{(-1)^(1+n)/n,n>=1}},0)"); + // an unrecognized symbolic-index form stays unevaluated + check("SeriesCoefficient(Exp(x^2), {x, 0, n})", // + "SeriesCoefficient(E^x^2,{x,0,n})"); + } } From 8afef5db8b975dbba656ae1a9bc093e3409420b0 Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 16:49:04 +0200 Subject: [PATCH 3/8] Defer RationalIntegration RootSum results until after Rubi RationalIntegration.integrate() ran before the Rubi rules in the Integrate cascade. That is form-safe for its polynomial, Horowitz-Ostrogradsky and degree-1/2 log/arctan parts, but its RootSum branch (irreducible denominator factors of degree >= 5) pre-empted Rubi's frequently simpler closed form and produced RootSum antiderivatives the downstream check often cannot verify. Add a RootSumMode {DEFER, EMIT} parameter to RationalIntegration.integrate(): - Pre-Rubi, Integrate now calls it with DEFER. When the antiderivative is essentially a bare RootSum, DEFER returns F.NIL so Rubi gets first crack (like the existing degree-3/4 deferral). A mixed Log(..)+RootSum(..) is a Plus and is still emitted pre-Rubi, preserving form-safety and avoiding a DerivativeDivides/cancelGCD ClassCastException seen when such integrands are routed into Rubi. - Post-Rubi, a new first stage calls it with EMIT, so a pure degree->=5 rational denominator that Rubi leaves unevaluated re-emits the identical RootSum. Output is unchanged except where Rubi has a simpler closed form. The 3-arg integrate() wrapper delegates with EMIT, so existing callers (RadicalSubstitution, IntegrateAlgorithmsTest) are unaffected. Full matheclipse-core suite green (3509 tests, 0 failures/errors). Co-Authored-By: Claude Opus 4.8 --- .../core/integrate/RationalIntegration.java | 55 ++++++++++++++++++- .../core/reflection/system/Integrate.java | 23 +++++--- 2 files changed, 70 insertions(+), 8 deletions(-) diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/integrate/RationalIntegration.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/integrate/RationalIntegration.java index d9c31fd247..37ae2f6979 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/integrate/RationalIntegration.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/integrate/RationalIntegration.java @@ -53,7 +53,22 @@ public class RationalIntegration { private RationalIntegration() {} /** - * Integrate a rational function in x. + * Controls whether the logarithmic part may emit a {@link F#RootSum} antiderivative for an + * irreducible denominator factor of degree >= 5. + */ + public enum RootSumMode { + /** + * Irreducible factors of degree >= 5 return {@link F#NIL} (defer to the Rubi rules), exactly + * like the degree 3/4 factors already do. + */ + DEFER, + /** Irreducible factors of degree >= 5 produce a {@link F#RootSum} antiderivative. */ + EMIT + } + + /** + * Integrate a rational function in x, emitting a {@link F#RootSum} antiderivative for + * irreducible denominator factors of degree >= 5 (i.e. {@link RootSumMode#EMIT}). * * @param integrand the integrand * @param x the integration variable @@ -61,6 +76,44 @@ private RationalIntegration() {} * @return the antiderivative or {@link F#NIL} */ public static IExpr integrate(IExpr integrand, IExpr x, EvalEngine engine) { + return integrate(integrand, x, engine, RootSumMode.EMIT); + } + + /** + * Integrate a rational function in x. + * + *

+ * In {@link RootSumMode#DEFER} mode, when the computed antiderivative is essentially a bare + * {@link F#RootSum} (a single top-level term containing a {@code RootSum}, with no additional + * closed-form {@code Log}/{@code ArcTan}/rational/polynomial terms), this returns {@link F#NIL} so + * the caller can hand the integrand to the Rubi rules first — those often have a far simpler + * closed form. A mixed antiderivative (e.g. {@code Log(x-1) + RootSum[...]}) is a {@code Plus} and + * is not deferred: it is a correct-by-construction closed form and deferring it would + * route the integrand into a fragile Rubi recursion. The RootSum is re-emitted by a post-Rubi + * fallback that calls this method with {@link RootSumMode#EMIT}. + * + * @param integrand the integrand + * @param x the integration variable + * @param engine the evaluation engine + * @param mode whether an antiderivative that is essentially a bare {@link F#RootSum} is emitted + * ({@link RootSumMode#EMIT}) or deferred to the Rubi rules ({@link RootSumMode#DEFER}) + * @return the antiderivative or {@link F#NIL} + */ + public static IExpr integrate(IExpr integrand, IExpr x, EvalEngine engine, RootSumMode mode) { + IExpr result = integrateRationalFunction(integrand, x, engine); + if (mode == RootSumMode.DEFER && result.isPresent() && !result.isPlus() + && !result.isFree(F.RootSum)) { + // Antiderivative is essentially a bare RootSum: defer to the Rubi rules. + return F.NIL; + } + return result; + } + + /** + * Compute the antiderivative of a rational function, always emitting a {@link F#RootSum} for + * irreducible denominator factors of degree >= 5. + */ + private static IExpr integrateRationalFunction(IExpr integrand, IExpr x, EvalEngine engine) { if (!Config.INTEGRATE_ALGORITHM_RATIONAL) { return F.NIL; } diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Integrate.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Integrate.java index 0cb8ac1f92..a81604d1ed 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Integrate.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/Integrate.java @@ -509,8 +509,13 @@ public IExpr evaluate(IAST holdallAST, final int argSize, final IExpr[] option, // un-wired (still reachable via the Method -> option and their direct tests). // Stage: native rational function integration - // (Hermite/Horowitz-Ostrogradsky reduction + Lazard-Rioboo-Trager logarithmic part) - result = RationalIntegration.integrate(fx, x, engine); + // (Hermite/Horowitz-Ostrogradsky reduction + Lazard-Rioboo-Trager logarithmic part). + // RootSumMode.DEFER: when the antiderivative is essentially a bare RootSum, defer the + // integrand to the Rubi rules (which often have a far simpler closed form) and re-emit + // the RootSum only as a post-Rubi fallback (see below). Closed-form results, including a + // mixed Log(..)+RootSum(..), are still produced here. + result = + RationalIntegration.integrate(fx, x, engine, RationalIntegration.RootSumMode.DEFER); if (result.isPresent()) { return result; } @@ -551,6 +556,15 @@ public IExpr evaluate(IAST holdallAST, final int argSize, final IExpr[] option, // unevaluated. Each self-verifies (D(result) == integrand). The deterministic, form-safe // stages (rational, radical, Chebyshev) run before Rubi (above). if (Config.INTEGRATE_ALGORITHMS) { + // RootSum fallback for rational functions whose denominator has an irreducible factor of + // degree >= 5, deferred from the pre-Rubi rational stage (RootSumMode.DEFER above). Runs + // only now that Rubi left the integral unevaluated, so Rubi's simpler closed form (when it + // has one) always wins. Correct-by-construction (Trager), reuses the full general logic. + result = + RationalIntegration.integrate(fx, x, engine, RationalIntegration.RootSumMode.EMIT); + if (result.isPresent()) { + return result; + } // Weierstrass t=Tan(x/2) substitution for rational trigonometric integrands. result = WeierstrassIntegration.integrate(fx, x, engine); if (result.isPresent()) { @@ -573,11 +587,6 @@ public IExpr evaluate(IAST holdallAST, final int argSize, final IExpr[] option, } } - // // RootSum fallback for 1/p(x), irreducible degree >= 5 --- - // IExpr rootSumResult = integrateOneOverPoly(fx, x, engine); - // if (rootSumResult.isPresent()) { - // return rootSumResult; - // } } return evaled ? ast : F.NIL; } finally { From 4d8e3022fbc2375140c992cff8eee5c73e5c83c3 Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 17:00:30 +0200 Subject: [PATCH 4/8] Improve AST1#equals() method --- .../org/matheclipse/core/expression/AST1.java | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/AST1.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/AST1.java index 2e9ed5ab49..a801f8c135 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/AST1.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/expression/AST1.java @@ -165,16 +165,18 @@ public boolean equals(final Object obj) { return false; } final IAST list = (IAST) rhs; - if (lhs.arg0 != list.head() && lhs.arg0 instanceof ISymbol) { - // compared with ISymbol object identity - return false; - } if (list.size() != SIZE) { return false; } - if (lhs.arg0 != list.head() && !(lhs.arg0 instanceof ISymbol) - && !lhs.arg0.equals(list.head())) { - return false; + if (lhs.arg0 instanceof ISymbol) { + if (lhs.arg0 != list.head()) { + // compared with ISymbol object identity + return false; + } + } else { + if (!lhs.arg0.equals(list.head())) { + return false; + } } final IExpr lhsArg1 = lhs.arg1; if (lhsArg1.getClass() != AST1.class) { From 8bc9074f3a98f095ece37c6ab2c47d7db8b61c52 Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 17:04:47 +0200 Subject: [PATCH 5/8] More `RootSum` JUnit tests --- .../core/reflection/system/RootSum.java | 5 + .../core/reflection/system/DTest.java | 6 + .../core/reflection/system/RootSumTest.java | 106 ++++++++++++++++++ 3 files changed, 117 insertions(+) create mode 100644 symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/RootSumTest.java diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/RootSum.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/RootSum.java index ab6e7d1f55..4fc3d1ca0c 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/RootSum.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/RootSum.java @@ -164,6 +164,11 @@ private static IExpr traceOverRoots(IExpr h, IExpr pMonic, ISymbol r, IExpr[] po * pMonic are not coprime (i.e. a root of the denominator is a root of the polynomial). */ private static IExpr modularInverse(IExpr b, IExpr pMonic, ISymbol r, EvalEngine engine) { + if (b.isFree(r)) { + // b is a unit (constant) with respect to r, so its inverse modulo pMonic is just 1/b. + // PolynomialExtendedGCD would return the degenerate cofactor s == 0 in this case. + return b.isZero() ? F.NIL : engine.evaluate(F.Power(b, F.CN1)); + } IExpr extendedGCD = S.PolynomialExtendedGCD.of(engine, b, pMonic, r); if (!extendedGCD.isList() || extendedGCD.size() != 3) { return F.NIL; diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/DTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/DTest.java index c409f8103c..222a083270 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/DTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/DTest.java @@ -807,6 +807,11 @@ public void testD002() { @Test public void testDRootSum() { + // #1^7/(8*#1^7) cancels to 1/8, so the summand is Log(x-#1)/8; its derivative telescopes + // over the roots of 1+#1^8 to the integrand x^7/(1+x^8). + check("D(RootSum(1+#1^8&,(Log(x-#1)*#1^7)/(8*#1^7)&),x)", // + "x^7/(1+x^8)"); + // A definite integral whose antiderivative is a RootSum(...Log(x-#1)...) is evaluated at the // bounds by substituting x->1 and x->3 into the summand (Limit of a RootSum), matching // -RootSum(...,Log(1-#1)/...&) + RootSum(...,Log(3-#1)/...&). @@ -819,5 +824,6 @@ public void testDRootSum() { check( "D(RootSum(#1^6-5*#1^4+5*#1^2+4&,((#1^4-3*#1^2+6)*Log(x-#1))/(6*#1^5-20*#1^3+10*#1)&) ,x)", // "(6-3*x^2+x^4)/(4+5*x^2-5*x^4+x^6)"); + } } diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/RootSumTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/RootSumTest.java new file mode 100644 index 0000000000..de6e7b4c46 --- /dev/null +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/reflection/system/RootSumTest.java @@ -0,0 +1,106 @@ +package org.matheclipse.core.reflection.system; + +import org.junit.jupiter.api.Test; +import org.matheclipse.core.system.ExprEvaluatorTestCase; + +/** + * Tests for {@link RootSum}. + * + */ +public class RootSumTest extends ExprEvaluatorTestCase { + + @Test + public void testSumOfRoots() { + // The sum of the roots of x^n-c has no x^(n-1) term, so it is 0. + check("RootSum(#^2 - 2 &, # &)", // + "0"); + check("RootSum(#^3 - 2 &, # &)", // + "0"); + check("RootSum(#^4 - 1 &, # &)", // + "0"); + // Vieta: the roots 2 and 3 sum to 5. + check("RootSum(#^2 - 5*# + 6 &, # &)", // + "5"); + } + + @Test + public void testPowerSumsOfRoots() { + // roots 1, 2: sum of squares = 5 + check("RootSum(#^2 - 3*# + 2 &, #^2 &)", // + "5"); + // roots 2, 3: sum of cubes = 8+27 = 35 + check("RootSum(#^2 - 5*# + 6 &, #^3 &)", // + "35"); + // roots +-Sqrt(2): sum of squares = 4 + check("RootSum(#^2 - 2 &, #^2 &)", // + "4"); + // roots +-I: sum of squares = -2 + check("RootSum(#^2 + 1 &, #^2 &)", // + "-2"); + // roots 1, -1, I, -I: sum of squares = 0, sum of 4th powers = 4 + check("RootSum(#^4 - 1 &, #^2 &)", // + "0"); + check("RootSum(#^4 - 1 &, #^4 &)", // + "4"); + } + + @Test + public void testUnsolvablePolynomials() { + // The power-sum method needs no explicit roots, so it also works for polynomials + // that have no closed-form (radical) roots. + check("RootSum(#^3 - # - 1 &, #^2 &)", // + "2"); + check("RootSum(#^3 + # + 1 &, #^2 &)", // + "-2"); + check("RootSum(#^5 - # - 1 &, #^2 &)", // + "0"); + } + + @Test + public void testNonMonicAndAffineForm() { + // non-monic f: 2*x^2-8 has roots +-2; sum of squares = 8 + check("RootSum(2*#^2 - 8 &, #^2 &)", // + "8"); + // affine form 3*#^2+1: 3*(sum of squares)+1*(root count) = 3*4+2 = 14 + check("RootSum(#^2 - 2 &, 3*#^2 + 1 &)", // + "14"); + // constant form sums to constant*(number of roots) + check("RootSum(#^2 - 2 &, 5 &)", // + "10"); + // linear polynomial: single root 3, squared = 9 + check("RootSum(# - 3 &, #^2 &)", // + "9"); + } + + @Test + public void testFormDegreeExceedingPolynomialDegree() { + // A form of degree >= deg(f) is first reduced modulo f. For x^3-x-1: + // #^3 -> #+1 gives sum = 3, and #^4 -> #^2+# gives sum = 2. + check("RootSum(#^3 - # - 1 &, #^3 &)", // + "3"); + check("RootSum(#^3 - # - 1 &, #^4 &)", // + "2"); + check("RootSum(#^2 - 2 &, #^4 &)", // + "8"); + } + + @Test + public void testInertForms() { + // A non-polynomial (Log) form is left unevaluated; this is the shape produced by + // Integrate for an irreducible high-degree denominator. + check("RootSum(1 + #1 + #1^2 + #1^3 + #1^4 & , Log(x + #1) & )", // + "RootSum(1+#1+#1^2+#1^3+#1^4&,Log(x+#1)&)"); + } + + @Test + public void testDerivativeIsRationalFunction() { + // D(RootSum(p&, Log(x-#1)/p'(#1) &), x) telescopes back to the integrand 1/p(x). + check("D(RootSum(#1^3 - 2 &, Log(x - #1)/(3*#1^2) &), x)", // + "1/(-2+x^3)"); + check("D(RootSum(#1^5 + #1 - 7 &, Log(x - #1)/(5*#1^4 + 1) &), x)", // + "1/(-7+x+x^5)"); + // The Integrate-produced antiderivative of x^7/(1+x^8) differentiates back to it. + check("D(RootSum(1+#1^8&,(Log(x-#1)*#1^7)/(8*#1^7)&),x)", // + "x^7/(1+x^8)"); + } +} From c3add885334e573ea873a9c6c0d6e4a444006068 Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 17:07:16 +0200 Subject: [PATCH 6/8] DifferenceDelta,DiscreteRatio improvements --- .../matheclipse/core/reflection/system/DiscreteRatio.java | 2 +- .../org/matheclipse/core/system/DifferenceDeltaTest.java | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java index 0b83f8801d..80cba22849 100644 --- a/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java +++ b/symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/reflection/system/DiscreteRatio.java @@ -64,7 +64,7 @@ public IExpr evaluate(final IAST ast, EvalEngine engine) { // Together cancels common polynomial factors (e.g. DiscreteRatio(n^2 + n, n) -> (2 + n)/n) // while keeping a single grouped fraction (Cancel would rewrite (1 + n)^2/n^2 as (1 + 1/n)^2). // Only adopt the FunctionExpand result when it fully reduces (no Gamma/Factorial residue left, - // e.g. a scaled argument like (2 n)! that woxi/Symja cannot reduce further). + // e.g. a scaled argument like (2 n)! that Symja cannot reduce further). IExpr expanded = engine.evaluate(F.Together(F.FunctionExpand(result))); if (expanded.isFree(S.Gamma, true) && expanded.isFree(S.Factorial, true)) { return expanded; diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java index fdc3839de0..4c900686d3 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/DifferenceDeltaTest.java @@ -5,7 +5,7 @@ public class DifferenceDeltaTest extends ExprEvaluatorTestCase { @Test - public void testDifferenceDelta() { + public void testDifferenceDelta001() { check("DifferenceDelta({f(i), g(i)}, i)", // "{-f(i)+f(1+i),-g(i)+g(1+i)}"); check("DifferenceDelta(Cosh(a*i+b),{i,2,h})", // @@ -73,14 +73,13 @@ public void testDifferenceDeltaIndefiniteSum() { } // ========================================================== - // Ported from the Woxi project (tests/interpreter_tests/calculus.rs, mod difference_delta). // The raw forward difference is now brought into canonical form: a numeric step is factored // (x^2 + x -> 2*(1 + x), and rational summands combine over a common denominator), a symbolic // step is expanded. Symja's form is sometimes cosmetically different but mathematically equal, // e.g. 1/((-1-n)*n) == -1/(n*(1+n)). // ========================================================== @Test - public void testDifferenceDeltaWoxi() { + public void testDifferenceDelta002() { check("DifferenceDelta(5, x)", "0"); check("DifferenceDelta(x, x)", "1"); check("DifferenceDelta(a*x + b, x)", "a"); From 1bce3ce00ecce144053f8fe5048c30a51a9bc196 Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 17:08:43 +0200 Subject: [PATCH 7/8] More JUnit test cases --- .../core/system/LaplaceTransformTest.java | 28 ++-- .../core/system/LowercaseTestCase.java | 11 +- .../matheclipse/core/system/SeriesTest.java | 25 ++-- .../core/rubi/AbstractRubiTestCase.java | 133 ++++++++++++++++-- .../core/rubi/InverseHyperbolicFunctions.java | 5 +- 5 files changed, 148 insertions(+), 54 deletions(-) diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java index 0ac6469840..0a32d10f11 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LaplaceTransformTest.java @@ -324,14 +324,8 @@ public void testLaplaceTransformEdgeCases() { "(a*b)/s"); } - // ========================================================== - // Ported from the Woxi project (laplace_transform module). Symja's canonical output - // form sometimes differs cosmetically from wolframscript / Woxi but is mathematically - // equal, e.g. 1/s == s^(-1), 1/(1+s^2) == (1 + s^2)^(-1), E^(-a*t) == E^(-(a*t)). - // ========================================================== - @Test - public void testLaplaceTransformWoxiMonomials() { + public void testLaplaceTransformMonomials() { // 1/s == s^(-1) check("LaplaceTransform(1, t, s)", // "1/s"); @@ -347,7 +341,7 @@ public void testLaplaceTransformWoxiMonomials() { } @Test - public void testLaplaceTransformWoxiTrig() { + public void testLaplaceTransformTrig() { // 1/(1+s^2) == (1 + s^2)^(-1) check("LaplaceTransform(Sin(t), t, s)", // "1/(1+s^2)"); @@ -358,7 +352,7 @@ public void testLaplaceTransformWoxiTrig() { } @Test - public void testLaplaceTransformWoxiExponential() { + public void testLaplaceTransformExponential() { // 1/(a+s) == (a + s)^(-1) check("LaplaceTransform(Exp(-a*t), t, s)", // "1/(a+s)"); @@ -368,24 +362,20 @@ public void testLaplaceTransformWoxiExponential() { } @Test - public void testLaplaceTransformWoxiLinearity() { + public void testLaplaceTransformLinearity() { check("LaplaceTransform(3*t^2 + 2*Sin(t), t, s)", // "6/s^3+2/(1+s^2)"); } @Test - public void testLaplaceTransformWoxiUnevaluated() { + public void testLaplaceTransformUnevaluated() { // Unknown functions should return unevaluated check("LaplaceTransform(BesselJ(0, t), t, s)", // "LaplaceTransform(BesselJ(0,t),t,s)"); } - // ========================================================== - // InverseLaplaceTransform, ported from the Woxi project. - // ========================================================== - @Test - public void testInverseLaplaceTransformWoxiMonomials() { + public void testInverseLaplaceTransformMonomials() { check("InverseLaplaceTransform(1/s, s, t)", // "1"); check("InverseLaplaceTransform(1/s^2, s, t)", // @@ -397,7 +387,7 @@ public void testInverseLaplaceTransformWoxiMonomials() { } @Test - public void testInverseLaplaceTransformWoxiTrig() { + public void testInverseLaplaceTransformTrig() { check("InverseLaplaceTransform(1/(s^2 + 1), s, t)", // "Sin(t)"); check("InverseLaplaceTransform(s/(s^2 + 1), s, t)", // @@ -409,7 +399,7 @@ public void testInverseLaplaceTransformWoxiTrig() { } @Test - public void testInverseLaplaceTransformWoxiExponential() { + public void testInverseLaplaceTransformExponential() { check("InverseLaplaceTransform(1/(s - a), s, t)", // "E^(a*t)"); // E^(-a*t) == E^(-(a*t)) @@ -418,7 +408,7 @@ public void testInverseLaplaceTransformWoxiExponential() { } @Test - public void testInverseLaplaceTransformWoxiUnevaluated() { + public void testInverseLaplaceTransformUnevaluated() { check("InverseLaplaceTransform(Log(s), s, t)", // "InverseLaplaceTransform(Log(s),s,t)"); } diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java index 5e8c99237c..b429ae11b5 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/LowercaseTestCase.java @@ -5517,8 +5517,7 @@ public void testDifferenceQuotient() { check("dq = DifferenceQuotient(f(x), {x, h})", // "(-x*Sin(x)+(h+x)*Sin(h+x))/h"); - // ported from the Woxi project (mod difference_quotient). The bare {x} form (no step) - // stays unevaluated; only the {x, h} step form evaluates. + // The bare {x} form (no step) stays unevaluated; only the {x, h} step form evaluates. check("DifferenceQuotient(x, x)", // "DifferenceQuotient(x,x)"); check("DifferenceQuotient(x^2, x)", // @@ -5669,8 +5668,7 @@ public void testDiscreteRatio() { check("DiscreteRatio(x^n * y^m, n, m)", // "1"); - // ported from the Woxi project (mod discrete_ratio). The ratio is reduced to a single - // grouped fraction; (1+2*n+n^2)/n^2 == Wolfram's (1+n)^2/n^2. + // The ratio is reduced to a single grouped fraction; (1+2*n+n^2)/n^2 == Wolfram's (1+n)^2/n^2. check("DiscreteRatio(n^2, n)", // "(1+2*n+n^2)/n^2"); // common polynomial factors cancel @@ -5732,9 +5730,8 @@ public void testDiscreteShift() { check("DiscreteShift(f(x), {x,2,h})", // "f(2*h+x)"); - // ported from the Woxi project (mod discrete_shift). A top-level Plus result is expanded; a - // single power/product is kept; an integer shift combines rational summands over a common - // denominator; a symbolic shift stays unfolded. + // A top-level Plus result is expanded; a single power/product is kept; an integer shift + // combines rational summands over a common denominator; a symbolic shift stays unfolded. check("DiscreteShift(n^2 + 3*n + 1, n)", // "5+5*n+n^2"); check("DiscreteShift(a*n^2 + b*n, n)", // diff --git a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java index b6891ca522..a39a0c7cb8 100644 --- a/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java +++ b/symja_android_library/matheclipse-core/src/test/java/org/matheclipse/core/system/SeriesTest.java @@ -62,7 +62,7 @@ public void testComposeSeries() { } @Test - public void testInverseSeries() { + public void testInverseSeries001() { check("InverseSeries(Sin(x), x)", // "x+x^3/6+3/40*x^5+O(x)^6"); @@ -655,9 +655,9 @@ public void testSeriesCoefficient() { } @Test - public void testSeriesCoefficientBinomialWoxi() { + public void testSeriesCoefficientBinomial() { // Finite binomial powers (a+b*x)^p with a symbolic index n now return a closed Binomial - // Piecewise instead of an opaque DifferenceRoot (matches Mathematica / WOXI). + // Piecewise instead of an opaque DifferenceRoot . check("SeriesCoefficient((1+x)^5, {x, 0, n})", // "Piecewise({{Binomial(5,n),0<=n<=5}},0)"); check("SeriesCoefficient((1+2*x)^3, {x, 0, n})", // @@ -1118,7 +1118,7 @@ public void testSeriesCoefficientHypergeometricPFQ() { } @Test - public void testResidue() { + public void testResidue001() { // Partial fraction - fastpath) check("Residue(1 / (x - 2), {x, 2})", // "1"); @@ -1815,7 +1815,6 @@ public void testHolonomicShiftedExpansion() { // PadeApproximant: functions whose Series is flagged Puiseux (because of a Sqrt) // but which still expand into integer powers of x. Previously these returned // unevaluated; now the approximant is computed from the normal Taylor polynomial. - // (gap found via the Woxi test suite) // ========================================================== @Test public void testPadeApproximantPuiseuxNormal() { @@ -1831,14 +1830,8 @@ public void testPadeApproximantPuiseuxNormal() { "(1-5/12*x)/(1+x/12)"); } - // ========================================================== - // Ported from the Woxi project (tests/interpreter_tests/calculus.rs). Symja's - // canonical output form sometimes differs cosmetically but is mathematically equal, - // e.g. -I*1/2 == -I/2, 1/n! == n!^(-1), (-1)^(1+n)/n == -((-1)^n/n). - // Series results print in Symja's normal series form rather than SeriesData[...]. - // ========================================================== @Test - public void testResidueWoxi() { + public void testResidue002() { check("Residue(1/z, {z, 0})", "1"); check("Residue(1/(z - 2), {z, 2})", "1"); check("Residue(1/(z^2 - 1), {z, 1})", "1/2"); @@ -1865,7 +1858,7 @@ public void testResidueWoxi() { } @Test - public void testInverseSeriesWoxi() { + public void testInverseSeries002() { // reversion of Sin gives the ArcSin series check("InverseSeries(Series(Sin(x), {x, 0, 5}))", // "x+x^3/6+3/40*x^5+O(x)^6"); @@ -1888,7 +1881,7 @@ public void testInverseSeriesWoxi() { } @Test - public void testComposeSeriesWoxi() { + public void testComposeSeries002() { check("ComposeSeries(Series(Exp(x), {x, 0, 3}), Series(Sin(x), {x, 0, 3}))", // "1+x+x^2/2+O(x)^4"); check("ComposeSeries(Series(Exp(y), {y, 0, 4}), Series(x + x^2, {x, 0, 4}))", // @@ -1907,7 +1900,7 @@ public void testComposeSeriesWoxi() { } @Test - public void testPadeApproximantWoxi() { + public void testPadeApproximant002() { check("PadeApproximant(Exp(x), {x, 0, {2, 2}})", // "(1+x/2+x^2/12)/(1-x/2+x^2/12)"); check("PadeApproximant(Exp(x), {x, 0, {1, 1}})", // @@ -1931,7 +1924,7 @@ public void testPadeApproximantWoxi() { } @Test - public void testSeriesCoefficientWoxi() { + public void testSeriesCoefficient002() { check("SeriesCoefficient(1/(1-x), {x, 0, 5})", "1"); check("SeriesCoefficient(Exp(x), {x, 0, 3})", "1/6"); check("SeriesCoefficient(Sin(x), {x, 0, 5})", "1/120"); diff --git a/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/AbstractRubiTestCase.java b/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/AbstractRubiTestCase.java index 6305da0d1b..42e3d5db58 100644 --- a/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/AbstractRubiTestCase.java +++ b/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/AbstractRubiTestCase.java @@ -1,6 +1,8 @@ package org.matheclipse.core.rubi; import java.io.StringWriter; +import java.util.Set; +import java.util.TreeSet; import java.util.concurrent.TimeUnit; import org.matheclipse.core.basic.Config; import org.matheclipse.core.eval.EvalControlledCallable; @@ -11,9 +13,10 @@ import org.matheclipse.core.eval.exception.Validate; import org.matheclipse.core.expression.F; import org.matheclipse.core.form.output.OutputFormFactory; +import org.matheclipse.core.interfaces.IAST; +import org.matheclipse.core.interfaces.IASTAppendable; import org.matheclipse.core.interfaces.IExpr; import org.matheclipse.core.parser.ExprParser; -import org.matheclipse.gpl.numbertheory.BigIntegerPrimality; import org.matheclipse.parser.client.ParserConfig; import org.matheclipse.parser.client.SyntaxError; import org.matheclipse.parser.client.math.MathException; @@ -81,15 +84,22 @@ private String printResult(IExpr integral, IExpr result, String expectedResult, // the expressions are structurally equal return expectedResult; } else { - IExpr diff = fEvaluator.eval(F.D(result, F.symbol("x"))); - temp = fEvaluator.eval(F.Subtract(diff, integral)); - // System.out.println(temp.toString()); - expected = fEvaluator.eval(F.PossibleZeroQ(temp)); - if (expected.isTrue()) { - // the expressions are structurally equal + // Form-independent correctness check: the antiderivative must differentiate back to the + // integrand, i.e. D(result, x) - integrand == 0. This is what lets a native algorithm + // stage produce a different-looking (but correct) antiderivative than Rubi/Mathematica. + // + // The two checks are OR-ed, so their order cannot change which results are accepted - + // only what it costs to decide. The numeric check runs in milliseconds, while the + // symbolic one calls Together/Simplify, which take tens of seconds on a big difference + // and can hang outright (JAS' subresultant GCD swells on complex-rational coefficients, + // and it never checks the engine deadline, so no timeout can stop it). So try the cheap + // one first and only fall back to the symbolic route when it cannot decide. + IExpr difference = fEvaluator.eval(F.Subtract(F.D(result, F.symbol("x")), integral)); + if (isFiniteDifferenceCorrect(result, integral, F.symbol("x")) + || isZeroAntiderivativeCheck(difference)) { return expectedResult; } else { - System.out.println("PossibleZeroQ[\n" + temp.toString() + " \n]"); + System.out.println("D(result) - integrand not provably zero:\n" + difference + "\n"); } } // IExpr resultTogether= F.Together.of(F.ExpandAll(result)); @@ -113,6 +123,108 @@ private String printResult(IExpr integral, IExpr result, String expectedResult, return buf.toString(); } + /** + * Symbolic zero test for {@code D(result, x) - integrand}: {@code PossibleZeroQ} escalated + * through {@code Together} and {@code Simplify}. Works when Symja can differentiate and simplify + * cleanly. + */ + private boolean isZeroAntiderivativeCheck(IExpr difference) { + if (difference.isZero() || fEvaluator.eval(F.PossibleZeroQ(difference)).isTrue()) { + return true; + } + IExpr together = fEvaluator.eval(F.Together(difference)); + if (together.isZero() || fEvaluator.eval(F.PossibleZeroQ(together)).isTrue()) { + return true; + } + IExpr simplified = fEvaluator.eval(F.Simplify(difference)); + return simplified.isZero() || fEvaluator.eval(F.PossibleZeroQ(simplified)).isTrue(); + } + + /** + * Numeric finite-difference correctness check: {@code (result(x0+h) - result(x0-h))/(2h)} must + * equal the integrand at {@code x0}, for several {@code x0} and random values of every parameter. + * This avoids the symbolic derivative entirely, so it verifies antiderivatives (e.g. + * {@code RootSum} or special-function forms) whose symbolic {@code D} Symja cannot simplify. + */ + private boolean isFiniteDifferenceCorrect(IExpr result, IExpr integrand, IExpr x) { + // Symbols come from the integrand; the antiderivative introduces no new free symbols (e.g. + // RootSum slots are bound). Note this must not use Variables(), which reports the *polynomial* + // variables: for a trigonometric integrand it answers {a, b, Cos[c+d*x]}, which neither + // contains x nor lets Cot[c+d*x]/Sec[c+d*x] in the antiderivative become numeric. + IAST symbols = freeSymbols(integrand); + if (!symbols.exists(symbol -> symbol.equals(x))) { + return false; // nothing to differentiate against + } + final double h = 1.0e-5; + int verified = 0; + for (int trial = 0; trial < 8 && verified < 3; trial++) { + double x0 = 0.4 + 0.37 * trial; + try { + IExpr resultPlus = numericAt(result, symbols, x, x0 + h, trial); + IExpr resultMinus = numericAt(result, symbols, x, x0 - h, trial); + IExpr integrandValue = numericAt(integrand, symbols, x, x0, trial); + if (!resultPlus.isNumber() || !resultMinus.isNumber() || !integrandValue.isNumber()) { + continue; // not numerically evaluable at this point; try another + } + IExpr derivative = + fEvaluator.eval(F.Divide(F.Subtract(resultPlus, resultMinus), F.num(2.0 * h))); + IExpr error = fEvaluator.eval(F.Abs(F.Subtract(derivative, integrandValue))); + IExpr tolerance = + fEvaluator.eval(F.Times(F.num(1.0e-3), F.Plus(F.C1, F.Abs(integrandValue)))); + if (!fEvaluator.eval(F.Less(error, tolerance)).isTrue()) { + return false; + } + verified++; + } catch (RuntimeException rex) { + // try another point + } + } + return verified >= 3; + } + + /** + * The free symbols of {@code expr}, sorted, as a {@code List(...)}. A symbol is free when it + * carries neither the {@code Constant} nor the {@code NumericFunction} attribute, which excludes + * both constants such as {@code Pi} and the heads of the expression such as {@code Cos} or + * {@code Plus}. + */ + private static IAST freeSymbols(IExpr expr) { + Set symbols = new TreeSet(); + collectFreeSymbols(expr, symbols); + IASTAppendable result = F.ListAlloc(symbols.size()); + result.appendAll(symbols); + return result; + } + + private static void collectFreeSymbols(IExpr expr, Set symbols) { + if (expr.isSymbol()) { + if (expr.isVariable()) { + symbols.add(expr); + } + return; + } + if (expr.isAST()) { + IAST ast = (IAST) expr; + // start at 0 to visit the head as well + for (int i = 0; i < ast.size(); i++) { + collectFreeSymbols(ast.get(i), symbols); + } + } + } + + /** + * Numeric value of {@code expr} with {@code x = xValue} and every other free symbol randomised. + */ + private IExpr numericAt(IExpr expr, IAST symbols, IExpr x, double xValue, int trial) { + IASTAppendable rules = F.ListAlloc(symbols.size()); + for (int i = 1; i < symbols.size(); i++) { + IExpr symbol = symbols.get(i); + double value = symbol.equals(x) ? xValue : (0.7 + 0.31 * i + 0.09 * trial); + rules.append(F.Rule(symbol, F.num(value))); + } + return fEvaluator.eval(F.N(F.ReplaceAll(expr, rules))); + } + /** Evaluates the given string-expression and returns the result in OutputForm */ public String interpreter(final String inputExpression, final String expectedResult, String manuallyCheckedResult) { @@ -189,7 +301,7 @@ public void checkLength(String evalString, String expectedResult, String manuall assertEquals(expectedResult, evaledResult); } } catch (AssertionError e) { - System.out.println(getName() + " - " + evalString); + // System.out.println(getName() + " - " + evalString); throw e; } catch (Exception e) { e.printStackTrace(); @@ -213,13 +325,12 @@ public void check(String evalString, String expectedResult, String manuallyCheck checkLength(evalString, expectedResult, manuallyCheckedResult, -1); } - /** The JUnit setup method */ @Override protected void setUp() { try { super.setUp(); Config.SHORTEN_STRING_LENGTH = 80; - Config.PRIME_FACTORS = new BigIntegerPrimality(); + // Config.PRIME_FACTORS = new BigIntegerPrimality(); F.await(); // start test with fresh instance EvalEngine engine = new EvalEngine(isRelaxedSyntax); diff --git a/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/InverseHyperbolicFunctions.java b/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/InverseHyperbolicFunctions.java index 76074cbb02..b6a77d6d15 100644 --- a/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/InverseHyperbolicFunctions.java +++ b/symja_android_library/matheclipse-io/src/test/java/org/matheclipse/core/rubi/InverseHyperbolicFunctions.java @@ -262,7 +262,10 @@ public void test0128() { public void test0129() { check( // "Integrate[(a + b*ArcCosh[-1 + d*x^2])^(-3/2), x]", // - "-((Sqrt[d*x^2]*Sqrt[-2 + d*x^2])/(b*d*x*Sqrt[a + b*ArcCosh[-1 + d*x^2]])) + (Sqrt[Pi/2]*Cosh[ArcCosh[-1 + d*x^2]/2]*Erfi[Sqrt[a + b*ArcCosh[-1 + d*x^2]]/(Sqrt[2]*Sqrt[b])]*(Cosh[a/(2*b)] - Sinh[a/(2*b)]))/(b^(3/2)*d*x) + (Sqrt[Pi/2]*Cosh[ArcCosh[-1 + d*x^2]/2]*Erf[Sqrt[a + b*ArcCosh[-1 + d*x^2]]/(Sqrt[2]*Sqrt[b])]*(Cosh[a/(2*b)] + Sinh[a/(2*b)]))/(b^(3/2)*d*x)", // + "(-Sqrt[d*x^2]*Sqrt[-2+d*x^2])/(b*d*x*Sqrt[a+b*ArcCosh[-1+d*x^2]])+(E^(a/(2*b))*Sqrt[Pi/\n" + + "2]*Cosh[ArcCosh[-1+d*x^2]/2]*Erf[Sqrt[a+b*ArcCosh[-1+d*x^2]]/(Sqrt[2]*Sqrt[b])])/(b^(\n" + + "3/2)*d*x)+(Sqrt[Pi/2]*Cosh[ArcCosh[-1+d*x^2]/2]*Erfi[Sqrt[a+b*ArcCosh[-1+d*x^2]]/(Sqrt[\n" + + "2]*Sqrt[b])])/(b^(3/2)*d*E^(a/(2*b))*x)", // 5886); } From 18208255ab42e93e6487daa6440c1dce2d64f83f Mon Sep 17 00:00:00 2001 From: axexlck Date: Sun, 19 Jul 2026 17:10:42 +0200 Subject: [PATCH 8/8] Disable spawning JUnit processes --- symja_android_library/pom.xml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/symja_android_library/pom.xml b/symja_android_library/pom.xml index ed33ec4514..d9ba481e90 100644 --- a/symja_android_library/pom.xml +++ b/symja_android_library/pom.xml @@ -557,9 +557,10 @@ **/*TestCase.java + core 1C true + -->