Read a factorial's logarithm outside the power form too (#765)

BREAKING-CHANGES.md

@@ -40,6 +40,7 @@ read first.
 | **silent** | a vanishing sine behind a constant factor | `NaN` | a value — but `sin(x)*ln(x)*2` loses its `0` |
 | **silent** | a limit assembling `oo^0` or `1^oo` | that form's value, `1` | the real answer, or unevaluated |
 | **silent** | a factorial under a vanishing exponent | unevaluated | a value, by Stirling |
+| **silent** | `ln(x!)` in a limit | unevaluated, and `ln(x!)/ln(x)` was `NaN` | a value |
 | **silent** | some integrals | not antiderivatives | closed forms |
 | **silent** | two integrals | answered correctly | unevaluated — a deliberate loss |
 | loud | `Compile` over a missing variable | `KeyNotFoundException` | `UncompilableNodeException` |
@@ -894,6 +895,39 @@ became unevaluated, 1 `NaN` became unevaluated, and 3 right values became uneval
 [#754](https://github.com/asc-community/AngouriMath/issues/754), PR
 [#760](https://github.com/asc-community/AngouriMath/pull/760).
 
+### A factorial's logarithm is read wherever it appears
+
+The expansion below arrived inside the rule that turns `f^g` into `e^(g * ln f)`, so it reached a
+factorial only under a vanishing exponent — the harder question answered and the easier one left:
+
+```
+lim x->+oo ln(x!) / (x * ln(x))   was  unevaluated   is  1
+lim x->+oo ln(x!) / x             was  unevaluated   is  +oo
+lim x->+oo ln(x!) / x^2           was  unevaluated   is  0
+lim x->+oo ln(x!) - x * ln(x)     was  unevaluated   is  -oo
+lim x->+oo ln((2*x)!) / (x*ln(x)) was  unevaluated   is  2
+lim x->+oo ln(x!) / ln(x)         was  NaN           is  +oo
+```
+
+The last is the one that was a **wrong answer** rather than a missing one: `NaN` claims the limit
+does not exist, and that quotient is asymptotic to `x` — 9.6e11 at `x = 1e12`.
+
+Applying the expansion wherever a factorial's logarithm appears is **not** unconditionally sound, so
+it is guarded. What Stirling drops is `1/(12f)`, and what that costs the answer is the rate at which
+the answer moves with the logarithm — so the coefficient is found by putting a variable where the
+logarithm is and differentiating, and the rewrite is refused unless `coefficient / f` vanishes.
+`x * (ln(x!) - (x*ln(x) - x + ln(2*pi*x)/2))` is `1/12`, built out of the dropped term itself, and an
+unguarded rewrite would answer it `0`.
+
+**What it costs.** `lim x->+oo ln((x+1)!) / (x*ln(x))` was unevaluated in about a second and now
+takes longer than 20 s to be unevaluated. The rewrite is right and the expression it hands over is
+right; the machinery cannot take that limit quickly, which is demonstrable without any of this by
+writing the expanded form out by hand. One of 70 in a generated factorial sweep, where 24 already
+timed out before this change.
+
+[#765](https://github.com/asc-community/AngouriMath/issues/765), PR
+[#766](https://github.com/asc-community/AngouriMath/pull/766).
+
 ### A factorial under a vanishing exponent has a limit
 
 A power whose base holds a factorial had none at all. Every route out of `ln(f)` runs through

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

@@ -326,6 +326,67 @@ private static Entity ApplySecondRemarkable(Entity expr, Variable x, Entity dest
             return (power * LogarithmExpanded(@base, x, dest)).InnerSimplified;
         }
 
+        /// <summary>
+        /// Stirling's expansion of <c>ln(f!)</c>, whose error is the vanishing
+        /// <c>1/(12f) + O(1/f^3)</c>.
+        /// </summary>
+        private static Entity StirlingSeries(Entity argument)
+            => argument * MathS.Ln(argument) - argument + MathS.Ln(2 * MathS.pi * argument) / 2;
+
+        /// <summary>
+        /// The expression with the logarithm of every diverging factorial in it replaced by
+        /// Stirling's expansion, or <see langword="null"/> where there is none or the error the
+        /// expansion drops would not vanish out of the answer.
+        /// </summary>
+        /// <remarks>
+        /// <see cref="StirlingExponent"/> can state its guard as <c>power / f -> 0</c> because
+        /// it knows the shape it is working in: the answer there is
+        /// <c>e^(power * ln(base))</c>, so <c>power</c> *is* the coefficient the dropped
+        /// <c>1/(12f)</c> gets multiplied by. Anywhere else that coefficient has to be found
+        /// rather than read off, and it is found by putting a variable where the logarithm is
+        /// and differentiating with respect to it.
+        /// <para/>
+        /// This is not a rewrite that may be applied wherever a factorial's logarithm appears.
+        /// <c>x * (ln(x!) - (x*ln(x) - x + ln(2*pi*x)/2))</c> is <c>1/12</c> -- an expression
+        /// built out of the dropped term itself -- and a rewrite that did not ask would answer
+        /// it <c>0</c>. There the coefficient is <c>x</c> and <c>x / x</c> does not vanish, so
+        /// it is refused. For <c>ln(x!) / x</c> the coefficient is <c>1/x</c> and the ratio is
+        /// <c>1/x^2</c>, so it is allowed.
+        /// <para/>
+        /// The coefficient is put back in terms of the logarithm before it is judged, since an
+        /// expression need not be linear in it: for <c>ln(x!)^2</c> the derivative is
+        /// <c>2*ln(x!)</c>, which grows like <c>2*x*ln(x)</c> and is refused -- correctly, as
+        /// the difference of the squares is <c>ln(x)/6</c> and diverges.
+        /// <a href="https://github.com/asc-community/AngouriMath/issues/765">#765</a>
+        /// </remarks>
+        private static Entity? StirlingRewritten(Entity expr, Variable x, Entity dest)
+        {
+            var logarithms = expr.Nodes
+                .OfType<Logf>()
+                .Where(logarithm => logarithm.Base == MathS.e
+                                    && logarithm.Antilogarithm is Factorialf { Argument: var argument }
+                                    && argument.ContainsNode(x)
+                                    && EvalAssumingContinuous(argument.Limit(x, dest)) == Real.PositiveInfinity)
+                .Distinct()
+                .ToList();
+            if (logarithms.Count == 0)
+                return null;
+            var rewritten = expr;
+            foreach (var logarithm in logarithms)
+            {
+                var argument = ((Factorialf)logarithm.Antilogarithm).Argument;
+                var probe = Variable.CreateTemp(expr.Vars);
+                var coefficient = expr.Substitute(logarithm, probe)
+                    .Differentiate(probe).InnerSimplified.Substitute(probe, logarithm);
+                if (coefficient.Nodes.Any(node => node is Derivativef || node == MathS.NaN))
+                    return null;
+                if (EvalAssumingContinuous((coefficient / argument).Limit(x, dest)) != 0)
+                    return null;
+                rewritten = rewritten.Substitute(logarithm, StirlingSeries(argument));
+            }
+            return rewritten;
+        }
+
         /// <summary>
         /// <c>ln(antilogarithm)</c> taken apart over products, quotients and powers, with
         /// Stirling's expansion written for the logarithm of a diverging factorial.
@@ -335,8 +396,7 @@ private static Entity LogarithmExpanded(Entity antilogarithm, Variable x, Entity
             {
                 Factorialf(var argument)
                     when EvalAssumingContinuous(argument.Limit(x, dest)) == Real.PositiveInfinity
-                        => argument * MathS.Ln(argument) - argument
-                           + MathS.Ln(2 * MathS.pi * argument) / 2,
+                        => StirlingSeries(argument),
                 Mulf(var a, var b) => LogarithmExpanded(a, x, dest) + LogarithmExpanded(b, x, dest),
                 Divf(var a, var b) => LogarithmExpanded(a, x, dest) - LogarithmExpanded(b, x, dest),
                 Powf(var b, var e) when !e.ContainsNode(x) || b.ContainsNode(x)
@@ -992,6 +1052,16 @@ private static bool IsIndeterminateQuotient(Entity numerator, Entity denominator
                     && Settled(extracted.ComputeLimitDivideEtImpera(x, toInfinity, ApproachFrom.Left)) is { } byGrowth)
                     return byGrowth;
 
+                // A factorial's logarithm has no reading anywhere else -- every route out of
+                // ln(f) differentiates f, and a factorial's derivative wants the digamma
+                // function. SolveAsIndeterminatePower reaches the ones sitting under a
+                // vanishing exponent; this reaches the rest, which is why ln(x!)/x had no
+                // answer while ((x!)/x^x)^(1/x) did.
+                // https://github.com/asc-community/AngouriMath/issues/765
+                if (StirlingRewritten(simplified, x, toInfinity) is { } byStirling
+                    && Settled(ComputeLimit(byStirling, x, toInfinity)) is { } byFactorial)
+                    return byFactorial;
+
                 return Settled(SolveAsDifferenceOfInfinities(simplified, x));
             }
             finally { rewriteDepth--; }

Sources/Tests/UnitTests/Calculus/StirlingLogarithmLimitTest.cs

@@ -0,0 +1,122 @@
+//
+// 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>
+    /// Stirling's expansion arrived with <c>SolveAsIndeterminatePower</c>, which is where
+    /// <c>f^g</c> becomes <c>e^(g * ln f)</c> -- so it found a factorial only under a vanishing
+    /// exponent. <c>((x!)/x^x)^(1/x)</c> was answered and <c>ln(x!)/x</c> was not, which is the
+    /// harder question answered and the easier one left.
+    ///
+    /// This applies it wherever a diverging factorial's logarithm appears, from the position
+    /// that is reached only once nothing else has a reading.
+    /// https://github.com/asc-community/AngouriMath/issues/765
+    /// </summary>
+    public sealed class StirlingLogarithmLimitTest
+    {
+        private static Entity LimitOf(string expression) =>
+            expression.ToEntity().Limit("x", Entity.Number.Real.PositiveInfinity).Simplify();
+
+        private static void AssertLimit(string expression, string expected)
+        {
+            var difference = (LimitOf(expression) - expected.ToEntity()).Simplify();
+            while (difference is Entity.Providedf(var inner, _)) difference = inner;
+            Assert.Equal(Entity.Number.Integer.Create(0), difference);
+        }
+
+        private static void AssertEquals(string expression, Entity expected) =>
+            Assert.Equal(expected.Evaled, LimitOf(expression).Evaled);
+
+        /// <summary>
+        /// <c>ln(n!) ~ n*ln(n)</c> is the standard statement of the result, and these are the
+        /// ordinary ways of writing it. Each value was checked numerically at up to
+        /// <c>x = 1e12</c> before being written here -- <c>ln(x!)/(x*ln(x))</c> approaches 1
+        /// from below and slowly, reaching only 0.9638 there.
+        /// </summary>
+        [Theory]
+        [InlineData("ln(x!) / (x * ln(x))", "1")]
+        [InlineData("ln(x!) / x^2", "0")]
+        [InlineData("ln((2*x)!) / (x * ln(x))", "2")]
+        [InlineData("ln((x+1)!) / x^2", "0")]
+        public void TheLogarithmOfAFactorialIsRead(string expression, string expected) =>
+            AssertLimit(expression, expected);
+
+        [Theory]
+        [InlineData("ln(x!) / x")]
+        [InlineData("ln((2*x)!) / x")]
+        [InlineData("ln((x+1)!) / x")]
+        public void TheLogarithmOverAPolynomialDiverges(string expression) =>
+            AssertEquals(expression, Entity.Number.Real.PositiveInfinity);
+
+        /// <summary>
+        /// <c>ln(x!)/ln(x)</c> was <c>NaN</c> -- the claim that the limit does not exist, where
+        /// it is <c>+oo</c>: the quotient is asymptotic to <c>x</c>, and is 9.6e11 at
+        /// <c>x = 1e12</c>. A wrong answer rather than a missing one, and the only one in this
+        /// family.
+        /// </summary>
+        [Theory]
+        [InlineData("ln(x!) / ln(x)")]
+        [InlineData("ln((2*x)!) / ln(x)")]
+        [InlineData("ln((x+1)!) / ln(x)")]
+        public void AQuotientOfLogarithmsIsNoLongerCalledNonExistent(string expression)
+        {
+            var limit = LimitOf(expression);
+            Assert.NotEqual(MathS.NaN, limit.Evaled);
+            Assert.Equal(Entity.Number.Real.PositiveInfinity, limit.Evaled);
+        }
+
+        /// <summary>
+        /// A difference, where the expansion's *additive* error is what makes this sound at
+        /// all: <c>ln(x!) - x*ln(x)</c> is <c>-x + ln(2*pi*x)/2</c> and diverges downwards.
+        /// The asymptotic for <c>x!</c> itself has a merely relative error and says nothing
+        /// here.
+        /// </summary>
+        [Fact]
+        public void ADifferenceAgainstTheLeadingTermDivergesDownwards() =>
+            AssertEquals("ln(x!) - x * ln(x)", Entity.Number.Real.NegativeInfinity);
+
+        /// <summary>
+        /// **The guard.** What Stirling drops is <c>1/(12f)</c>, and what that costs the answer
+        /// is the rate at which the answer moves with the logarithm -- so the coefficient is
+        /// read off by putting a variable where the logarithm is and differentiating, and the
+        /// rewrite is refused unless <c>coefficient / f</c> vanishes. Here the coefficient is
+        /// <c>x</c> and the ratio is 1, so it is refused and the limit is left unsettled.
+        /// <para/>
+        /// The guard is sufficient rather than necessary, and this expression shows the gap:
+        /// its limit is <c>-oo</c>, which the dropped <c>1/12</c> would not have changed. It is
+        /// refused because the same coefficient on
+        /// <c>x * (ln(x!) - (x*ln(x) - x + ln(2*pi*x)/2))</c> gives <c>1/12</c> exactly -- an
+        /// expression built out of the dropped term -- and nothing available here tells the two
+        /// apart.
+        /// </summary>
+        [Fact]
+        public void ACoefficientThatDoesNotVanishAgainstTheFactorialIsRefused()
+        {
+            var limit = "ln(x!) * x - x^2 * ln(x)".ToEntity()
+                .Limit("x", Entity.Number.Real.PositiveInfinity);
+            Assert.True(limit.Evaled is Entity.Limitf,
+                $"the expansion should be refused here, and the limit came back {limit.Evaled}");
+        }
+
+        /// <summary>
+        /// Expressions with no factorial in them, which must be untouched -- the rewrite is
+        /// reached only where a diverging factorial's logarithm is actually present, and only
+        /// once nothing else has answered.
+        /// </summary>
+        [Theory]
+        [InlineData("ln(x) / x", "0")]
+        [InlineData("ln(x^2) / ln(x)", "2")]
+        [InlineData("(ln(x) + 1) / ln(x)", "1")]
+        public void AnExpressionWithoutAFactorialIsUntouched(string expression, string expected) =>
+            AssertLimit(expression, expected);
+    }
+}