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

[Rafael-SOWNet]

Closes #754 — its remaining half. The wrong answer it also reported was fixed in #760.

What was wrong

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, nothing behind it had a reading either, and lim x->+oo ((x!)/x^x)^(1/x) was the last lim:factorial miss in the corpus.

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

Why the expansion goes on the exponent

ln(f!) is 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 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, which is pinned as a test.

The factorial's logarithm is not visible until the base's logarithm 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 #760 are updated, not loosened

#760 pinned (x!)^(1/x), (x!)^(1/ln x) and (x!)^(1/x^2) as left unsettled. Each is answered now. The third is the one #760 recorded as the thing it cost — right by luck rather than by reading, since the same substitution gave 1 for (x!)^(1/x) where the answer is +oo. It comes back as 1 by reading, and BREAKING-CHANGES.md is corrected accordingly.

Measured

Every value checked numerically with mpmath at 50 digits before being claimed, at up to x = 1e9:

x        (x!/x^x)^(1/x)      (x!)^(1/x^2)
1e6      0.36788232          1.000012816
1e9      0.36787945          1.000000020
         1/e = 0.36787944    -> 1

This matters here: SymPy 1.14.0 answers limit((factorial(x)/x**x)**(1/x), x, oo) with 0, which is wrong, so it is not a usable oracle for this one.

Over the 225 generated powers from #760, diffed against master: six results differ, all six are one of the limits above, and every one is an unevaluated node becoming a value. No answer changed into a different answer.

• casbench 112/117 → 113/117, 0 wrong, 0 error, 0 timeout — the lim:factorial problem is solved
• unit tests 5307 passed, 0 failed (5291 on master; 16 added)
• F# wrapper tests 130 passed
• propcheck 0 failures, rootcheck 596/596, simpsweep 0 disagreements of 10463