(With exact remainder and optional phase parameter.) EulMac[fn,var lo,hi,trms,phase] Here are three terms of Stirling's formula: In[305]:= Assuming[n > 1, EulMac[Log[z], z, 1, n, 3]] Out[305]= 1 - n + n*Log[n] - Inactive[Sum][Log[z], {z, 2, n}] == (-(1/24))*Inactive[Sum][ ConditionalExpression[ (-1 + 3*¢*(1 + ¢)* (-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢]))/(15*¢^3*(1 + ¢)^3), Re[¢] > 0 || Re[¢] < -1 || NotElement[¢, Reals]], {¢, 1, -1 + n}] + Inactive[Sum][((-1)^¢*BernoulliB[1 + ¢]* Derivative[¢][-Log[1 + #1] + Log[n + #1] & ][0])/ (1 + ¢)!, {¢, 0, 3}] The self-evidently superfluous ConditionalExpression surrounds the remainder summand. It is typically textually cumbersome but numerically small. The Sums are inactive to discourage Mathematica from struggling to "simplify" them. Activating gives Log[Pochhammer] instead of Log[Factoria] if you don't call n a positive integer. In[311]:= Assuming[n > 1 && Element[n,Integers], FullSimplify /@ Activate[Out[305]]]; In[314]:= Out[311] /. ConditionalExpression[ x_, __] :> FullSimplify[x]; In[315]:= Solve[%, Log[n!]] Out[315]= {{Log[n!] -> (1/(360*n^3))* (-1 + 30*n^2 + 331*n^3 - 360*n^4 + 180*n^3*Log[n] + 360*n^4*Log[n] + 15*n^3*Sum[(1/(15*¢^3* (1 + ¢)^3))*(-1 + 3*¢*(1 + ¢)*(-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢])), {¢, 1, -1 + n}])}} In[318]:= Expand /@ %315[[1, 1]] Out[318]= Log[n!] -> 331/360 - 1/(360*n^3) + 1/(12*n) - n + Log[n]/2 + n*Log[n] + (1/24)* Sum[(1/(15*¢^3*(1 + ¢)^3))* (-1 + 3*¢*(1 + ¢)* (-1 + 10*¢*(1 + ¢)* (1 + 12*¢*(1 + ¢))) + 180*¢^3*(1 + ¢)^3* (1 + 2*¢)*(Log[¢] - Log[1 + ¢])), {¢, 1, -1 + n}] Note 331/360 vs ln (rt2 pi): In[44]:= N[{Log[Sqrt[2*Pi]], 331/360}] Out[44]= {0.918938533, 0.919444444} Testing on 8! In[319]:= %318 /. n -> 8; In[320]:= Expand /@ % Out[320]= Log[40320] -> Log[2] + Log[3] + Log[4] + Log[5] + Log[6] + Log[7] + Log[8] (because the remainder sum is exact). Punting it, In[321]:= %318 /. Sum[__] -> 0 Out[321]= Log[n!] -> 331/360 - 1/(360 n^3) + 1/(12 n) - n + Log[n]/2 + n Log[n] Testing with n=8, In[322]:= % /. n -> 8 Out[322]= Log[40320] -> -(1303169/184320) + (17 Log[8])/2 In[323]:= N[%] Out[323]= 10.6046029 -> 10.6051088 Repeating the analysis starting with Assuming[n > 1, EulMac[Log[z], z, 1, n, 5, 1/2]] (i.e., one more Bernoulli term and phase = 1/2) Out[46]= Log[n!] -> 628991/612360 - 31/(1260 (1 + 2 n)^5) + 7/(360 (1 + 2 n)^3) - 1/( 12 + 24 n) + Log[2] - n (1 + Log[2]) - (3 Log[3])/ 2 + (1/2 + n) Log[1 + 2 n] In[47]:= % /. n -> 8. Out[47]= 10.6046029 -> 10.6046256 with the ^-5 term contributing only In[48]:= %%[[2, 2]] /. {{n -> 8}, {n -> 8.}} Out[48]= {-31/1789019820, -1.73279243*10^-8} so most of the .00002 error resides in the (discarded) remainder term: Out[49]= -433349316889/62104545180 - 8 Log[2] - 5 Log[3]/2 - Log[4] - Log[5] - Log[6] - Log[7] - Log[8] + 17 Log[17]/2 In[50]:= N[%] Out[50]= 0.0000226652532 The function is EulMac[xp_, v_, lo_, hi_, deg_, y_: 0] := Integrate[xp /. v -> y + v, {v, lo, hi}] - Inactive[Sum][xp, {v, lo + 1, hi}] == (-1)^deg*Inactive[Sum][ Integrate[(deg + 1)*v^deg* D[xp /. v -> v*y + ¢ + 1, {v, deg + 1}] - BernoulliB[deg + 1, y + v]* D[xp /. v :> y + v + ¢, {v, deg + 1}], {v, 0, 1}], {¢, lo, hi - 1}]/(deg + 1)! + Inactive[Sum][(-1)^¢* Derivative[¢][Evaluate[(xp /. v -> # + hi) - (xp /. v -> # + lo)] &][y]* BernoulliB[¢ + 1, y]/(¢ + 1)!, {¢, 0, deg}] The motivation for porting was NSum's seeming inability to get so much as four digits of Sum 2 to oo of 1/Log[n!] - 1/Log[(2 n)!] - 1/Log[(-1 + 2 n)!] apparently due to NSum's internal Euler-Maclaurin failing to try Method -> "DoubleExponential" on the corresponding NIntegrate. I now believe the sum ~ 1.4424630245831935 . This is still not accurate enough for PSLQ or the like, and ISC humorously proposes 1/Log[2 + 1/(300 Khinchin Gamma[1/6])]. --rwg