Download or view Stirling.frink in plain text format
/** 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]
}
Download or view Stirling.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19963 days, 19 hours, 51 minutes ago.