Answer the limits that come out as infinity minus infinity (#682)

lim x -> +oo of e^x - x, of x - ln(x) and of sqrt(x^2 + x) - x were all left
unevaluated. Every solver reads either a polynomial or a substitution of +oo into
the expression, and a difference of two parts that each grow without bound is
neither: substituting gives oo - oo, which says nothing.

Two readings settle most of them. Where one part outgrows the other the answer is
that part's, which is what the ratio of the two says -- x / e^x tends to 0, so
e^x - x tends to +oo. Where the two grow alike the ratio tends to 1 and says
nothing, and then a - b = (a^2 - b^2) / (a + b), an identity wherever a + b is
not zero, which on the way to +oo it is not. It is worth using where squaring
removes a root: the leading terms of the numerator cancel and what is left is an
ordinary quotient, so sqrt(x^2 + x) - x becomes x / (sqrt(x^2 + x) + x).

Neither reading works on a root of a polynomial, because no solver can say
anything about one. sqrt(x^2 + x) is now written x * sqrt(1 + 1/x), which puts
its growth in plain sight and leaves a factor that +oo can simply be substituted
into. That is sound in the one direction this is used in: P = x^d * (P / x^d),
and (uv)^r = u^r v^r wants u positive, which x^d is on the way to +oo, and a
destination of -oo has already been turned into +oo by substituting -x for x.
It also answers sqrt(x^2 + x) / x, x / sqrt(x^2 + 1) and sqrt(x^2 - x) / x, none
of which had an answer before. An odd degree is not covered: the growth of
sqrt(x^3 + x) is x^(3/2), and l'Hopital's rule is stopped from working on a
quotient containing that.

All of this is reached only once the expression as written has been found to have
no answer, so nothing that had one is asked twice, and it sits beside l'Hopital's
rule rather than in the descent that visits every part of every expression --
putting it there instead cost the test suite four times its running time.

Simplification is asked again as well, since it can hand back an expression of
another kind altogether and which parts a limit is broken into is decided by that
kind: sqrt(x^2 + 1) / sqrt(x^2 + 3x) is broken up as a quotient, while what
simplifying gives is the single root sqrt((x^2 + 1) / (x^2 + 3x)), whose argument
can be read straight off. That one took 2.2 seconds to reach no answer and now
takes 14 ms to reach 1.

3937 tests, corpus 93/117 from 91 with nothing wrong, in error or timing out, and
the oo - oo family complete. 18 of the 27 new tests fail without the change; the
9 that pass are the ones pinning what already worked.

Author: Rafael-SOWNet

Sources/AngouriMath/Functions/Continuous/Limits/Solvers/Solvers.Definition.cs

@@ -75,6 +75,14 @@ private static Entity ExpandLogarithm(Entity expr)
                     // reciprocal, and a NaN from it would claim the limit does not exist.
                     if (ApplylHopitalRule(expr, x, dest) is { } lhopital && lhopital.Evaled != MathS.NaN)
                         return lhopital;
+                    // The rewrites are worth their cost only where there is no answer without
+                    // them, and each of them costs an expansion or a simplification of the
+                    // whole expression. The descent visits every part of every expression, so
+                    // rather than pay for them at each of those parts on every limit ever
+                    // taken, the descent is walked a second time with them turned on, and only
+                    // for an expression that has just been found to have no answer at all.
+                    if (SolveByRewriting(expr, x, dest) is { } rewritten && rewritten.Evaled != MathS.NaN)
+                        return rewritten;
                     return atInfinity;
                 }
                 else if (expr.ComputeLimitDivideEtImpera(x, dest, ApproachFrom.Left) is { } fromLeft

Sources/AngouriMath/Functions/Continuous/Limits/Transformations.cs

@@ -6,9 +6,11 @@
 //
 
 using AngouriMath.Core.Multithreading;
+using PeterO.Numbers;
 using System;
 using System.Linq;
 using static AngouriMath.Entity;
+using static AngouriMath.Entity.Number;
 
 namespace AngouriMath.Functions.Algebra
 {
@@ -129,6 +131,8 @@ private static bool AlreadyBeingDifferentiated(Entity quotient)
         private static bool GrewTooMuch(Entity quotient, Entity applied)
             => applied.Nodes.Count() > quotient.Nodes.Count() + 8;
 
+        [ThreadStatic] private static bool suppresslHopital;
+
         /// <summary>
         /// A product that has a reciprocal factor in it, rewritten as a quotient, or
         /// <see langword="null"/> if it has none. The rule below only reads quotients, so
@@ -169,7 +173,7 @@ private static bool GrewTooMuch(Entity quotient, Entity applied)
         {
             if (lHopitalDepth == 0)
                 lHopitalApplications = 0;
-            if (lHopitalDepth >= MaxlHopitalDepth || lHopitalApplications >= MaxlHopitalApplications)
+            if (lHopitalDepth >= MaxlHopitalDepth || lHopitalApplications >= MaxlHopitalApplications || suppresslHopital)
                 return null;
             // Held for the whole of the rule and not just for the recursive call below, since
             // asking what the two parts tend to is itself a limit that the rule may be applied
@@ -217,6 +221,184 @@ private static bool GrewTooMuch(Entity quotient, Entity applied)
             return null;
         }
 
+        /// <summary>
+        /// How deep one limit may go into rewriting itself. Each of the rewrites below hands
+        /// back another limit to take, and that one is entitled to be rewritten in turn, so
+        /// without a bound the work would multiply.
+        /// </summary>
+        private const int MaxRewriteDepth = 2;
+
+        [ThreadStatic] private static int rewriteDepth;
+
+        /// <summary>
+        /// The readings that need the expression rewritten before any solver can see anything
+        /// in it. Every one of them costs an expansion or a simplification of the whole
+        /// expression, which is why this sits here rather than in the descent that visits each
+        /// of its parts, and why it is reached only once it is settled that the expression as
+        /// written has no answer.
+        /// </summary>
+        private static Entity? SolveByRewriting(Entity expr, Variable x, Entity dest)
+        {
+            if (rewriteDepth >= MaxRewriteDepth)
+                return null;
+            // -oo is +oo with -x written for x, the same substitution the solvers are handed.
+            // Reading the growth of a root off x^d depends on it: x^d is positive in the one
+            // direction only.
+            if (dest.Evaled is Real { IsNegative: true })
+                expr = expr.Substitute(x, -x);
+            var toInfinity = Real.PositiveInfinity;
+            var simplified = expr.Simplify();
+            if (simplified is Providedf(var body, _))
+                simplified = body;
+
+            rewriteDepth++;
+            try
+            {
+                // Simplification can hand back an expression of another kind altogether, and
+                // which parts a limit is broken into is decided by that kind:
+                // sqrt(x^2 + 1) / sqrt(x^2 + 3x) is broken up as a quotient, while what
+                // simplifying gives is the single root sqrt((x^2 + 1) / (x^2 + 3x)), whose
+                // argument can be read straight off.
+                if (simplified.GetType() != expr.GetType()
+                    && Settled(simplified.ComputeLimitDivideEtImpera(x, toInfinity, ApproachFrom.Left)) is { } byShape)
+                    return byShape;
+
+                if (ExtractRadicalGrowth(simplified, x) is { } extracted
+                    && Settled(extracted.ComputeLimitDivideEtImpera(x, toInfinity, ApproachFrom.Left)) is { } byGrowth)
+                    return byGrowth;
+
+                return Settled(SolveAsDifferenceOfInfinities(simplified, x));
+            }
+            finally { rewriteDepth--; }
+
+            // Rewriting brings domain conditions with it -- dividing by x^d is only the same
+            // expression where x is not zero -- and a limit is taken of a continuous
+            // expression regardless of the points where it is undefined, as the solvers
+            // themselves do with the same shape.
+            static Entity? Settled(Entity? limit)
+            {
+                while (limit is Providedf(var inner, _)) limit = inner;
+                return limit is null || limit.Evaled == MathS.NaN ? null : limit;
+            }
+        }
+
+        /// <summary>
+        /// The whole expression with every root of a polynomial rewritten so that its growth is
+        /// a factor of its own -- <c>sqrt(x^2 + x)</c> becomes <c>x * sqrt(1 + 1/x)</c> -- or
+        /// <see langword="null"/> if it has no such root. Every solver reads a polynomial or a
+        /// substitution of +oo, and neither can say anything about a root of a sum, so
+        /// <c>lim x -&gt; +oo sqrt(x^2 + x) / x</c> was left unevaluated while the rewritten
+        /// <c>sqrt(1 + 1/x)</c> is settled by substitution alone.
+        /// </summary>
+        /// <remarks>
+        /// <c>P = x^d * (P / x^d)</c>, and <c>(uv)^r = u^r v^r</c> needs <c>u</c> to be positive,
+        /// which <c>x^d</c> is for every x past some point on the way to +oo. That is the only
+        /// direction this is used in: a destination of -oo has already been turned into +oo by
+        /// substituting -x for x before any of this runs.
+        /// </remarks>
+        private static Entity? ExtractRadicalGrowth(Entity expr, Variable x)
+        {
+            if (!expr.Nodes.Any(IsRootOfASum))
+                return null;
+            var extracted = expr.Replace(RewriteRoot);
+            return extracted == expr ? null : extracted.Simplify();
+
+            bool IsRootOfASum(Entity node)
+                => node is Powf(Sumf or Minusf, Number.Rational and not Number.Integer);
+
+            Entity RewriteRoot(Entity node)
+            {
+                if (!IsRootOfASum(node)
+                    || node is not Powf(var @base, var power)
+                    || !TreeAnalyzer.TryGetPolynomial(@base, x, out var monomials))
+                    return node;
+                var degree = monomials.Keys.Aggregate(EInteger.Zero, EInteger.Max);
+                if (degree.CompareTo(EInteger.One) < 0)
+                    return node;
+                var growth = MathS.Pow(x, Number.Integer.Create(degree));
+                return MathS.Pow(x, Number.Integer.Create(degree) * power) * MathS.Pow((@base / growth).Simplify(), power);
+            }
+        }
+
+        /// <summary>
+        /// How many times a single limit may be broken down as a difference of two divergent
+        /// parts. The conjugate below turns one difference into another, and every part of it
+        /// is a limit in its own right, so without a bound the work would multiply.
+        /// </summary>
+        private const int MaxDifferenceDepth = 2;
+
+        [ThreadStatic] private static int differenceDepth;
+
+        /// <summary>
+        /// Which infinity an already computed limit is, or 0 if it is finite or not a number.
+        /// </summary>
+        private static int InfiniteSign(Entity? limit)
+            => limit?.Evaled is Real { IsFinite: false, IsNaN: false } real ? (real.IsNegative ? -1 : 1) : 0;
+
+        /// <summary>
+        /// oo - oo, which says nothing on its own: whichever of the two grows faster decides
+        /// the answer, and if neither does the difference can still be finite. Two of the
+        /// standard readings are covered -- one part outgrowing the other, and the conjugate
+        /// for a difference containing a root.
+        /// </summary>
+        private static Entity? SolveAsDifferenceOfInfinities(Entity expr, Variable x)
+        {
+            if (differenceDepth >= MaxDifferenceDepth || expr is not (Sumf or Minusf))
+                return null;
+            var terms = Sumf.LinearChildren(expr).ToArray();
+            if (terms.Length != 2)
+                return null;
+            var dest = Real.PositiveInfinity;
+            differenceDepth++;
+            try
+            {
+                var (firstSign, secondSign) =
+                    (InfiniteSign(ComputeLimit(terms[0], x, dest)), InfiniteSign(ComputeLimit(terms[1], x, dest)));
+                if (firstSign == 0 || secondSign != -firstSign)
+                    return null;
+                // Written as minuend - subtrahend with both parts tending to +oo, so that the
+                // reading below does not have to carry the signs around with it.
+                var (minuend, subtrahend) = firstSign > 0
+                    ? (terms[0], -terms[1])
+                    : (terms[1], -terms[0]);
+
+                // The faster growing part decides the answer whenever there is one, which is
+                // what the ratio of the two says: lim x -> +oo e^x - x is +oo because x / e^x
+                // tends to 0, and -oo the other way round. This settles nothing when the two
+                // grow alike, and the ratio then tends to 1 rather than to 0 or to infinity.
+                if (ComputeLimit((subtrahend / minuend).Simplify(), x, dest) is { } ratio)
+                {
+                    if (ratio.Evaled == Number.Integer.Zero)
+                        return Real.PositiveInfinity;
+                    if (InfiniteSign(ratio) > 0)
+                        return Real.NegativeInfinity;
+                }
+
+                // a - b = (a^2 - b^2) / (a + b), an identity wherever a + b is not zero, which
+                // it is not on the way to +oo. It is only an improvement when squaring removes
+                // a root, and then the leading terms cancel in the numerator and what is left
+                // is an ordinary quotient: sqrt(x^2 + x) - x becomes x / (sqrt(x^2 + x) + x).
+                if (!ContainsRadical(minuend, x) && !ContainsRadical(subtrahend, x))
+                    return null;
+                MultithreadingFunctional.ExitIfCancelled();
+                var conjugate = ((minuend * minuend - subtrahend * subtrahend) / (minuend + subtrahend)).Simplify();
+                if (conjugate == expr)
+                    return null;
+                // Without the roots the numerator is an ordinary polynomial and the denominator
+                // is what it was, so the quotient either falls to the solvers directly or it is
+                // no better than what it replaced. Differentiating it repeatedly is a long way
+                // round to the same nothing, and it is most of the cost here.
+                suppresslHopital = true;
+                try { return ComputeLimit(conjugate, x, dest); }
+                finally { suppresslHopital = false; }
+            }
+            finally { differenceDepth--; }
+        }
+
+        private static bool ContainsRadical(Entity expr, Variable x)
+            => expr.Nodes.Any(node =>
+                node is Powf(var @base, Number.Rational and not Number.Integer) && @base.ContainsNode(x));
+
         private static Entity ApplyTrivialTransformations(Entity expr, Variable x, Entity dest, Func<Entity, Entity, Entity> transformation)
             => expr switch
             {

Sources/Tests/UnitTests/Calculus/DifferenceAtInfinityTest.cs

@@ -0,0 +1,98 @@
+//
+// 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 System;
+using System.Threading.Tasks;
+using AngouriMath;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Calculus
+{
+    /// <summary>
+    /// oo - oo, and the roots of polynomials that so often produce it. Every solver reads
+    /// either a polynomial or a substitution of +oo, and a root of a sum is neither, so a
+    /// difference of two of them was left unevaluated however plainly it converged.
+    /// </summary>
+    public sealed class DifferenceAtInfinityTest
+    {
+        // Compared as numbers, since "1/2" parses as a division of two integers rather than as
+        // the rational the limit answers with.
+        private static void AssertLimit(string expression, string destination, string expected) =>
+            Assert.Equal(
+                expected.ToEntity().Evaled,
+                expression.ToEntity().Limit("x", destination.ToEntity()).Evaled);
+
+        // One part outgrowing the other decides the answer, which is what the ratio of the
+        // two says.
+        [Theory]
+        [InlineData("e ^ x - x", "+oo", "+oo")]
+        [InlineData("x - e ^ x", "+oo", "-oo")]
+        [InlineData("x - ln(x)", "+oo", "+oo")]
+        [InlineData("ln(x) - x", "+oo", "-oo")]
+        [InlineData("x ^ 2 - x", "+oo", "+oo")]
+        [InlineData("x - x ^ 2", "+oo", "-oo")]
+        [InlineData("sqrt(x) - ln(x)", "+oo", "+oo")]
+        [InlineData("e ^ x - x ^ 10", "+oo", "+oo")]
+        public void TheFasterGrowingPartDecides(string expression, string destination, string expected) =>
+            AssertLimit(expression, destination, expected);
+
+        // a - b = (a^2 - b^2) / (a + b) removes the roots from the numerator, and what is
+        // left of it after the leading terms cancel is an ordinary quotient.
+        [Theory]
+        [InlineData("sqrt(x ^ 2 + x) - x", "+oo", "1/2")]
+        [InlineData("sqrt(x ^ 2 + 1) - x", "+oo", "0")]
+        [InlineData("x - sqrt(x ^ 2 - x)", "+oo", "1/2")]
+        [InlineData("sqrt(x + 1) - sqrt(x)", "+oo", "0")]
+        [InlineData("sqrt(x ^ 2 + 4 * x) - x", "+oo", "2")]
+        [InlineData("sqrt(x ^ 2 + x) + x", "-oo", "-1/2")]
+        public void ADifferenceOfRootsGoesThroughItsConjugate(string expression, string destination, string expected) =>
+            AssertLimit(expression, destination, expected);
+
+        // sqrt(x^2 + x) written as x * sqrt(1 + 1/x) says its own growth, and then +oo can
+        // simply be substituted into what is left. An odd degree is not covered: the growth of
+        // sqrt(x^3 + x) is x^(3/2), and a quotient with that in it is one l'Hopital's rule is
+        // stopped from working on, so sqrt(x^3 + x) / x^2 is still left unevaluated.
+        [Theory]
+        [InlineData("sqrt(x ^ 2 + x) / x", "+oo", "1")]
+        [InlineData("x / sqrt(x ^ 2 + 1)", "+oo", "1")]
+        [InlineData("sqrt(4 * x ^ 2 + 1) / x", "+oo", "2")]
+        [InlineData("sqrt(x ^ 2 + 1) / sqrt(x ^ 2 + 3 * x)", "+oo", "1")]
+        [InlineData("sqrt(x ^ 2 - x) / x", "+oo", "1")]
+        [InlineData("sqrt(x ^ 4 + x) / x ^ 2", "+oo", "1")]
+        [InlineData("sqrt(x ^ 6 + 1) / x ^ 3", "+oo", "1")]
+        public void ARootOfAPolynomialSaysHowFastItGrows(string expression, string destination, string expected) =>
+            AssertLimit(expression, destination, expected);
+
+        // The forms that already had an answer must keep it, including the ones where the two
+        // parts tend to the same infinity and so are not indeterminate at all.
+        [Theory]
+        [InlineData("x + sqrt(x)", "+oo", "+oo")]
+        [InlineData("sqrt(x ^ 2 + x) + x", "+oo", "+oo")]
+        [InlineData("x - x", "+oo", "0")]
+        [InlineData("(x ^ 2 + 1) / (x ^ 2 - 1)", "+oo", "1")]
+        [InlineData("1 / x - 1 / x ^ 2", "+oo", "0")]
+        [InlineData("(1 + 1/x) ^ x", "+oo", "e")]
+        public void EstablishedLimitsAreUnaffected(string expression, string destination, string expected) =>
+            AssertLimit(expression, destination, expected);
+
+        /// <summary>
+        /// Neither reading settles this one: the two parts grow alike, so the ratio says
+        /// nothing, and the conjugate is a quotient the solvers cannot read either. Leaving it
+        /// unevaluated is honest; taking a long time over it is not, and the conjugate used to
+        /// be handed to l'Hopital's rule, which spent seconds arriving at the same nothing.
+        /// </summary>
+        [Fact]
+        public void AFormNeitherReadingSettlesTerminatesWithoutClaimingAnAnswer()
+        {
+            var task = Task.Run(() =>
+                "sqrt(x ^ 2 + 3 * x) - sqrt(x ^ 2 + 1)".ToEntity().Limit("x", "+oo".ToEntity()));
+            Assert.True(task.Wait(TimeSpan.FromSeconds(30)), "the limit did not terminate");
+            Assert.IsType<Entity.Limitf>(task.Result);
+        }
+    }
+}

Sources/Tests/UnitTests/Calculus/LimitAtInfinityTest.cs

@@ -66,25 +66,22 @@ public void EstablishedLimitsAreUnaffected(string expression, string destination
             Assert.Equal(expected.ToEntity(), Limit(expression, destination).Simplify());
 
         /// <summary>
-        /// Differentiating both parts of x / sqrt(x^2 + 1) gives its own reciprocal, so the rule
-        /// never settles and has to be stopped; x^20 / e^x would settle but only after more
+        /// Differentiating both parts of (x + sin(x)) / x gives back something no closer to an
+        /// answer, so the rule has to be stopped; x^20 / e^x would settle but only after more
         /// steps than the bound allows. Either way the limit is left unevaluated, which is
         /// honest. What must not happen is a hang, and the answer must not become NaN either,
         /// since NaN asserts that the limit does not exist.
         /// </summary>
         /// <remarks>
-        /// The last two are the forms that showed a bound on the number of steps is not by
-        /// itself a bound on the work. Each step asks what the two parts of its quotient tend
-        /// to, and those are limits the rule may be applied to in turn, so sixteen steps deep
-        /// is far more than sixteen steps: both of these ran for over twenty seconds before the
-        /// rule was stopped from going round a cycle or handing on a quotient bigger than the
-        /// one it came from.
+        /// The last is the form that showed a bound on the number of steps is not by itself a
+        /// bound on the work. Each step asks what the two parts of its quotient tend to, and
+        /// those are limits the rule may be applied to in turn, so sixteen steps deep is far
+        /// more than sixteen steps: it ran for over twenty seconds before the rule was stopped
+        /// from going round a cycle or handing on a quotient bigger than the one it came from.
         /// </remarks>
         [Theory]
-        [InlineData("x / sqrt(x ^ 2 + 1)")]
         [InlineData("(x + sin(x)) / x")]
         [InlineData("x ^ 20 / e ^ x")]
-        [InlineData("sqrt(x ^ 2 - x) / x")]
         [InlineData("x ^ (3/2) * sqrt(1 + 1 / x ^ 2) / x ^ 2")]
         public void AFormTheRuleCannotSettleTerminatesWithoutClaimingAnAnswer(string expression)
         {