---------- Forwarded message ---------
From: Bill Gosper <billgosper(a)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
