Merge two powers of one base in a product (#781), and guard a zero integration rate (#785)

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

@@ -11,8 +11,53 @@ namespace AngouriMath.Functions.Algebra
 {
     internal static class IntegralPatterns
     {
+        /// <summary>
+        /// The argument through which each rule below reads its integrand's dependence on
+        /// <paramref name="x"/>, or <see langword="null"/> for a shape none of them matches.
+        /// </summary>
+        /// <remarks>
+        /// The conditions on the two two-argument nodes are the ones their own rules carry:
+        /// a base that mentions <c>x</c> is a different integral, not this one.
+        /// </remarks>
+        private static Entity? RateBearingArgument(Entity expr, Entity.Variable x) => expr switch
+        {
+            Entity.Sinf(var arg) => arg,
+            Entity.Cosf(var arg) => arg,
+            Entity.Secantf(var arg) => arg,
+            Entity.Cosecantf(var arg) => arg,
+            Entity.Tanf(var arg) => arg,
+            Entity.Cotanf(var arg) => arg,
+            Entity.Absf(var arg) => arg,
+            Entity.Signumf(var arg) => arg,
+            Entity.Arcsinf(var arg) => arg,
+            Entity.Arccosf(var arg) => arg,
+            Entity.Arctanf(var arg) => arg,
+            Entity.Arccotanf(var arg) => arg,
+            Entity.Logf(var @base, var arg) when !@base.ContainsNode(x) => arg,
+            Entity.Powf(var @base, var power) when !@base.ContainsNode(x) => power,
+            _ => null
+        };
+
         internal static Entity? TryStandardIntegrals(Entity expr, Entity.Variable x) => expr switch
         {
+            // Every rule below divides by the linear rate of its integrand's argument, and
+            // an argument that mentions x without depending on it has a rate of zero. That
+            // put a literal division by zero into the answer, so `e ^ (x + -x)` integrated
+            // to `e ^ (x + -x) / (0 * ln(e))`, which evaluates to NaN. Each of those
+            // integrands depends on x only through that argument, so a zero rate means the
+            // integrand is a constant and integrates to itself times x.
+            //
+            // The rate has to be decidably zero: a symbolic one, as in sin(a * x), is not,
+            // and answering sin(a * x) * x there would be wrong for every non-zero a.
+            //
+            // Reachable from ordinary input once two powers of one base are gathered:
+            // e^x * e^(-x) becomes e^(x + -x), whose rate is zero while `x + -x` is still
+            // written out. https://github.com/asc-community/AngouriMath/issues/785
+            _ when RateBearingArgument(expr, x) is { } constant
+                   && TreeAnalyzer.TryGetPolyLinear(constant, x, out var rate, out _)
+                   && rate.Evaled is Entity.Number.Complex { IsZero: true }
+                => expr * x,
+
             Entity.Sinf(var arg) when
                 TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _)  =>
                     -MathS.Cos(arg) / a,

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

@@ -46,9 +46,33 @@ namespace AngouriMath.Functions.Algebra
 {
     internal static partial class Integration
     {
+        /// <summary>
+        /// Brings the two powers of one base in a product together, so that an integrand
+        /// which is a power in disguise is recognised as one.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// This has to happen on every recursive call rather than once on the way in. The
+        /// shape is produced by the integrator itself -- distributing a product over a sum
+        /// turns <c>sin(x)^4 * (5 - 6*sin(x)^2)</c> into a term <c>sin(x)^4 * (-6) * sin(x)^2</c>
+        /// -- so normalising only the caller's input would miss every case the integrator
+        /// generates for itself. https://github.com/asc-community/AngouriMath/issues/781
+        /// </para>
+        /// <para>
+        /// Deliberately not followed by <c>InnerSimplified</c>: that rewrites <c>x^(-2)</c>
+        /// back into <c>1/x^2</c>, and <see cref="IndefiniteIntegralSolver.SolveAsPolynomialTerm"/>
+        /// rewrites a <c>1/x^n</c> it is handed into <c>Pow(x, -n)</c>, so the pair recurse
+        /// into each other until the stack runs out. The gathering folds the exponents it
+        /// builds and leaves the rest of the tree alone.
+        /// </para>
+        /// </remarks>
+        private static Entity Normalized(Entity expr) =>
+            expr.Replace(Patterns.GatherPowersOfOneBase);
+
         /// <summary>Does not add the constant of integration because this is called recursively.</summary>
         internal static Entity? ComputeIndefiniteIntegral(Entity expr, Entity.Variable x, bool integrateByParts = true)
         {
+            expr = Normalized(expr);
             if (!expr.ContainsNode(x)) return expr * x; // base case, handle here
             if ((IntegralPatterns.TryStandardIntegrals(expr, x)) is { } answer) return answer;
             // The flag has to be handed on. Every one of these recurses, and a solver that

Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.Power.cs

@@ -143,6 +143,110 @@ when ReduceRadical(radicand, power) is { } reduced => reduced,
             _ => x
         };
 
+        /// <summary>
+        /// Gathers the factors of a product that are powers of one base, wherever they sit
+        /// in it: <c>a^n * c * a^m</c> becomes <c>a^(n+m) * c</c>.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// <see cref="PowerRules"/> already has <c>{}^n * {}^m = {}^(n+m)</c>, but it pairs
+        /// two sibling nodes, and a product is a tree rather than a list -- in
+        /// <c>sin(x)^4 * (-6) * sin(x)^2</c> the constant sits between the two powers, so
+        /// they are never siblings and the rule never fires. Full <see cref="Entity.Simplify"/>
+        /// gets these because it reassociates and sorts the factors first; nothing that
+        /// normalises without sorting does.
+        /// </para>
+        /// <para>
+        /// <c>a^n * a^m = a^(n+m)</c> needs no condition: with <c>a^n</c> read as
+        /// <c>e^(n Log a)</c> on the principal branch, the two sides are
+        /// <c>e^(n Log a) * e^(m Log a)</c> and <c>e^((n+m) Log a)</c>, which are equal for
+        /// every complex <c>n</c> and <c>m</c>. This is what makes it unlike
+        /// <c>(a^b)^c = a^(b*c)</c>, which moves the branch and is guarded above.
+        /// The one point it does not cover is <c>a = 0</c> with a negative exponent, where
+        /// <c>0^2 * 0^(-1)</c> is undefined and <c>0^1</c> is 0 -- so this is applied to
+        /// integrands, where an antiderivative may differ on a measure-zero set, and not in
+        /// <see cref="PowerRules"/>. https://github.com/asc-community/AngouriMath/issues/781
+        /// </para>
+        /// </remarks>
+        internal static Entity GatherPowersOfOneBase(Entity x)
+        {
+            if (x is not (Mulf or Divf))
+                return x;
+
+            // Exponent null marks a factor that is carried through untouched.
+            var factors = new List<(Entity Base, Entity? Exponent)>();
+            var merged = false;
+            foreach (var factor in Mulf.LinearChildren(x))
+            {
+                var (@base, exponent) = Decompose(factor);
+                // Numeric factors are left to evaluation, which folds 2 * 2 into 4. Gathering
+                // them here would write it as 2^2 and call that progress.
+                if (@base is Number)
+                {
+                    factors.Add((factor, null));
+                    continue;
+                }
+                var found = false;
+                for (var i = 0; i < factors.Count; i++)
+                    if (factors[i].Exponent is { } already && factors[i].Base == @base)
+                    {
+                        // Only the exponent is folded, never the product around it.
+                        // InnerSimplified on the whole expression rewrites x^(-2) back into
+                        // 1/x^2, and SolveAsPolynomialTerm turns a 1/x^n it is handed into
+                        // Pow(x, -n) again -- the two normalisations chase each other until
+                        // the stack runs out.
+                        factors[i] = (@base, (already + exponent).InnerSimplified);
+                        merged = found = true;
+                        break;
+                    }
+                if (!found)
+                    factors.Add((@base, exponent));
+            }
+
+            // Rebuilding unconditionally would rewrite every quotient in the tree as a
+            // negative power, since LinearChildren flattens a / b into a * b^(-1).
+            if (!merged)
+                return x;
+
+            Entity? result = null;
+            foreach (var (@base, exponent) in factors)
+            {
+                var factor = exponent switch
+                {
+                    null => @base,
+                    Integer(1) => @base,
+                    _ => new Powf(@base, exponent)
+                };
+                result = result is null ? factor : result * factor;
+            }
+            return result ?? x;
+        }
+
+        /// <summary>
+        /// Reads a factor of a product as a base and an exponent.
+        /// </summary>
+        /// <remarks>
+        /// The nested case is not cosmetic: <see cref="Mulf.LinearChildren"/> writes a
+        /// divisor as <c>(...)^(-1)</c>, so <c>x^3 / x^2</c> arrives as
+        /// <c>x^3 * (x^2)^(-1)</c>, whose second factor has base <c>x^2</c> rather than
+        /// <c>x</c> and would be read as an unrelated base. Unwrapping it is
+        /// <c>(a^b)^n = a^(b*n)</c>, which holds for a whole <c>n</c> whatever the sign of
+        /// <c>a</c> -- the same guard the <c>({}^{})^{}</c> rule above carries, and for the
+        /// same reason. https://github.com/asc-community/AngouriMath/issues/752
+        /// </remarks>
+        private static (Entity Base, Entity Exponent) Decompose(Entity factor)
+        {
+            if (factor is not Powf(var @base, var exponent))
+                return (factor, 1);
+            while (@base is Powf(var inner, var innerExponent)
+                   && (exponent is Integer || inner.Evaled is Real { IsPositive: true }))
+            {
+                exponent = innerExponent * exponent;
+                @base = inner;
+            }
+            return (@base, exponent);
+        }
+
         /// <summary>
         /// Largest divisor tried when reducing a radical. Reducing sqrt(n) exactly would
         /// mean factoring n, which is not something a simplification pass can afford to

Sources/Tests/UnitTests/Calculus/ConstantRateIntegrandTest.cs

@@ -0,0 +1,97 @@
+//
+// Copyright (c) 2019-2022 Angouri.
+// AngouriMath is licensed under MIT.
+// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
+// Website: https://am.angouri.org.
+//
+
+using AngouriMath;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Calculus
+{
+    /// <summary>
+    /// Every standard-integral rule reads a linear rate out of its integrand's argument and
+    /// divides by it. An argument that mentions <c>x</c> without depending on it has a rate
+    /// of zero, and the division was written into the answer unguarded:
+    /// <code>
+    ///     e ^ (x + -x)  ->  e ^ (x + -x) / (0 * ln(e))    -- NaN, where it is x
+    /// </code>
+    /// https://github.com/asc-community/AngouriMath/issues/785
+    /// Found while fixing https://github.com/asc-community/AngouriMath/issues/781, whose
+    /// gathering of two powers of one base turns <c>e^x * e^(-x)</c> into exactly this
+    /// shape -- but the defect predates it and is reachable by writing the exponent out.
+    /// </summary>
+    public sealed class ConstantRateIntegrandTest
+    {
+        /// <summary>
+        /// A zero rate means the integrand is a constant, so the antiderivative is that
+        /// constant times x. Checked by differentiating back rather than by its printed
+        /// form, which still carries the unreduced argument.
+        /// </summary>
+        [Theory]
+        [InlineData("e ^ (x + -x)", "1")]
+        [InlineData("sin(x + -x)", "0")]
+        [InlineData("cos(x - x)", "1")]
+        [InlineData("2 ^ (x - x)", "1")]
+        public void AZeroRateIntegrandIsAConstant(string integrand, string value)
+        {
+            var antiderivative = integrand.Integrate("x");
+            Assert.False(antiderivative is Entity.Integralf,
+                $"{integrand} was declined: {antiderivative.Stringize()}");
+            var difference = (antiderivative.Differentiate("x") - value.ToEntity()).Simplify();
+            while (difference is Entity.Providedf(var inner, _)) difference = inner;
+            Assert.Equal(Entity.Number.Integer.Create(0), difference);
+        }
+
+        /// <summary>
+        /// The guard must fire on a zero rate and on nothing else -- a rate of one or of
+        /// any other constant is what the rules below it are for, and reading those as
+        /// constants would replace every one of these answers with the integrand times x.
+        /// </summary>
+        [Theory]
+        [InlineData("e ^ (2 * x)")]
+        [InlineData("sin(3 * x + 1)")]
+        [InlineData("cos(x)")]
+        [InlineData("ln(2 * x)")]
+        public void ANonZeroRateStillIntegratesByItsRule(string integrand)
+        {
+            var antiderivative = integrand.Integrate("x");
+            Assert.False(antiderivative is Entity.Integralf,
+                $"{integrand} was declined: {antiderivative.Stringize()}");
+            var difference = (antiderivative.Differentiate("x") - integrand.ToEntity()).Simplify();
+            while (difference is Entity.Providedf(var inner, _)) difference = inner;
+            Assert.Equal(Entity.Number.Integer.Create(0), difference);
+        }
+
+        /// <summary>
+        /// <c>abs</c> takes the guard too, and is checked on the antiderivative rather than
+        /// by differentiating it back: the library differentiates <c>abs(u)</c> as
+        /// <c>sgn(u) * u'</c> and answers <c>NaN</c> at <c>u = 0</c>, so the round trip
+        /// says nothing about this rule either way.
+        /// </summary>
+        [Fact]
+        public void AZeroRateAbsIntegrandIsZero()
+        {
+            var antiderivative = "abs(x + -x)".Integrate("x");
+            Assert.False(antiderivative is Entity.Integralf);
+            Assert.Equal(Entity.Number.Integer.Create(0),
+                (antiderivative - "C".ToEntity()).Simplify());
+        }
+
+        /// <summary>
+        /// A symbolic rate is not decidably zero, so the guard must leave it to the rules --
+        /// answering <c>sin(a*x) * x</c> here would be wrong for every non-zero a.
+        /// </summary>
+        [Fact]
+        public void ASymbolicRateIsLeftToTheRules()
+        {
+            var antiderivative = "sin(a * x)".Integrate("x");
+            Assert.False(antiderivative is Entity.Integralf);
+            var difference = (antiderivative.Differentiate("x") - "sin(a * x)".ToEntity()).Simplify();
+            while (difference is Entity.Providedf(var inner, _)) difference = inner;
+            Assert.Equal(Entity.Number.Integer.Create(0), difference);
+        }
+    }
+}