From c64cfe61a39e4964a65cd1414c1cfa3b0801ff3c Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 6 Aug 2026 13:33:11 +0000 Subject: [PATCH] Read a factorial's logarithm by Stirling's expansion (#754) A power whose base holds a factorial had no limit at all. f^g is e^(g * ln f), and SolveAsIndeterminatePower computes the limit of that exponent -- but every route out of ln(f) runs through differentiating f, and a factorial's derivative wants the digamma function, which this library does not have. The rule declined and nothing behind it had a reading either, so lim x->+oo ((x!)/x^x)^(1/x) was the last lim:factorial miss in the corpus. Stirling's expansion is stated for exactly that logarithm: ln(f!) = f*ln(f) - f + ln(2*pi*f)/2 + 1/(12f) + O(1/f^3). Applied to the exponent rather than substituted for the factorial in the base, and that is the whole of what makes it sound -- what is dropped here *vanishes*, where the asymptotic for f! itself has an error that is merely relative and survives being raised to a power. Vanishing is still not sufficient, because the dropped term is multiplied by the exponent the rewrite sits under: an error of 1/(12f) in the logarithm contributes power/(12f) to the exponent, so power/f -> 0 is required. For ((x!)/x^x)^(1/x) that ratio is 1/x^2. (x!)^x fails it and is left alone. The factorial's own logarithm is not visible until the logarithm of the base is taken apart -- ln(x!/x^x) is one node and nothing simplifies it -- so ln is split over products, quotients and powers, confined to logarithms that actually hold a diverging factorial. That split assumes the parts are positive on the approach, which is the same assumption the simplifier's ln(a) + ln(b) = ln(a*b) already makes, and it is reached only by expressions that have no answer at all without it. Three tests from PR #760 pinned these as unsettled and are updated rather than loosened: each is now answered, and (x!)^(1/x^2) -- recorded there as the one thing that change cost, right by luck rather than by reading -- comes back as 1 by reading. Every value checked numerically at up to x = 1e9 before being claimed, since SymPy 1.14.0 answers (x!/x^x)^(1/x) with 0 and is not a usable oracle here. Measured over 225 generated powers: six results differ from master and all six are an unevaluated node becoming a value, with no answer changed. casbench 112/117 -> 113/117 with 0 wrong. Unit tests 5291 -> 5307 passing with 16 added, F# 130, propcheck 0 failures, rootcheck 596/596, simpsweep 0 disagreements. Co-Authored-By: Claude Opus 5 --- BREAKING-CHANGES.md | 46 ++++++- .../Continuous/Limits/Transformations.cs | 99 +++++++++++++-- .../IndeterminatePowerSubstitutionTest.cs | 31 +++-- .../Calculus/StirlingFactorialLimitTest.cs | 113 ++++++++++++++++++ 4 files changed, 254 insertions(+), 35 deletions(-) create mode 100644 Sources/Tests/UnitTests/Calculus/StirlingFactorialLimitTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index d3c7c4a8c..19bba9395 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -39,6 +39,7 @@ read first. | **silent** | many limits | `NaN`, or unevaluated | a value | | **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** | some integrals | not antiderivatives | closed forms | | **silent** | two integrals | answered correctly | unevaluated — a deliberate loss | | loud | `Compile` over a missing variable | `KeyNotFoundException` | `UncompilableNodeException` | @@ -871,18 +872,18 @@ answer off one. `0^0` is deliberately untouched for the opposite reason: its `Na says `lim x->0 x^x` does not exist, which `LimitTest.TestNoLimit` pins, and declining it would turn a considered "does not exist" into "not settled". -**What it costs.** Three limits that were answered `1` are now unevaluated, because the substitution +**What it costs.** Two limits that were answered `1` are now unevaluated, because the substitution that used to land on them is gone and nothing else has a reading: ``` -lim x->+oo (x!) ^ (1/x^2) was 1 is unevaluated (it is 1) lim x->0 (1/x) ^ x was 1 is unevaluated lim x->0 (1/x) ^ (x^2) was 1 is unevaluated ``` -The first was right by luck: the same substitution gave `1` for `(x!)^(1/x)`, where the answer is -`+oo`, and the two cannot be told apart at the point where the value is read off. Answering it -properly wants Stirling's expansion of `ln(x!)`, which is the second half of #754. The other two are +A third, `lim x->+oo (x!)^(1/x^2)`, was withdrawn here and is answered `1` again by the entry +below — it was right by luck when this landed, since the same substitution gave `1` for +`(x!)^(1/x)` where the answer is `+oo`, and the two cannot be told apart at the point where the +value is read off. Stirling's expansion answers it by reading instead. The other two are two-sided limits at `0` whose base has no two-sided limit at all and which are not real to the left of `0`, so the `1` came from the complex continuation; **their one-sided readings are unchanged** — `lim x->0+ (1/x)^x` is still `1`. @@ -893,6 +894,41 @@ 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 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 +differentiating `f`, and a factorial's derivative wants the digamma function, which this library does +not have — so the rule that reads `f^g` as `e^(g * ln f)` declined and nothing behind it had a +reading either: + +``` +lim x->+oo (x! / x^x) ^ (1/x) was unevaluated is 1/e +lim x->+oo (x!) ^ (1/x) was unevaluated is +oo +lim x->+oo (x!) ^ (1/ln(x)) was unevaluated is +oo +lim x->+oo (x!) ^ (1/x^2) was unevaluated is 1 +lim x->+oo (x! / e^x) ^ (1/x) was unevaluated is +oo +lim x->+oo ((x+1)! / x!) ^ (1/x) was unevaluated is 1 +``` + +Stirling's expansion is stated for exactly that logarithm — `ln(f!)` is +`f*ln(f) - f + ln(2*pi*f)/2 + 1/(12f) + O(1/f^3)` — and it is applied to the **exponent** rather than +substituted for the factorial in the base. That is what makes it sound: what is dropped here +*vanishes*, where the asymptotic for `f!` itself has an error that is merely relative and survives +being raised to a power. Vanishing is still not enough on its own, since the dropped term is +multiplied by the exponent it sits under, so `power / f -> 0` is required; `(x!)^x` fails that and is +left alone. + +Every one of these values was checked numerically before being claimed — `(x!/x^x)^(1/x)` is +`0.3678794453` at `x = 1e9` against `1/e = 0.3678794412`. Note that **SymPy 1.14.0 answers that one +`0`**, so it is not a usable oracle here. + +Nothing without a factorial in it changes: over 225 generated powers, six results differ and all six +are one of these, every one of them an unevaluated node becoming a value. `casbench` 112/117 → +113/117. + +[#754](https://github.com/asc-community/AngouriMath/issues/754), PR +[#764](https://github.com/asc-community/AngouriMath/pull/764). + ### `lim x->2 signum(x)` no longer kills the process `Signumf` was the one node whose limit override handed back an unevaluated limit of the very diff --git a/Sources/AngouriMath/Functions/Continuous/Limits/Transformations.cs b/Sources/AngouriMath/Functions/Continuous/Limits/Transformations.cs index d12762fe3..ef37c3983 100644 --- a/Sources/AngouriMath/Functions/Continuous/Limits/Transformations.cs +++ b/Sources/AngouriMath/Functions/Continuous/Limits/Transformations.cs @@ -247,23 +247,37 @@ private static Entity ApplySecondRemarkable(Entity expr, Variable x, Entity dest return null; if (EvalAssumingContinuous(power.Limit(x, dest, side)) != 0) return null; - var baseLimit = EvalAssumingContinuous(@base.Limit(x, dest, side)); - if (baseLimit != 0 && !IsInfiniteNode(baseLimit)) - return null; - // Every route out of ln(f) runs through differentiating f, so a base this library - // cannot differentiate is one the rewrite cannot finish on: it would only hand the - // rules an expression with a hole in it and let them work at it. A factorial is the - // case that matters -- its derivative wants the digamma function, which is not here, - // and comes back as NaN -- and lim x->+oo ((x!) / x^x)^(1/x) is the expression. It - // has no answer either way, and without this it takes a long time not to find one. - var derivative = @base.Differentiate(x).InnerSimplified; - if (derivative.Nodes.Any(node => node is Derivativef || node == MathS.NaN)) - return null; + + // f^g is e^(g * ln f), and this rule computes the limit of the exponent. Where the + // base holds a diverging factorial, that logarithm is what Stirling's expansion is + // stated for -- so the expansion is applied here, to the exponent, rather than to + // the base, where it would have to reproduce the factorial itself and its merely + // *relative* error. https://github.com/asc-community/AngouriMath/issues/754 + var byStirling = StirlingExponent(@base, power, x, dest); + + if (byStirling is null) + { + var baseLimit = EvalAssumingContinuous(@base.Limit(x, dest, side)); + if (baseLimit != 0 && !IsInfiniteNode(baseLimit)) + return null; + // Every route out of ln(f) runs through differentiating f, so a base this + // library cannot differentiate is one the rewrite cannot finish on: it would + // only hand the rules an expression with a hole in it and let them work at it. + // A factorial is the case that matters -- its derivative wants the digamma + // function, which is not here, and comes back as NaN. That is why the expansion + // above is tried first: where it applies there is no factorial left to + // differentiate, and this guard is asking about an expression the rule is no + // longer going to use. + var derivative = @base.Differentiate(x).InnerSimplified; + if (derivative.Nodes.Any(node => node is Derivativef || node == MathS.NaN)) + return null; + } indeterminatePowerDepth++; try { - if (ComputeLimit((power * MathS.Ln(@base)).InnerSimplified, x, dest, side) is not { } exponent + var exponentExpr = byStirling ?? (power * MathS.Ln(@base)).InnerSimplified; + if (ComputeLimit(exponentExpr, x, dest, side) is not { } exponent || exponent.Evaled == MathS.NaN) return null; return MathS.e.Pow(exponent).InnerSimplified; @@ -271,6 +285,65 @@ private static Entity ApplySecondRemarkable(Entity expr, Variable x, Entity dest finally { indeterminatePowerDepth--; } } + /// + /// power * ln(base) with the logarithm of every diverging factorial in it + /// replaced by Stirling's expansion, or where there is no such + /// factorial or the expansion would not be sound here. + /// + /// + /// ln(f!) is f*ln(f) - f + ln(2*pi*f)/2 + 1/(12f) + O(1/f^3), and what is + /// dropped **vanishes** -- where the asymptotic for f! itself has an error that + /// is merely relative. That is why the expansion is written for the logarithm and + /// applied to this exponent rather than substituted for the factorial in the base. + /// + /// Vanishing is still not sufficient, because the dropped term is multiplied by the + /// exponent the rewrite sits under: the answer is e^(power * ln(base)), so an + /// error of 1/(12f) in the logarithm contributes power/(12f) to the + /// exponent. Requiring power / f -> 0 is what makes it disappear. For + /// ((x!) / x^x)^(1/x) that ratio is 1/x^2. + /// + /// The logarithm has to be taken apart before the factorial's own is visible: + /// ln(x!/x^x) is one node, and nothing here simplifies it. Splitting it over + /// products, quotients and powers assumes the parts are positive on the approach, which + /// is the same assumption the simplifier's ln(a) + ln(b) = ln(a*b) already + /// makes; it is confined to logarithms that actually hold a diverging factorial, so it + /// is reached only by expressions that have no answer at all without it. + /// #754 + /// + private static Entity? StirlingExponent(Entity @base, Entity power, Variable x, Entity dest) + { + var factorials = @base.Nodes.OfType() + .Where(f => f.Argument.ContainsNode(x) + && EvalAssumingContinuous(f.Argument.Limit(x, dest)) == Real.PositiveInfinity) + .ToList(); + if (factorials.Count == 0) + return null; + // The dropped 1/(12f) is multiplied by the exponent this sits under, so it only + // disappears where power/f does. + foreach (var factorial in factorials) + if (EvalAssumingContinuous((power / factorial.Argument).Limit(x, dest)) != 0) + return null; + return (power * LogarithmExpanded(@base, x, dest)).InnerSimplified; + } + + /// + /// ln(antilogarithm) taken apart over products, quotients and powers, with + /// Stirling's expansion written for the logarithm of a diverging factorial. + /// + private static Entity LogarithmExpanded(Entity antilogarithm, Variable x, Entity dest) + => antilogarithm switch + { + Factorialf(var argument) + when EvalAssumingContinuous(argument.Limit(x, dest)) == Real.PositiveInfinity + => argument * MathS.Ln(argument) - argument + + MathS.Ln(2 * MathS.pi * argument) / 2, + 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) + => e * LogarithmExpanded(b, x, dest), + _ => MathS.Ln(antilogarithm) + }; + /// /// How deep the rewriting of one power into another may go. The exponent it asks about /// is a limit in its own right and may hold a power of the same shape, so without a diff --git a/Sources/Tests/UnitTests/Calculus/IndeterminatePowerSubstitutionTest.cs b/Sources/Tests/UnitTests/Calculus/IndeterminatePowerSubstitutionTest.cs index 8042facb2..8abd5d1a9 100644 --- a/Sources/Tests/UnitTests/Calculus/IndeterminatePowerSubstitutionTest.cs +++ b/Sources/Tests/UnitTests/Calculus/IndeterminatePowerSubstitutionTest.cs @@ -94,31 +94,28 @@ public void ALogarithmicExponentOverAnExponentialDiverges() => AssertDiverges("(e^x) ^ (1/ln(x))", "+oo"); /// - /// The factorial cases from the report. A factorial's derivative wants the digamma - /// function, which this library does not have, so SolveAsIndeterminatePower - /// declines and nothing else has a reading either -- they are left unsettled rather - /// than answered wrongly. (x!)^(1/x) is +oo: by Stirling it grows like - /// x/e, and (100!)^(1/100) is already 37.99. + /// The factorial cases from the report, which were left unsettled when this guard + /// landed and are answered now that Stirling's expansion of ln(f!) reaches them. + /// (x!)^(1/x) grows like x/e -- (100!)^(1/100) is already 37.99 -- + /// and it is the value the substitution used to read off as 1. /// [Theory] [InlineData("(x!) ^ (1/x)")] [InlineData("(x!) ^ (1/ln(x))")] - public void AFactorialUnderAVanishingExponentIsLeftUnsettled(string expression) => - AssertNotSettled(expression, "+oo"); + public void AFactorialUnderAVanishingExponentDiverges(string expression) => + AssertDiverges(expression, "+oo"); /// - /// **What this costs.** (x!)^(1/x^2) is 1 -- its logarithm is - /// ln(x!)/x^2 ~ ln(x)/x, which vanishes -- and the old substitution happened to - /// land on that 1 for the same reason it landed on 1 for (x!)^(1/x), where the - /// answer is +oo. It was right by luck rather than by reading, and the two - /// cannot be told apart at the point where the value is read off. - /// - /// Answering it properly wants Stirling's expansion of ln(x!), which is the - /// second half of #754 and is not attempted here. + /// **What this used to cost, and no longer does.** (x!)^(1/x^2) is 1 -- its + /// logarithm is ln(x!)/x^2 ~ ln(x)/x, which vanishes -- and the old substitution + /// landed on that 1 for the same reason it landed on 1 for (x!)^(1/x), where the + /// answer is +oo. It was right by luck rather than by reading, so it was + /// withdrawn along with the case that was wrong. Stirling's expansion answers it by + /// reading, which is what makes keeping it worth anything. /// [Fact] - public void TheFactorialCaseThatWasAccidentallyRightIsAlsoLeftUnsettled() => - AssertNotSettled("(x!) ^ (1/x^2)", "+oo"); + public void TheFactorialCaseThatWasAccidentallyRightIsAnsweredByReading() => + AssertLimit("(x!) ^ (1/x^2)", "+oo", "1"); /// /// 0^0 is deliberately **not** guarded. It evaluates to NaN, and that NaN is diff --git a/Sources/Tests/UnitTests/Calculus/StirlingFactorialLimitTest.cs b/Sources/Tests/UnitTests/Calculus/StirlingFactorialLimitTest.cs new file mode 100644 index 000000000..c3589af76 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/StirlingFactorialLimitTest.cs @@ -0,0 +1,113 @@ +// +// 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 +{ + /// + /// A power whose base holds a factorial had no limit at all, because every route out of + /// ln(f) runs through differentiating f and a factorial's derivative wants the + /// digamma function, which this library does not have. + /// + /// f^g is e^(g * ln f), and Stirling's expansion is stated for exactly that + /// logarithm: ln(f!) = f*ln(f) - f + ln(2*pi*f)/2 + 1/(12f) + O(1/f^3). Applying it + /// to the exponent rather than substituting for the factorial in the base is what makes it + /// sound -- what is dropped here **vanishes**, where the asymptotic for f! itself has + /// an error that is merely relative and survives being raised to a power. + /// + /// Vanishing is still not enough on its own: the dropped term is multiplied by the exponent + /// the rewrite sits under, so power / f -> 0 is required. For (x!/x^x)^(1/x) + /// that ratio is 1/x^2. + /// + /// https://github.com/asc-community/AngouriMath/issues/754 + /// + public sealed class StirlingFactorialLimitTest + { + 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 AssertDiverges(string expression) => + Assert.Equal(Entity.Number.Real.PositiveInfinity, LimitOf(expression).Evaled); + + /// + /// The reported case, and the last lim:factorial miss in the corpus. + /// x!/x^x is sqrt(2*pi*x) * e^(-x) to leading order, so its x-th root is + /// 1/e -- checked numerically at x = 1e9, where the quotient's root is + /// 0.3678794453 against 1/e = 0.3678794412. + /// + [Fact] + public void TheCorpusMissIsAnswered() => + AssertLimit("(x! / x^x) ^ (1/x)", "1/e"); + + /// + /// The wrong answer the same issue reported, which was 1 and then unevaluated, + /// and is a value now. (x!)^(1/x) is asymptotic to x/e: at x = 1e9 + /// it is 367879445, against x/e = 367879441. + /// + [Theory] + [InlineData("(x!) ^ (1/x)")] + [InlineData("(x!) ^ (2/x)")] + [InlineData("(x!) ^ (1/ln(x))")] + [InlineData("(x! * x) ^ (1/x)")] + [InlineData("(x! / e^x) ^ (1/x)")] + public void AFactorialUnderAVanishingExponentDiverges(string expression) => + AssertDiverges(expression); + + /// + /// Where the exponent vanishes fast enough to hold the factorial's growth down. Each of + /// these was checked numerically before being written here. + /// + [Theory] + [InlineData("(x!) ^ (1/x^2)", "1")] + [InlineData("(x! / x^x) ^ (1/x^2)", "1")] + [InlineData("((x+1)! / x!) ^ (1/x)", "1")] + public void AnExponentThatVanishesFasterHoldsItDown(string expression, string expected) => + AssertLimit(expression, expected); + + /// + /// The guard, and the reason the expansion is not simply applied wherever a factorial + /// appears. What Stirling drops is 1/(12f), and the answer is + /// e^(power * ln(base)) -- so that error contributes power/(12f) to the + /// exponent and only disappears where power/f does. (x!)^x fails that + /// (the ratio is x/x, which is 1), and it is left unanswered rather than + /// answered from an expansion that does not hold there. + /// + [Theory] + [InlineData("(x!) ^ x")] + [InlineData("(x!) ^ (x^2)")] + public void AnExponentThatDoesNotVanishAgainstTheFactorialIsNotRewritten(string expression) + { + var limit = expression.ToEntity().Limit("x", Entity.Number.Real.PositiveInfinity); + Assert.True(limit.Evaled is Entity.Limitf or Entity.Number.Real { IsFinite: false }, + $"{expression} should be left unsettled or diverge, and came back {limit.Evaled}"); + } + + /// + /// The powers with no factorial in them, which must reach the same answers by the same + /// route they always did -- the expansion is only reached where a diverging factorial + /// is actually present. + /// + [Theory] + [InlineData("x ^ (1/x)", "1")] + [InlineData("(1 + 1/x) ^ x", "e")] + [InlineData("(x - 5) ^ x / x ^ x", "e^(-5)")] + [InlineData("(x^2) ^ (1/ln(x))", "e^2")] + [InlineData("(1 - 1/x^2) ^ (x^2)", "1/e")] + public void APowerWithoutAFactorialIsUntouched(string expression, string expected) => + AssertLimit(expression, expected); + } +}