Read a factorial where it has no logarithm at all (#754)

BREAKING-CHANGES.md

@@ -895,6 +895,40 @@ 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 itself is read where it has no logarithm at all
+
+The two readings below both need the factorial to be under a logarithm already — one supplies the
+logarithm, by turning `f^g` into `e^(g * ln f)`, and the other rewrites a logarithm that is written
+down. Neither reaches `x! / x^x`, which has none:
+
+```
+lim x->+oo x! / x^x           was  unevaluated   is  0
+lim x->+oo x^x / x!           was  unevaluated   is  +oo
+lim x->+oo x! / e^x           was  unevaluated   is  +oo
+lim x->+oo ln(x! / x^x) / x   was  a 20 s+ hang  is  -1
+lim x->+oo ln(x! / x^x)       was  unevaluated   is  -oo
+```
+
+The logarithm is supplied instead: for a positive expression `lim H` is `e^(lim ln H)`, and `ln H` is
+where the expansion applies. The guard is the power of the factorial the expression depends on, which
+is exactly the coefficient `ln(f!)` carries in `ln H`, so the dropped `1/(12f)` has to vanish against
+it — `(x!)^x` gives `x`, and `x/(12x)` does not, so it is refused.
+
+Not by substituting `e^(Stirling(f))` for the factorial, which is the obvious move and measured much
+worse: it puts an `e` to a large exponent into the expression, the machinery evaluates that constant
+to a hundred-digit decimal, and everything downstream carries it. `lim x->+oo (x!/e^x)^(1/x)` went
+from half a second to over a minute that way, on an expression the library already answered.
+
+The rewrite that reads a *logarithm* is broadened at the same time, from `ln(f!)` to any logarithm
+holding a diverging factorial, so `ln(x! / x^x)` is taken apart rather than left as one opaque node.
+
+Over 70 generated factorial expressions, against the release before these three changes: **33
+answered → 66 answered, 25 timeouts → 1, and nothing lost**. Values spot-checked numerically at up to
+`x = 1e6`.
+
+[#754](https://github.com/asc-community/AngouriMath/issues/754), PR
+[#767](https://github.com/asc-community/AngouriMath/pull/767).
+
 ### 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

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

@@ -364,29 +364,135 @@ private static Entity StirlingSeries(Entity argument)
             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)
+                                    && DivergingFactorials(logarithm.Antilogarithm, x, dest).Count > 0)
                 .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));
+                // The logarithm may hold several factorials, and what the expansion drops from
+                // it is the sum of their 1/(12f). Each has to vanish against the coefficient.
+                foreach (var factorial in DivergingFactorials(logarithm.Antilogarithm, x, dest))
+                    if (EvalAssumingContinuous((coefficient / factorial.Argument).Limit(x, dest)) != 0)
+                        return null;
+                rewritten = rewritten.Substitute(logarithm,
+                    LogarithmExpanded(logarithm.Antilogarithm, x, dest));
             }
             return rewritten;
         }
 
