The binomial differential, in the case Chebyshev's criterion makes a polynomial

Sources/.editorconfig

@@ -345,3 +345,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed
 
 [Tests/UnitTests/Calculus/TrigonometricReductionTest.cs]
 file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n
+
+[Tests/UnitTests/Calculus/BinomialDifferentialTest.cs]
+file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n

Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs

@@ -909,6 +909,208 @@ Entity Carried(int numerator, int denominator, int nextSine, int nextCosine)
                  + Carried(sinePower - 1, sum, sinePower - 2, cosinePower);
         }
 
+        /// <summary>
+        /// A <b>binomial differential</b> <c>x^m (a + b x^n)^(p/q)</c>, in the case Chebyshev's
+        /// criterion makes a polynomial: where <c>(m + 1)/n</c> is a whole number of at least one.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// <c>x^2 sqrt(1 + x^3)</c> came out and <c>x^5 sqrt(1 + x^3)</c> did not, and the two are
+        /// the same substitution. Under <c>u = (a + b x^n)^(1/q)</c> the first leaves a monomial
+        /// in <c>u</c>, which the general substitution finds because <c>x^(n-1)</c> is <c>du</c>
+        /// up to a constant; the second leaves a <b>polynomial</b>, which it does not, because
+        /// what multiplies <c>du</c> there is <c>x^3</c> rather than a constant.
+        /// </para>
+        /// <para>
+        /// Writing <c>s = (m + 1)/n</c> and substituting:
+        /// </para>
+        /// <code>
+        ///     x^n = (u^q - a)/b       x^m dx = (q/(n b)) ((u^q - a)/b)^(s-1) u^(q-1) du
+        ///     int x^m (a + b x^n)^(p/q) dx = (q/(n b)) int ((u^q - a)/b)^(s-1) u^(p+q-1) du
+        /// </code>
+        /// <para>
+        /// which is a polynomial in <c>u</c> exactly when <c>s</c> is a whole number of at least
+        /// one — the first of Chebyshev's three cases. Expanded by the binomial theorem and
+        /// integrated term by term, so the rule is <b>closed</b> and asks the integrator nothing.
+        /// </para>
+        /// <para>
+        /// <b>The other two cases are not here.</b> <c>p/q</c> whole is a whole power of a
+        /// polynomial, which expanding already answers; and <c>s + p/q</c> whole wants
+        /// <c>u = ((a + b x^n)/x^n)^(1/q)</c>, which leaves a polynomial over a power of
+        /// <c>u^q - b</c> — a rational function rather than a polynomial, so a search rather than
+        /// a closed answer. Chebyshev also proved there is no fourth case: outside those three the
+        /// integrand has no elementary antiderivative at all, which is worth knowing before
+        /// anyone goes looking.
+        /// </para>
+        /// <para>
+        /// <b>Asked, not volunteered</b>:
+        /// <a href="https://github.com/asc-community/AngouriMath/issues/1265">#1265</a>.
+        /// </para>
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveABinomialDifferential(Entity expr, Entity.Variable x)
+        {
+            if (!Integration.AnsweringTheQuestionAsked)
+                return null;
+            if (!TryReadABinomialDifferential(expr, x, out var power, out var exponent,
+                    out var inner, out var free, out var leading, out var factor))
+                return null;
+
+            // s = (m + 1)/n, whole and at least one, which is what makes the expansion finite.
+            if ((power + 1) % inner != 0)
+                return null;
+            var s = (power + 1) / inner;
+            if (s < 1)
+                return null;
+
+            var q = exponent.Denominator.ToInt32Checked();
+            var p = exponent.Numerator.ToInt32Checked();
+            var a = Number.Rational.Create(free);
+            var b = Number.Rational.Create(leading);
+
+            // (q/(n b^s)) * int (u^q - a)^(s-1) u^(p+q-1) du, the binomial expanded.
+            Entity total = 0;
+            var binomial = EInteger.One;
+            var u = Variable.CreateUnique(expr, "u_binom");
+            for (var i = 0; i <= s - 1; i++)
+            {
+                // Term i of (u^q - a)^(s-1) is C(s-1, i) u^(q i) (-a)^(s-1-i).
+                var raised = q * i + p + q;
+                if (raised == 0)
+                    return null;   // the power rule would divide by zero; not a polynomial after all
+                total += Number.Integer.Create(binomial)
+                       * MathS.Pow(-a, Number.Integer.Create(s - 1 - i))
+                       * MathS.Pow(u, Number.Integer.Create(raised)) / Number.Integer.Create(raised);
+                binomial = binomial * (s - 1 - i) / (i + 1);
+            }
+
+            var outside = Number.Integer.Create(q)
+                        / (Number.Integer.Create(inner) * MathS.Pow(b, Number.Integer.Create(s)));
+            var bracket = (a + b * MathS.Pow(x, Number.Integer.Create(inner))).InnerSimplified;
+            var back = MathS.Pow(bracket, Number.Rational.Create(EInteger.One, exponent.Denominator));
+            return (factor * outside * total.Substitute(u, back)).InnerSimplified;
+        }
+
+        /// <summary>
+        /// Reads <paramref name="expr"/> as a rational multiple of <c>x^m (a + b x^n)^(p/q)</c>,
+        /// with <c>q</c> above one and <c>n</c> at least two.
+        /// </summary>
+        /// <remarks>
+        /// A whole exponent on the bracket is a polynomial and wants expanding rather than this;
+        /// <c>n = 1</c> is a radical of something linear, which
+        /// <see cref="SolveByLinearRadicalSubstitution"/> answers in its own terms and more
+        /// shortly.
+        /// </remarks>
+        private static bool TryReadABinomialDifferential(
+            Entity expr, Entity.Variable x, out int power, out ERational exponent,
+            out int inner, out ERational free, out ERational leading, out Entity factor)
+        {
+            power = 0;
+            exponent = ERational.Zero;
+            inner = 0;
+            free = ERational.Zero;
+            leading = ERational.Zero;
+            factor = 1;
+
+            Entity? bracket = null;
+            var exponentFound = ERational.Zero;
+            var powerFound = 0;
+            Entity constantFactor = 1;
+            if (!Read(expr, 1) || bracket is null)
+                return false;
+
+            var lowest = exponentFound.ToLowestTerms();
+            if (lowest.Denominator.Equals(EInteger.One)
+                || !lowest.Numerator.CanFitInInt32() || !lowest.Denominator.CanFitInInt32())
+                return false;
+
+            // `a + b x^n`: one term free of the variable and one monomial in it.
+            Entity? constantPart = null;
+            Entity? monomial = null;
+            foreach (var term in Sumf.LinearChildren(bracket))
+                if (!term.ContainsNode(x))
+                    constantPart = constantPart is null ? term : constantPart + term;
+                else if (monomial is null)
+                    monomial = term;
+                else
+                    return false;
+            if (constantPart is null || monomial is null)
+                return false;
+            if (!TryReadAMonomial(monomial, x, out var degree, out var coefficient) || degree < 2)
+                return false;
+            if (constantPart.Evaled is not Number.Rational constantValue
+                || coefficient.Evaled is not Number.Rational coefficientValue
+                || coefficientValue.ERational.IsZero)
+                return false;
+
+            power = powerFound;
+            exponent = lowest;
+            inner = degree;
+            free = constantValue.ERational;
+            leading = coefficientValue.ERational;
+            factor = constantFactor;
+            return true;
+
+            bool Read(Entity node, int multiplicity)
+            {
+                switch (node)
+                {
+                    case Variable v when v == x:
+                        powerFound += multiplicity;
+                        return true;
+                    case Mulf(var left, var right):
+                        return Read(left, multiplicity) && Read(right, multiplicity);
+                    case Divf(var above, var below):
+                        return Read(above, multiplicity) && Read(below, -multiplicity);
+                    case Powf(var @base, Number.Integer whole)
+                        when @base == x && whole.EInteger.CanFitInInt32():
+                        powerFound += multiplicity * whole.EInteger.ToInt32Checked();
+                        return true;
+                    case Powf(var @base, Number.Rational raised) when @base.ContainsNode(x):
+                        if (bracket is not null && bracket != @base)
+                            return false;
+                        bracket = @base;
+                        exponentFound += ERational.FromInt32(multiplicity) * raised.ERational;
+                        return true;
+                    case Number.Rational rational when node is not Powf:
+                        constantFactor = multiplicity > 0
+                            ? constantFactor * MathS.Pow(rational, multiplicity)
+                            : constantFactor / MathS.Pow(rational, -multiplicity);
+                        return true;
+                    default:
+                        return false;
+                }
+            }
+        }
+
+        /// <summary>Reads <c>c * x^n</c>, giving the degree and the coefficient.</summary>
+        private static bool TryReadAMonomial(Entity term, Entity.Variable x, out int degree, out Entity coefficient)
+        {
+            degree = 0;
+            coefficient = 1;
+            switch (term)
+            {
+                case Variable v when v == x:
+                    degree = 1;
+                    return true;
+                case Powf(var @base, Number.Integer whole) when @base == x && whole.EInteger.CanFitInInt32():
+                    degree = whole.EInteger.ToInt32Checked();
+                    return true;
+                case Mulf(var left, var right) when !left.ContainsNode(x):
+                    if (!TryReadAMonomial(right, x, out degree, out var fromRight))
+                        return false;
+                    coefficient = left * fromRight;
+                    return true;
+                case Mulf(var left, var right) when !right.ContainsNode(x):
+                    if (!TryReadAMonomial(left, x, out degree, out var fromLeft))
+                        return false;
+                    coefficient = right * fromLeft;
+                    return true;
+                default:
+                    return false;
+            }
+        }
+
         /// <summary>
         /// A power of the variable times an odd half-power of a quadratic without a linear term —
         /// <c>x^m (a + b x^2)^(k/2)</c> with <c>k</c> odd — turned into a power of the sine times

Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs

