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