+        /// <summary>
+        /// Whether an expression is already being read through its own logarithm further up,
+        /// in which case doing it again buys nothing and costs a great deal.
+        /// </summary>
+        [ThreadStatic] private static bool substitutingFactorial;
+
+        /// <summary>
+        /// The fixed power of <paramref name="target"/> that <paramref name="expr"/> depends
+        /// on, or <see langword="null"/> where there is no such power.
+        /// </summary>
+        /// <remarks>
+        /// Once the expression is read through its logarithm, <c>ln(expr)</c> holds
+        /// <c>ln(target)</c> multiplied by exactly this power -- so it is the coefficient the
+        /// dropped <c>1/(12f)</c> arrives with, and the same guard as everywhere else applies
+        /// to it. Only multiplication, division and a constant power keep a factor a factor; a
+        /// sum does not (<c>x! + x</c> is not a fixed power of the factorial at all), and
+        /// neither does an exponent that moves -- for <c>(x!)^x</c> the power is <c>x</c>,
+        /// <c>x/(12x)</c> does not vanish, and it is refused.
+        /// <para/>
+        /// Read off the shape rather than by differentiating and simplifying. That is the same
+        /// answer at a fraction of the cost: the symbolic route needs a full
+        /// <see cref="Entity.Simplify"/> to cancel the target back out -- <c>InnerSimplified</c>
+        /// leaves <c>p * (1/x^x) / (p/x^x)</c> standing -- and it is asked on every expression
+        /// holding a factorial.
+        /// </remarks>
+        private static int? PowerItAppearsTo(Entity expr, Entity target, Variable x)
+        {
+            if (expr == target)
+                return 1;
+            if (!expr.ContainsNode(target))
+                return 0;
+            switch (expr)
+            {
+                case Mulf(var multiplier, var multiplicand):
+                    return PowerItAppearsTo(multiplier, target, x)
+                         + PowerItAppearsTo(multiplicand, target, x);
+                case Divf(var dividend, var divisor):
+                    return PowerItAppearsTo(dividend, target, x)
+                         - PowerItAppearsTo(divisor, target, x);
+                case Powf(var @base, var exponent)
+                    when !exponent.ContainsNode(x) && !exponent.ContainsNode(target)
+                         && exponent.Evaled is Integer power:
+                    return PowerItAppearsTo(@base, target, x) * (int)power;
+                default:
+                    return null;
+            }
+        }
+
+        /// <summary>
+        /// The distinct factorials in an expression whose argument runs off to infinity, which
+        /// are the ones Stirling's expansion is stated for.
+        /// </summary>
+        private static List<Factorialf> DivergingFactorials(Entity expr, Variable x, Entity dest)
+            => expr.Nodes
+                .OfType<Factorialf>()
+                .Where(factorial => factorial.Argument.ContainsNode(x)
+                                    && EvalAssumingContinuous(factorial.Argument.Limit(x, dest))
+                                       == Real.PositiveInfinity)
+                .Distinct()
+                .ToList();
+
+        /// <summary>
+        /// The logarithm of an expression holding a diverging factorial, with Stirling's
+        /// expansion written into it -- so that the limit is <c>e</c> to the limit of this --
+        /// or <see langword="null"/> where there is no such factorial or the expansion would
+        /// not be sound.
+        /// </summary>
+        /// <remarks>
+        /// <see cref="StirlingRewritten"/> reaches a factorial only where a logarithm is
+        /// already written down, and <c>x! / x^x</c> has none, so that shape had no limit at
+        /// all. It is read here by supplying the logarithm: for a positive expression
+        /// <c>lim H</c> is <c>e^(lim ln H)</c>, and <c>ln H</c> is where the expansion applies.
+        /// <para/>
+        /// **Not by substituting <c>e^(Stirling(f))</c> for the factorial**, which is the
+        /// obvious move and was measured to be much worse. It puts an <c>e</c> to a large
+        /// exponent into the expression, the machinery evaluates that constant to a
+        /// hundred-digit decimal, and everything downstream carries it:
+        /// <c>lim x-&gt;+oo (x!/e^x)^(1/x)</c> went from half a second to over a minute, on an
+        /// expression <see cref="StirlingExponent"/> already answers. Going through the
+        /// logarithm keeps the only <c>e</c> in the final answer.
+        /// <para/>
+        /// The guard is <see cref="PowerItAppearsTo"/>: it is the power of the factorial the
+        /// expression depends on, which is exactly the coefficient <c>ln(f!)</c> carries in
+        /// <c>ln H</c>, so the dropped <c>1/(12f)</c> has to vanish against it as everywhere
+        /// else. <c>(x!)^x</c> gives <c>x</c>, and <c>x/(12x)</c> does not vanish.
+        /// </remarks>
+        private static Entity? StirlingByItsOwnLogarithm(Entity expr, Variable x, Entity dest)
+        {
+            var factorials = DivergingFactorials(expr, x, dest);
+            if (factorials.Count == 0)
+                return null;
+            foreach (var factorial in factorials)
+            {
+                // The power the expression depends on the factorial through is the coefficient
+                // the dropped 1/(12f) arrives with, once the whole thing is read through its
+                // logarithm -- so it is the same guard as everywhere else, and it is refused
+                // where it does not vanish.
+                if (PowerItAppearsTo(expr, factorial, x) is not { } power || power == 0)
+                    return null;
+                if (EvalAssumingContinuous(((Entity)power / factorial.Argument).Limit(x, dest)) != 0)
+                    return null;
+            }
+            return LogarithmExpanded(expr, x, dest);
+        }
+
         /// <summary>
         /// <c>ln(antilogarithm)</c> taken apart over products, quotients and powers, with
         /// Stirling's expansion written for the logarithm of a diverging factorial.
@@ -1062,6 +1168,23 @@ private static bool IsIndeterminateQuotient(Entity numerator, Entity denominator
                     && Settled(ComputeLimit(byStirling, x, toInfinity)) is { } byFactorial)
                     return byFactorial;
 
+                // And the factorials that are not under a logarithm at all, where what may be
+                // dropped is a relative error rather than an additive one. Second, because the
+                // logarithm above is the cheaper reading and the sounder one -- an additive
+                // error that vanishes needs less of the expression to be true of it than a
+                // relative one that tends to 1.
+                if (!substitutingFactorial
+                    && StirlingByItsOwnLogarithm(simplified, x, toInfinity) is { } logarithm)
+                {
+                    substitutingFactorial = true;
+                    try
+                    {
+                        if (Settled(ComputeLimit(logarithm, x, toInfinity)) is { } exponent)
+                            return MathS.e.Pow(exponent).InnerSimplified;
+                    }
+                    finally { substitutingFactorial = false; }
+                }
+
                 return Settled(SolveAsDifferenceOfInfinities(simplified, x));
             }
             finally { rewriteDepth--; }

