Answer the limits that come out as infinity minus infinity

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