Skip to content

Read a factorial's logarithm by Stirling's expansion (#754) #764

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Rafael-SOWNet merged 1 commit into from
Aug 6, 2026
Merged

Read a factorial's logarithm by Stirling's expansion (#754) #764

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 41 additions & 5 deletions BREAKING-CHANGES.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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` |
Expand Down Expand Up @@ -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`.
Expand All @@ -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
Expand Down
99 changes: 86 additions & 13 deletions Sources/AngouriMath/Functions/Continuous/Limits/Transformations.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -247,30 +247,103 @@ 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;
}
finally { indeterminatePowerDepth--; }
}

/// <summary>
/// <c>power * ln(base)</c> with the logarithm of every diverging factorial in it
/// replaced by Stirling's expansion, or <see langword="null"/> where there is no such
/// factorial or the expansion would not be sound here.
/// </summary>
/// <remarks>
/// <c>ln(f!)</c> is <c>f*ln(f) - f + ln(2*pi*f)/2 + 1/(12f) + O(1/f^3)</c>, and what is
/// dropped **vanishes** -- where the asymptotic for <c>f!</c> 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.
/// <para/>
/// Vanishing is still not sufficient, because the dropped term is multiplied by the
/// exponent the rewrite sits under: the answer is <c>e^(power * ln(base))</c>, so an
/// error of <c>1/(12f)</c> in the logarithm contributes <c>power/(12f)</c> to the
/// exponent. Requiring <c>power / f -> 0</c> is what makes it disappear. For
/// <c>((x!) / x^x)^(1/x)</c> that ratio is <c>1/x^2</c>.
/// <para/>
/// The logarithm has to be taken apart before the factorial's own is visible:
/// <c>ln(x!/x^x)</c> 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 <c>ln(a) + ln(b) = ln(a*b)</c> 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.
/// <a href="https://github.com/asc-community/AngouriMath/issues/754">#754</a>
/// </remarks>
private static Entity? StirlingExponent(Entity @base, Entity power, Variable x, Entity dest)
{
var factorials = @base.Nodes.OfType<Factorialf>()
.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;
}

/// <summary>
/// <c>ln(antilogarithm)</c> taken apart over products, quotients and powers, with
/// Stirling's expansion written for the logarithm of a diverging factorial.
/// </summary>
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)
};

/// <summary>
/// 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
Expand Down
31 changes: 14 additions & 17 deletions Sources/Tests/UnitTests/Calculus/IndeterminatePowerSubstitutionTest.cs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,31 +94,28 @@ public void ALogarithmicExponentOverAnExponentialDiverges() =>
AssertDiverges("(e^x) ^ (1/ln(x))", "+oo");

/// <summary>
/// The factorial cases from the report. A factorial's derivative wants the digamma
/// function, which this library does not have, so <c>SolveAsIndeterminatePower</c>
/// declines and nothing else has a reading either -- they are left unsettled rather
/// than answered wrongly. <c>(x!)^(1/x)</c> is <c>+oo</c>: by Stirling it grows like
/// <c>x/e</c>, and <c>(100!)^(1/100)</c> 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 <c>ln(f!)</c> reaches them.
/// <c>(x!)^(1/x)</c> grows like <c>x/e</c> -- <c>(100!)^(1/100)</c> is already 37.99 --
/// and it is the value the substitution used to read off as <c>1</c>.
/// </summary>
[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");

/// <summary>
/// **What this costs.** <c>(x!)^(1/x^2)</c> is 1 -- its logarithm is
/// <c>ln(x!)/x^2 ~ ln(x)/x</c>, which vanishes -- and the old substitution happened to
/// land on that 1 for the same reason it landed on 1 for <c>(x!)^(1/x)</c>, where the
/// answer is <c>+oo</c>. 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.
/// <para/>
/// Answering it properly wants Stirling's expansion of <c>ln(x!)</c>, which is the
/// second half of #754 and is not attempted here.
/// **What this used to cost, and no longer does.** <c>(x!)^(1/x^2)</c> is 1 -- its
/// logarithm is <c>ln(x!)/x^2 ~ ln(x)/x</c>, which vanishes -- and the old substitution
/// landed on that 1 for the same reason it landed on 1 for <c>(x!)^(1/x)</c>, where the
/// answer is <c>+oo</c>. 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.
/// </summary>
[Fact]
public void TheFactorialCaseThatWasAccidentallyRightIsAlsoLeftUnsettled() =>
AssertNotSettled("(x!) ^ (1/x^2)", "+oo");
public void TheFactorialCaseThatWasAccidentallyRightIsAnsweredByReading() =>
AssertLimit("(x!) ^ (1/x^2)", "+oo", "1");

/// <summary>
/// <c>0^0</c> is deliberately **not** guarded. It evaluates to NaN, and that NaN is
Expand Down
Loading
Loading