/** This reproduce the additional Stirling coefficients for Stirling's
    approximation to the factorial funcion from:

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

    This requires use of the generalized binomial theorem with rational
    numbers (instead of integers) as arguments.  This can be obtained using
    the result that binomial[r, k] = fallingFactorial[r, k] / k!

    Note that, from
    https://en.wikipedia.org/wiki/Stirling%27s_approximation

    "As n → ∞, the error in the truncated series is asymptotically equal to the
     first omitted term. This is an example of an asymptotic expansion. It is
     not a convergent series; for any particular value of n there are only so
     many terms of the series that improve accuracy, after which accuracy
     worsens."
*/

// Used for calculating falling factorial.
use Pochhammer.frink

/** This calculates StirlingCoeffient a_k.  This will be the term in the sum
    a_k / n^k

    For example, the coefficients are
    (1 + 1/(12 n) + 1/(288 n^2) ... )

    This implements equation 5 in Nemes */
StirlingCoefficient[k] :=
{
   ak = (2k)!/(2^k k!)
   sum1 = 0
   for i = 0 to 2k
   {
      sum2 = 0
      sum1a = generalizedBinomial[k+i-1/2, i] generalizedBinomial[3k + 1/2, 2k-i] 2^i 
      for j = 0 to i
         sum2 = sum2 + binomial[i,j] (-1)^j j! S[2k + i + j, j] / (2k + i + j)!

      sum1 = sum1 + sum1a * sum2
   }

   return sum1 * ak
}

/** This is the equation after (5) in Nemes. */
S[p,q] :=
{
   s = 1/q!
   sum = 0
   for l = 0 to q
      sum = sum + (-1)^l binomial[q, l] (q-l)^p

   return s * sum
}

/** Testing */
/*for k = 0 to 8
   println["$k\t" + StirlingCoefficient[k]]
*/