@@ -430,6 +430,10 @@ private static Entity Normalized(Entity expr) =>
             // cosine rather than by rationalising it. After the linear radical, which answers a
             // root of something linear in its own terms and more shortly.
             if ((answer = IndefiniteIntegralSolver.SolveARadicalOfAQuadraticAsTrigonometric(expr, x)) is { }) return answer;
+            // And the binomial differential, a power of the variable beside a fractional power of
+            // `a + b x^n`. After the quadratic, which is the case `n = 2` and answers it through
+            // the trigonometric substitution -- a shorter answer than a root of a root.
+            if ((answer = IndefiniteIntegralSolver.SolveABinomialDifferential(expr, x)) is { }) return answer;
             // Product-to-sum among the rewrites rather than before them, because a product of
             // trigonometric functions of *equal* arguments is a power and wants a different tool;
             // this only fires where the arguments differ, which is exactly what every substitution

Sources/Tests/UnitTests/Calculus/BinomialDifferentialTest.cs

@@ -0,0 +1,136 @@
+//
+// Copyright (c) 2019-2026 Angouri.
+// AngouriMath is licensed under MIT.
+// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
+// Website: https://am.angouri.org.
+//
+
+using System;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Calculus
+{
+    /// <summary>
+    /// The binomial differential <c>x^m (a + b x^n)^(p/q)</c>, in the case Chebyshev's criterion
+    /// makes a polynomial.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    /// <remarks>
+    /// <para>
+    /// <c>x^2 sqrt(1 + x^3)</c> came out and <c>x^5 sqrt(1 + x^3)</c> did not, under the same
+    /// substitution. <c>u = (a + b x^n)^(1/q)</c> leaves a monomial in <c>u</c> for the first —
+    /// which the general substitution finds, because <c>x^(n-1)</c> is <c>du</c> up to a constant
+    /// — and a polynomial for the second, which it does not.
+    /// </para>
+    /// <para>
+    /// The condition is <c>s = (m + 1)/n</c> being a whole number of at least one. Checked by
+    /// differentiating back and comparing at points.
+    /// </para>
+    /// </remarks>
+    [Trait("Area", "Calculus")]
+    public sealed class BinomialDifferentialTest
+    {
+        private static readonly double[] Points = { 0.21, 0.35, 0.52, 0.71, 0.83 };
+
+        private static void DifferentiatesBack(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("NaN", integral.Stringize());
+            Assert.DoesNotContain("integral(", integral.Stringize());
+
+            var derivative = integral.Substitute("C", 0).Differentiate("x");
+            var original = integrand.ToEntity();
+            var compared = 0;
+            foreach (var at in Points)
+            {
+                var got = derivative.Substitute("x", at).EvalNumerical();
+                var want = original.Substitute("x", at).EvalNumerical();
+                if (got.IsNaN || want.IsNaN)
+                    continue;
+                compared++;
+                var difference = Math.Abs((double)(got - want).RealPart)
+                               + Math.Abs((double)(got - want).ImaginaryPart);
+                var scale = Math.Max(1.0, Math.Abs((double)want.RealPart));
+                Assert.True(difference / scale < 1e-9,
+                    $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, "
+                    + $"where the integrand is {want}");
+            }
+            Assert.True(compared >= 4,
+                $"only {compared} of {Points.Length} points were comparable for {integrand}, "
+                + "so this asserts almost nothing");
+        }
+
+        /// <summary>
+        /// <c>s = 2</c> and above, where the substitution leaves a polynomial rather than a
+        /// monomial. These are the ones that were declined.
+        /// </summary>
+        [Theory]
+        [InlineData("x^5*sqrt(1 + x^3)")]
+        [InlineData("x^5*(1 + x^3)^(2/3)")]
+        [InlineData("x^5/(1 + x^3)^(1/3)")]
+        [InlineData("x^3*(1 + x^2)^(1/3)")]
+        [InlineData("x^7*(1 + x^4)^(1/2)")]
+        [InlineData("x^8*sqrt(1 + x^3)")]
+        public void WhereTheSubstitutionLeavesAPolynomial(string integrand)
+            => DifferentiatesBack(integrand);
+
+        /// <summary>
+        /// <c>s = 1</c>, the monomial case, which the general substitution already answered. It
+        /// reaches this rule first now, so the answers have to be right from here too.
+        /// </summary>
+        [Theory]
+        [InlineData("x^2*sqrt(1 + x^3)")]
+        [InlineData("x^2*(1 + x^3)^(1/3)")]
+        [InlineData("x^2/(1 + x^3)^(1/3)")]
+        [InlineData("x*sqrt(1 + x^2)")]
+        [InlineData("x*(1 + x^2)^(1/3)")]
+        [InlineData("x^3/sqrt(1 + x^4)")]
+        [InlineData("x^2*(4 - x^3)^(1/3)")]
+        public void TheMonomialCaseIsUnchangedInSubstance(string integrand)
+            => DifferentiatesBack(integrand);
+
+        /// <summary>
+        /// A rational factor in front and a negative coefficient inside, both of which the read
+        /// carries rather than declining over.
+        /// </summary>
+        [Theory]
+        [InlineData("3*x^5*sqrt(1 + x^3)")]
+        [InlineData("x^5*sqrt(2 - 3*x^3)/2")]
+        [InlineData("x^5*(1 + 2*x^3)^(1/3)")]
+        public void AFactorAndACoefficient(string integrand) => DifferentiatesBack(integrand);
+
+        /// <summary>
+        /// Outside Chebyshev's first case. <c>s</c> not whole is either his second or third case
+        /// or — and this is the part worth knowing before anyone goes looking — no elementary
+        /// antiderivative at all, which he proved. These must be declined or answered by another
+        /// rule, never wrong.
+        /// </summary>
+        [Theory]
+        [InlineData("sqrt(1 + x^3)")]
+        [InlineData("x*sqrt(1 + x^3)")]
+        [InlineData("x^3*sqrt(1 + x^3)")]
+        [InlineData("1/(x*(1 + x^3)^(1/3))")]
+        public void OutsideTheFirstCase(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("NaN", integral.Stringize());
+            if (integral.Stringize().Contains("integral("))
+                return;   // unanswered is a legitimate verdict; a wrong answer is not
+            DifferentiatesBack(integrand);
+        }
+
+        /// <summary>
+        /// What the read must refuse: a whole exponent on the bracket, which is a polynomial and
+        /// wants expanding; a linear bracket, which the linear radical rule answers more shortly;
+        /// and a quadratic one, which the trigonometric substitution answers more shortly still.
+        /// </summary>
+        [Theory]
+        [InlineData("x^2*(1 + x^3)^2")]
+        [InlineData("x*sqrt(1 + x)")]
+        [InlineData("sqrt(1 + 2*x)")]
+        [InlineData("x^3*sqrt(1 + x^2)")]
+        [InlineData("1/(1 + x^2)^(3/2)")]
+        public void TheNeighboursAreUntouched(string integrand) => DifferentiatesBack(integrand);
+    }
+}