Defer Evaluating RootSums

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;
@@ -1000,6 +1010,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 +1196,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 <code>(a + b*x)^p</code> with a symbolic index
+     * <code>n</code>, expanded around <code>x==0</code>.
+     *
+     * <p>
+     * The coefficient of <code>x^n</code> in <code>(a + b*x)^p</code> is
+     * <code>b^n * a^(p-n) * Binomial(p, n)</code>. When <code>p</code> is a non-negative integer the
+     * expansion is finite, so the index is bounded by <code>0 &lt;= n &lt;= p</code>; otherwise the
+     * series runs for all <code>n &gt;= 0</code>. Returning this closed Binomial form is a
+     * readability improvement over the equivalent (but opaque) <code>DifferenceRoot</code> the
+     * holonomic engine would otherwise emit for a finite binomial power.
+     *
+     * @return a {@link S#Piecewise} expression, or {@link F#NIL} if <code>function</code> 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

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) {

symja_android_library/matheclipse-core/src/main/java/org/matheclipse/core/integrate/RationalIntegration.java

@@ -53,14 +53,67 @@ public class RationalIntegration {
   private RationalIntegration() {}
 
   /**
-   * Integrate a rational function in <code>x</code>.
+   * Controls whether the logarithmic part may emit a {@link F#RootSum} antiderivative for an
+   * irreducible denominator factor of degree &gt;= 5.
+   */
+  public enum RootSumMode {
+    /**
+     * Irreducible factors of degree &gt;= 5 return {@link F#NIL} (defer to the Rubi rules), exactly
+     * like the degree 3/4 factors already do.
+     */
+    DEFER,
+    /** Irreducible factors of degree &gt;= 5 produce a {@link F#RootSum} antiderivative. */
+    EMIT
+  }
+
+  /**
+   * Integrate a rational function in <code>x</code>, emitting a {@link F#RootSum} antiderivative for
+   * irreducible denominator factors of degree &gt;= 5 (i.e. {@link RootSumMode#EMIT}).
    *
    * @param integrand the integrand
    * @param x the integration variable
    * @param engine the evaluation engine
    * @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 <code>x</code>.
+   *
+   * <p>
+   * In {@link RootSumMode#DEFER} mode, when the computed antiderivative is <em>essentially a bare
+   * {@link F#RootSum}</em> (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 <em>not</em> 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 &gt;= 5.
+   */
+  private static IExpr integrateRationalFunction(IExpr integrand, IExpr x, EvalEngine engine) {
     if (!Config.INTEGRATE_ALGORITHM_RATIONAL) {
       return F.NIL;
     }

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 <code>x</code> is canonicalized: with a numeric step it is factored
+   * (matching wolframscript, e.g. <code>x^2 + x -&gt; 2*(1 + x)</code>, and rational summands combine
+   * over a common denominator), with a symbolic step it is only expanded. A difference of unknown
+   * functions such as <code>f(x + 1) - f(x)</code> 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.

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 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 <code>f(k + step) / f(k)</code>.
    */
-  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

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