/** 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] }