Stirling.frink

View or download 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]
}


View or download 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 18664 days, 10 hours, 17 minutes ago.