---------- Forwarded message --------- From: Bill Gosper <billgosper@gmail.com> Date: Fri, Jun 19, 2020 at 10:41 AM Subject: Decently convergent gamma (factorial) formula To: Wolfram Tech Support Assuming you precompute the small Zetas, z! == E^((-EulerGamma)*z + Sum[((-z)^j*Zeta[j])/j, {j, 2, k}])*Product[ 1/(E^Inactive[Sum][(-(z/n))^j/j, {j, k}]*(1 + z/n)), {n, k}] Error ~ (z/k)^(1 + k)/(1 + k) , E.g., with k=9 terms, In[349]:= List @@ %304 /. z -> 1/4 /. k -> 9 // Activate Out[349]= {(1/4)!, (1/65211575) 33554432 E^( 6424827001954153772735623651396190479/ 9669127754221773701410455552000000000 - EulerGamma/4 + \[Pi]^2/ 192 + \[Pi]^4/92160 + \[Pi]^6/23224320 + \[Pi]^8/4954521600 - Zeta[3]/192 - Zeta[5]/5120 - Zeta[7]/114688 - Zeta[9]/2359296)} In[353]:= N@%349 Out[353]= {0.906402477055477, 0.906402477055477} The Inactive was to suppress the LerchPhi junk function. In[360]:= %304 /. k -> ∞ // Activate Out[360]= z! == Gamma[1 + z] Not great, but usable for complex z. I should investigate why it gets as many digits for z = ½ as for z = ¼. This is a good time to mention the (generalized) Remez algoritthm—an ingenious way to optimize an approximation over a fixed interval. E.g., for real Factorial you only need z ϵ [0,1] or [-½,½], e.g.. Introduce a new unknown, traditionally named 𝜇, and use all your numerical coefficients as first approximations to the solution of a system of simultaneous nonlinear equations stating that the error (absolute or relative, your choice) is alternately 𝜇 and -𝜇 at a bunch of turning points that you inititally imagine are distributed along the interval like the turning points of a Chebychev polynomial mapped onto that same interval. Notice that nowhere do we minimize anything. Just equalize the ripple amplitude. Depending on the form of the approximation, it can be hard as heck to get this process started. But once you get a solution, compute where the ripples really are, and solve the new system for new coefficients. Then the process will converge like a beast (quadratically). It may happen that the form of your approximation fortuitously has an extra turning point, and you have to manually choose which ripples to equalize. But when you succeed, the final |𝜇| will probably be dramatically small. And then explode as soon as you step ouside the interval. —rwg