# Stirling.frink

``` /** This program explores Stirling's approximation to the factorial of n, also     written "n!", including to arbitrary precision.  This also includes Bill     Gosper's approximation which is probably better for small numbers.     See:      https://mathworld.wolfram.com/StirlingsApproximation.html     https://en.wikipedia.org/wiki/Stirling's_approximation      Nemes,Gergő, On the coefficients of the asymptotic expansion of n!, Journal of Integer Sequences 13 (2010), no. 6, Article 10.6.6, 5 pp.      http://www.cs.uwaterloo.ca/journals/JIS/VOL13/Nemes/nemes2.pdf */ use pi2.frink use e.frink use ArbitraryPrecision.frink /** Calculates a rough approximation of n! using Stirling's approximation,     but using arbitrary-precision functions so you know that the error is in     Stirling's approximation, and not in the numerics.     Stirling's approximation can be improved by adding additional terms.     See:      Nemes, Gergő, On the coefficients of the asymptotic expansion of n!, Journal of Integer Sequences 13 (2010), no. 6, Article 10.6.6, 5 pp.      http://www.cs.uwaterloo.ca/journals/JIS/VOL13/Nemes/nemes2.pdf     The coefficents in the Nemes paper are somewhat complicated to derive, but     the first 4 terms usually suffice for improved accuracy when n is large.     A program that derives the coefficients is found in     StirlingCoefficients.frink */ factorialStirling[n, digits, additionalTerms=true] := {    prec = getPrecision[]    try    {       setPrecision[digits]       pi2 = Pi.get2Pi[digits]              if additionalTerms == true          prod = 1 + 1/(12 n) + 1/(288 n^2) - 139/(51840 n^3) - 571/(2488320 n^4)       else          prod = 1              return sqrt[pi2, digits] arbitraryPow[n, n+1/2, digits] arbitraryExp[-n] prod    }    finally    {       setPrecision[prec]    } } /** Calculates interval bounds of n! using Stirling's approximation,     using arbitrary-precision functions so you know that the error is in     Stirling's approximation, and not in the numerics.  This version gives you     an interval that provides the upper and lower bounds of the result.     TODO:  Increase the bounds with additional terms? */ factorialStirlingBounds[n, digits] := {    pi2 = Pi.get2Pi[digits]    pow = sqrt[pi2, digits] arbitraryPow[n, n+1/2, digits]    lower = arbitraryExp[-n + 1/(12 n+1)]    upper = arbitraryExp[-n + 1/(12 n)]    return new interval[pow lower, pow upper] } /** Calculates Gosper's formula for n! which works better for small numbers     especially. */ factorialGosper[n, digits] := {    pi = Pi.getPi[digits]    return sqrt[(2n + 1/3) pi, digits] arbitraryPow[n,n] arbitraryExp[-n] } ```