Sources/Tests/UnitTests/Calculus/StirlingFactorialItselfTest.cs

@@ -0,0 +1,112 @@
+//
+// 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>
+    /// The two earlier readings of Stirling's expansion both need the factorial to be under a
+    /// logarithm already -- one supplies the logarithm itself, by turning <c>f^g</c> into
+    /// <c>e^(g * ln f)</c>, and the other rewrites a logarithm that is written down. Neither
+    /// reaches <c>x! / x^x</c>, which has no logarithm anywhere in it.
+    ///
+    /// Here the logarithm is supplied instead: for a positive expression <c>lim H</c> is
+    /// <c>e^(lim ln H)</c>, and <c>ln H</c> is where the expansion applies. Not by substituting
+    /// <c>e^(Stirling(f))</c> for the factorial, which is the obvious move and measured much
+    /// worse -- it puts an <c>e</c> to a large exponent into the expression, the machinery
+    /// evaluates that constant to a hundred-digit decimal, and <c>(x!/e^x)^(1/x)</c> went from
+    /// half a second to over a minute.
+    ///
+    /// The guard is the power of the factorial the expression depends on, which is exactly the
+    /// coefficient <c>ln(f!)</c> carries in <c>ln H</c>, so the dropped <c>1/(12f)</c> has to
+    /// vanish against it. For <c>(x!)^x</c> that power is <c>x</c>, and <c>x/(12x)</c> does
+    /// not vanish.
+    ///
+    /// https://github.com/asc-community/AngouriMath/issues/754
+    /// </summary>
+    public sealed class StirlingFactorialItselfTest
+    {
+        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>
+        /// The factorial compared against the things it is usually compared against. Each was
+        /// checked numerically before being written here: <c>x!/x^x</c> is 4.0e-433 at
+        /// <c>x = 1000</c>, and <c>x!/e^x</c> is 2.0e+2133.
+        /// </summary>
+        [Fact]
+        public void AFactorialAgainstItsOwnPowerVanishes() =>
+            AssertLimit("x! / x^x", "0");
+
+        [Theory]
+        [InlineData("x^x / x!")]
+        [InlineData("x! / e^x")]
+        [InlineData("x! / (x-1)!")]
+        public void AFactorialAgainstAWeakerGrowthDiverges(string expression) =>
+            AssertEquals(expression, Entity.Number.Real.PositiveInfinity);
+
+        /// <summary>
+        /// The logarithm of a quotient holding a factorial, which the rewrite that reads a
+        /// logarithm now takes apart rather than leaving as one opaque node. Substituting for
+        /// the factorial alone would leave <c>ln(e^S / x^x)</c>, which nothing here simplifies
+        /// -- so these have to go through the logarithm's expansion and not this one.
+        /// </summary>
+        [Theory]
+        [InlineData("ln(x! / x^x) / x", "-1")]
+        [InlineData("ln(x! / x^x) / (x * ln(x))", "0")]
+        [InlineData("ln(x^x / x!) / x", "1")]
+        public void ALogarithmOfAQuotientHoldingAFactorialIsTakenApart(string expression, string expected) =>
+            AssertLimit(expression, expected);
+
+        [Theory]
+        [InlineData("ln(x! / x^x)")]
+        [InlineData("x * ln(x! / x^x)")]
+        [InlineData("ln(x! / x^x) / ln(x)")]
+        public void ALogarithmOfAVanishingQuotientDivergesDownwards(string expression) =>
+            AssertEquals(expression, Entity.Number.Real.NegativeInfinity);
+
+        /// <summary>
+        /// **The guard.** <c>(x!)^x</c> depends on its factorial through the power <c>x</c>,
+        /// so the dropped <c>1/(12f)</c> would arrive as <c>x/(12x)</c>, which does not
+        /// vanish. It is refused there, and the answer is still <c>+oo</c>, reached without
+        /// any of this.
+        /// </summary>
+        [Theory]
+        [InlineData("(x!) ^ x")]
+        [InlineData("(x!) ^ (x^2)")]
+        public void APowerThatGrowsAgainstTheFactorialIsRefused(string expression) =>
+            AssertEquals(expression, Entity.Number.Real.PositiveInfinity);
+
+        /// <summary>
+        /// Nothing without a factorial in it is touched -- the rewrite is reached only where a
+        /// diverging factorial is present, and only once nothing else has answered.
+        /// </summary>
+        [Theory]
+        [InlineData("x^x / e^x")]
+        [InlineData("e^x / x^x")]
+        [InlineData("x^20 / e^x")]
+        public void AnExpressionWithoutAFactorialIsUntouched(string expression)
+        {
+            var limit = LimitOf(expression).Evaled;
+            Assert.True(limit == Entity.Number.Real.PositiveInfinity || limit == Entity.Number.Integer.Create(0),
+                $"{expression} came back {limit}");
+        }
+    }
+}