Re: [math-fun] Gamma function (formulas collection)
wds>many formulas about gamma function are in my paper:http://rangevoting.org/WarrenSmithPages/homepage/gammprox.pdf _______________________________________________ Interesting stuff. Random remarks: P4 middle is garbled, at least in my Preview. Did you mean asymptosy vs asymptopia? P6, In[16]:= F[0, c] == Assuming[c > 0,Integrate[ x^(2*n + 1)/(1 + x^2)^(n + 1)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]] Out[16]= F[0,c] == (-Cos[c])*CosIntegral[c] + (1/2)*Sin[c]*(Pi - 2*SinIntegral[c]) In[12]:= G[0, c] == Assuming[c > 0, Integrate[ x^(2*n)/(1 + x^2)^(n + 1)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]] Out[12]= G[0, c] == CosIntegral[c]*Sin[c] + (1/2)*Cos[c]*(Pi - 2*SinIntegral[c]) In[14]:= H[0, c] == Assuming[c > 0, Integrate[ x^(2*n)/(1 + x^2)^(n)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]] Out[14]= H[0, c] == 1/c Eliminating G and H from the F recurrence, 4*(n-1)*n*F(n)=2*(n-1)*(6*n-5)*F(n-1)-(12*n^2-32*n+c^2+22)*F(n-2)+2*(n-2)*(2*n-3)*F(n-3) Mma gets only MeijerG for n>0 in all of these. P7 (40), see http://www.tweedledum.com/rwg/idents.htm near bottom of 1st screen. (42) and (43) raise the interesting question of evaluating Product[((1/(12*n)) + 1)*((n - 1/2)*(E^(1/n) - 1))^n,{n, 1, Infinity}] (i.e., prod(((n-1/2)*(%e^(1/n)-1))^n*(1+1/12/n),n,1,inf) --"cwg"
Sho nuff, Mma doesn't know this class of MeijerG special values: {F[1, c] -> (1/4)*(-2 + Cos[c]*(-4*CosIntegral[c] + c*(Pi - 2*SinIntegral[c])) + 2*Sin[c]*(Pi + c*CosIntegral[c] - 2*SinIntegral[c])), F[2, c] -> (1/16)*(-12 + 2*CosIntegral[c]*((-8 + c^2)*Cos[c] + 7*c*Sin[c]) + 7*c*Cos[c]*(Pi - 2*SinIntegral[c]) - (-8 + c^2)*Sin[c]*(Pi - 2*SinIntegral[c])), F[3, c] -> (1/96)*(2*(-44 + c^2) + 6*(-16 + 5*c^2)*Cos[c]*CosIntegral[c] - 2*c*(-57 + c^2)*CosIntegral[c]*Sin[c] - c*(-57 + c^2)*Cos[c]*(Pi - 2*SinIntegral[c]) - 3*(-16 + 5*c^2)*Sin[c]*(Pi - 2*SinIntegral[c]))} (These are sufficient to initialize the 3rd order recurrence, below.) As in Warren's paper, F[n,c] is "just" a MeijerG/(2 n! sqrt(π)): {F[0, c] -> (1/2)*(-2*Cos[c]*CosIntegral[c] + Sin[c]*(Pi - 2*SinIntegral[c])), F[1, c] -> MeijerG[{{-1}, {}}, {{0, 0, 1/2}, {}}, c^2/4]/(2*Sqrt[Pi]), F[2, c] -> MeijerG[{{-2}, {}}, {{0, 0, 1/2}, {}}, c^2/4]/(4*Sqrt[Pi]), F[3, c] -> MeijerG[{{-3}, {}}, {{0, 0, 1/2}, {}}, c^2/4]/(12*Sqrt[Pi])} F[n,c] is a polynomial of first degree in Sin[c], Cos[c], SinIntegral[c], CosIntegral[c], and π, and (for n>1) degree 2*Floor[n/2] in c itself, with a messy profusion of crossterms. Maybe we could do better with a finite Abel-Plana using the function log(n). E.g., Out[128]= Log[m!] == -m - 2*Integrate[ArcTan[(m*y)/(1 + m + y^2)]/ (-1 + E^(2*Pi*y)), {y, 0, Infinity}] + (1/2 + m)*Log[1 + m] In[129]:= % /. m -> -1/2 Out[129]= Log[Pi]/2 == 1/2 - 2*Integrate[-(ArcTan[y/(2*(1/2 + y^2))]/(-1 + E^(2*Pi*y))), {y, 0, Infinity}] In[130]:= N[List @@ %] Out[130]= {0.572365, 0.572365} But I guess that just brings us pretty well back to Binet. The argument of the atan is <1 for the whole 0<y<oo, so integrating the usual atan powerseries produces nice convergence, but the termwise integrals are mysterious, even for m=1. Peachy source of definite integral identities, though: Log[(1/4)!] == (1/4)*(-1 - 8*Integrate[ArcTan[y/(5 + 4*y^2)]/(-1 + E^(2*Pi*y)), {y, 0, Infinity}] + 3*Log[5/4]) On Fri, Dec 30, 2011 at 1:57 PM, Bill Gosper <billgosper@gmail.com> wrote:
wds>many formulas about gamma function are in my paper:http://rangevoting.org/WarrenSmithPages/homepage/gammprox.pdf
_______________________________________________
Interesting stuff. Random remarks: P4 middle is garbled, at least in my Preview. Did you mean asymptosy vs asymptopia?
P6, In[16]:= F[0, c] == Assuming[c > 0,Integrate[
x^(2*n + 1)/(1 + x^2)^(n + 1)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]]
Out[16]= F[0,c] == (-Cos[c])*CosIntegral[c] + (1/2)*Sin[c]*(Pi - 2*SinIntegral[c])
In[12]:= G[0, c] == Assuming[c > 0, Integrate[
x^(2*n)/(1 + x^2)^(n + 1)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]]
Out[12]= G[0, c] == CosIntegral[c]*Sin[c] + (1/2)*Cos[c]*(Pi - 2*SinIntegral[c])
In[14]:= H[0, c] == Assuming[c > 0, Integrate[
x^(2*n)/(1 + x^2)^(n)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]]
Out[14]= H[0, c] == 1/c
Eliminating G and H from the F recurrence, 4*(n-1)*n*F(n)=2*(n-1)*(6*n-5)*F(n-1)-(12*n^2-32*n+c^2+22)*F(n-2)+2*(n-2)*(2*n-3)*F(n-3)
Mma gets only MeijerG for n>0 in all of these.
P7 (40), see http://www.tweedledum.com/rwg/idents.htm near bottom of 1st screen.
(42) and (43) raise the interesting question of evaluating
Product[((1/(12*n)) + 1)*((n - 1/2)*(E^(1/n) - 1))^n,{n, 1, Infinity}]
(i.e., prod(((n-1/2)*(%e^(1/n)-1))^n*(1+1/12/n),n,1,inf)
--"cwg"
Man, it's hard to get even a riesable value. .90058064316485891020 is probably only good to 12 digits. --rwg
DanA> Some other values might be checked. E.g., since Gamma(z) Gamma(1-z) == pi/sin(pi z), setting z = 1/2 + i gives |Gamma(1/2 + i)|^2 = 2pi /(e^pi + e^-pi). --Dan So again, the big mystery is the arg. Plot[Arg[Gamma[1/2 + t*I]], {t, -4, 4}] shows a nice smooth curve vaguely resembling x^3-7x. What's the slope at the mystery point t=1? In[260]:= D[Arg[Gamma[1/2 + t*I]], t] Out[260]= I*Gamma[1/2 + I*t]*PolyGamma[0, 1/2 + I*t]* Derivative[1][Arg][Gamma[1/2 + I*t]] In[261]:= N[% /. t -> 1] Out[261]= -0.772374 + 1.01843 I A perfectly smooth, perfectly real function has a non-real derivative. Is this Jan 1 or April 1? Macsyma does this easily: (c111) RATSIMP(NEGFACTORIAL(SUBST(1,T,DIFF(%,T)))); 2 %i + 1 2 %i - 1 psi (--------) + psi (- --------) 0 2 0 2 (d111) --------------------------------- 2 (c112) DFLOAT(%); (d112) - 0.05176165099441d0 --rwg In a different reply I said, (42) and (43) raise the interesting question of evaluating Product[((1/(12*n)) + 1)*((n - 1/2)*(E^(1/n) - 1))^n,{n, 1, Infinity}] (i.e., prod(((n-1/2)*(%e^(1/n)-1))^n*(1+1/12/n),n,1,inf) This is similar to the age-old question of Product[Cos[x/n], {n, 1, \[Infinity]}] and will probably require at least one new special function.
[...] (42) and (43) raise the interesting question of evaluating
Product[((1/(12*n)) + 1)*((n - 1/2)*(E^(1/n) - 1))^n,{n, 1, Infinity}]
(i.e., prod(((n-1/2)*(%e^(1/n)-1))^n*(1+1/12/n),n,1,inf)
--"cwg"
Man, it's hard to get even a riesable value. .90058064316485891020 is
probably only good to 12 digits. --rwg
Foo, it's easy. Expand the prodand at oo (= Foo/F): In[259]:= Series[(1/(12*n) + 1)*((n - 1/2)*(E^(1/n) - 1))^n, {n, Infinity, 3}] Out[259]= (E^ SeriesData[n, Infinity, {Log[2], 0, -1/12, -1/24}, -1, 3, 1]* SeriesData[n, Infinity, {1, 1/12}, 0, 4, 1])/2^n That's not a series! In[260]:= Series[Normal[%], {n, Infinity, 3}] Out[260]= SeriesData[n, Infinity, {1, 0, -13/288, 1/5184}, 0, 4, 1] That's a series, and the n^-3 term is wrong! In[261]:= Series[Normal[Series[(1/(12*n) + 1)*((n - 1/2)*(E^(1/n) - 1))^n, {n, Infinity, 4}]], {n, Infinity, 3}] Out[261]= SeriesData[n, Infinity, {1, 0, -13/288, -409/25920}, 0, 4, 1] In[262]:= Normal[%] Out[262]= 1 - 409/(25920*n^3) - 13/(288*n^2) Construct two factors from this: In[263]:= %[[1 ;; 2]]*%[[{1, 3}]] Out[263]= (1 - 409/(25920*n^3))*(1 - 13/(288*n^2)) Dividing these out of the prodand In[264]:= Out[261]/% Out[264]= SeriesData[n, Infinity, {1}, 0, 4, 1] raises the convergence rate to n^-4. And we have a closed form for what we divided out: In[265]:= Product[(1 + 1/12/n)*((n - 1/2)*(E^(1/n) - 1))^n, {n, Infinity}] == Product[(1 + 1/12/n)*(((n - 1/2)*(E^(1/n) - 1))^n/%%), {n, Infinity}]* Product[%%, {n, Infinity}] Out[265]= Product[(1 + 1/(12*n))*((-1 + E^(1/n))*(-(1/2) + n))^n, {n, Infinity}] == (12*Sqrt[2/13]* Product[((1 + 1/(12*n))*((-1 + E^(1/n))*(-(1/2) + n))^n)/ ((1 - 409/(25920*n^3))*(1 - 13/(288*n^2))), {n, Infinity}]*Sin[(1/12)*Sqrt[13/2]*Pi])/ (Pi*Gamma[1 - (1/12)*(409/15)^(1/3)]* Gamma[1 + (1/24)*(409/15)^(1/3) - (I*(409/5)^(1/3))/(8*3^(5/6))]* Gamma[1 + (1/24)*(409/15)^(1/3) + (I*(409/5)^(1/3))/(8*3^(5/6))]) But o(n^-4) isn't quite enough for Mma to get 22 digits, say. Dividing out binomials for the next few terms got Julian 0.900580643164907990153824318444289111926734003669821962 + 0.*10^-55 I fairly quickly. Unfortunately, this only agrees with the results below to about 30D. Alternatively, you can (theoretically) eliminate arbitrarily many terms just by fooling with the first order term! E.g., with parameters Out[219]= {e -> (45240 - 2153*Sqrt[3770])/221676, f -> (45240 + 2153*Sqrt[3770])/221676, b -> (1/340)*(-110 + Sqrt[3770]), a -> (1/340)*(-110 - Sqrt[3770])} In[220]:= FullSimplify[ Series[Normal[ Series[(1 + a/n)^e*(1 + b/n)^f*((n - 1/2)*(E^(1/n) - 1))^n /. %, {n, \[Infinity], 6}]], {n, \[Infinity], 5}]] Out[220]= 16633 SeriesData[n, Infinity, {1, 0, 0, 0, 0, -(---------)}, 0, 6, 1] 308448000 is O(n^-5). We have Product[(1 + 1/12/n)*((n - 1/2)*(E^(1/n) - 1))^n, {n, \[Infinity]}]-> Product[(1 + a/n)^e*(1 + b/n)^f*((n - 1/2)*(E^(1/n) - 1))^n, {n, \[Infinity]}]/ Product[(1 + a/n)^e*(1 + b/n)^f/(1 + 1/12/n), {n, \[Infinity]}] Plugging in a,b,f,g, the rational function product becomes a pile of Gammas, and In[268]:= N[%[[2]], 33] Out[268]= 0.900580643164907990153824318444627 took just a few seconds. Here are two further runs with timings: Timing[N[%%[[2]],44]] {19.5682, 0.9005806431649079901538243184443488363691376176930701107} Timing[N[%%%[[2]],69]] {20.0368, 0.900580643164907990153824318444348836369137617068042914805385390880213} NProduct seems to be suffering from the cancellation in e^(1/n)-1. For integer k>1, Julian worked out the formula for Product[1-a/n^k,{n, 1, Infinity}] This would be especially interesting if we could find among the profusion of ak, ek, and fk satisfying Product[1-ak/n^ek)^fk,{k,oo}] = ((n - 1/2)*(E^(1/n) - 1))^n any which are predicable functions of k. --rwg
I missed where Warren's paper probably says that, for "unrestricted" (arbitrary precision) numerics, the divergence of the log Gamma and polygamma series is rather irrelevant. A plausible algorithm is simply to add some positive integer n to z so that the Stirling series z+n swoops within the desired accuracy. Then divide by (z+1)(z+2)...(z+n). Windschitl's and Warren's improvements make this even more attractive. Looking to provide Gamma numerics for Robert's RIES, I accidentally overkilled the problem by Remezing up f[z]:=(0.011987912746518014561506488044010 + z^2 (-0.028678142684262583619372798888834 + z^2 (-0.10376325088888750555218325558841 + z^2 (-0.034612435282031905677869242849747 - 0.0024691358024691358024691358024691 z^2))))/(z^7 * (45.984169004371332413538940188341 + z^2 (54.768197084214553421696753890699 + z^2 (14.960893432079977923321408051895 + z^2)))) Then plotting z! /(Sqrt[2 Pi] (E^(-f[z] - 4 z) z^(2 + 6 z) Sinh[1/z]^(2 z))^(1/4) -1 9 < z < oo shows maximum error 6*10^-26 (at only one ripple because I neglected the weight function). So up to nine divides could be required for this method. But RIES only needs 10^-17, for which f[z] could be a lot simpler and n smaller than 9. To squeeze some actual fun out of this: Why is f an odd function? --rwg PS, I wasn't confused in that 2001 Remez posting--just confusing. In the next sentences I explained why truncating the Chebys doesn't minimax. I was just using "Chebychef optimal" to mean the best you could do by subtracting off Chebys. On Fri, Dec 30, 2011 at 1:57 PM, Bill Gosper <billgosper@gmail.com> wrote:
wds>many formulas about gamma function are
in my paper:http://rangevoting.org/WarrenSmithPages/homepage/gammprox.pdf
_______________________________________________
Interesting stuff. Random remarks: P4 middle is garbled, at least in my Preview. Did you mean asymptosy vs asymptopia?
P6, In[16]:= F[0, c] == Assuming[c > 0,Integrate[
x^(2*n + 1)/(1 + x^2)^(n + 1)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]]
Out[16]= F[0,c] == (-Cos[c])*CosIntegral[c] + (1/2)*Sin[c]*(Pi - 2*SinIntegral[c])
In[12]:= G[0, c] == Assuming[c > 0, Integrate[
x^(2*n)/(1 + x^2)^(n + 1)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]]
Out[12]= G[0, c] == CosIntegral[c]*Sin[c] + (1/2)*Cos[c]*(Pi - 2*SinIntegral[c])
In[14]:= H[0, c] == Assuming[c > 0, Integrate[
x^(2*n)/(1 + x^2)^(n)*E^(-c*x) /. n -> 0, {x, 0, \[Infinity]}]]
Out[14]= H[0, c] == 1/c
Eliminating G and H from the F recurrence, 4*(n-1)*n*F(n)=2*(n-1)*(6*n-5)*F(n-1)-(12*n^2-32*n+c^2+22)*F(n-2)+2*(n-2)*(2*n-3)*F(n-3)
Mma gets only MeijerG for n>0 in all of these.
P7 (40), see http://www.tweedledum.com/rwg/idents.htm near bottom of 1st screen.
(42) and (43) raise the interesting question of evaluating
Product[((1/(12*n)) + 1)*((n - 1/2)*(E^(1/n) - 1))^n,{n, 1, Infinity}]
(i.e., prod(((n-1/2)*(%e^(1/n)-1))^n*(1+1/12/n),n,1,inf)
--"cwg"
I don't know how many folks consider practical numerical computation to be "fun", but I do (-: I've been looking for a numeric Gamma method that gives 34 digits of precision suitable for quad-precision floating-point [2] although I would settle for 31 digits suitable for "double-double". On Mon, Jan 9, 2012 at 02:24, Bill Gosper <billgosper@gmail.com> wrote:
I missed where Warren's paper probably says that, for "unrestricted" (arbitrary precision) numerics, the divergence of the log Gamma and polygamma series is rather irrelevant. A plausible algorithm is simply to add some positive integer n to z so that the Stirling series z+n swoops within the desired accuracy. Then divide by (z+1)(z+2)...(z+n). Windschitl's and Warren's improvements make this even more attractive. [...]
In the 2006 paper from Warren Smith [1], (which is very good by the way, now that I've read it more carefully) I am only seeing Windschitl in one place, credited with the approximation: z! = z Gamma[z] = Sqrt[2*Pi*z] * [z/E]^z * [z*Sinh[1/z]]^[z/2] * [1+O(|z|^-5)] which (if the O(|z|^-5) part is taken at face value) means that for 34 digits, I'd need to ensure z>6300000 (6.3 million) and that clearly won't work. With the Stirling series for log(x!), the biggest term is x log(x), and the smallest is of magnitude BernoulliB[2*n]/[2*n*[2*n-1]*x^[2*n-1]] Back in 1998 I studied this for the case of 15 digits of precision, and I worked out that I could I sum 8 terms, with the first omitted term being the "+43867/244188*x^17" term, which gives 15 digits provided my argument x is at least 7. No problem, quite reasonable. But if I want near 34 digits of precision, it seems I need to sum the terms at least as far as the x^-23 term, with the residual being of order 657931/(300 x^25). With that number of terms, x needs to be at least 24. The worst case is when computing something like Gamma[0.6]. In that case I need to: - Multiply out the product P=(x+1)(x+2)... which has 23 terms (22 multiply plus 23 add) - Sum the terms of the series (12 terms: 12 add, 24 multiply and 12 divide) - compute Log[n] and Log[P] as always, plus the other lead terms (a few more primitive operations) That's an awful lot of operations, which is why I still want to use Lanczos' method instead, or Gosper's latest work if that turns out to give me the digits I need. The Lanczos method has the same requirement for two calls to Log[], but it appears I can get 34 digits of precision with about 17 terms, using 2 adds and one divide per term. For example, at [3] and [4] is a 15-digit precision version using 6 terms, which is 12 adds and 6 divides and Viktor Toth (see [5]) uses 80 terms to get about 165 digits of precision. I haven't multiplied out the matrices to find coefficients for 34-digit precision (not yet). Maybe someone can set me straight on this... - Robert More references and more flexible source code at: http://mrob.com/pub/ries/lanczos-gamma.html [1] http://rangevoting.org/WarrenSmithPages/homepage/gammprox.pdf [2] http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format [3] William Press et al., Numerical Recipes in C, section 6.1, function "gammln" (pp. 213-216 in the hardcover 2nd edition) [4] http://paulbourke.net/miscellaneous/functions/ (which gives the same code as [3]) [5] http://www.vttoth.com/CMS/projects/41-the-lanczos-approximation -- Robert Munafo -- mrob.com Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 - mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com
Hello mr Munafo, maybe this source code is of interest to you ? It is about the computation of the gamma function to quadruple precision, if I recall this problem has been solved, I mean there are well established routines for doing that in Fortran or C, you can look at the code here : http://people.sc.fsu.edu/~jburkardt/f77_src/specfun/specfun.html best regards, simon plouffe Le 09/01/2012 16:35, Robert Munafo a écrit :
I don't know how many folks consider practical numerical computation to be "fun", but I do (-:
I've been looking for a numeric Gamma method that gives 34 digits of precision suitable for quad-precision floating-point [2] although I would settle for 31 digits suitable for "double-double".
On Mon, Jan 9, 2012 at 02:24, Bill Gosper<billgosper@gmail.com> wrote:
I missed where Warren's paper probably says that, for "unrestricted" (arbitrary precision) numerics, the divergence of the log Gamma and polygamma series is rather irrelevant. A plausible algorithm is simply to add some positive integer n to z so that the Stirling series z+n swoops within the desired accuracy. Then divide by (z+1)(z+2)...(z+n). Windschitl's and Warren's improvements make this even more attractive. [...]
In the 2006 paper from Warren Smith [1], (which is very good by the way, now that I've read it more carefully) I am only seeing Windschitl in one place, credited with the approximation:
z! = z Gamma[z] = Sqrt[2*Pi*z] * [z/E]^z * [z*Sinh[1/z]]^[z/2] * [1+O(|z|^-5)]
which (if the O(|z|^-5) part is taken at face value) means that for 34 digits, I'd need to ensure z>6300000 (6.3 million) and that clearly won't work.
With the Stirling series for log(x!), the biggest term is x log(x), and the smallest is of magnitude BernoulliB[2*n]/[2*n*[2*n-1]*x^[2*n-1]]
Back in 1998 I studied this for the case of 15 digits of precision, and I worked out that I could I sum 8 terms, with the first omitted term being the "+43867/244188*x^17" term, which gives 15 digits provided my argument x is at least 7. No problem, quite reasonable.
But if I want near 34 digits of precision, it seems I need to sum the terms at least as far as the x^-23 term, with the residual being of order 657931/(300 x^25). With that number of terms, x needs to be at least 24.
The worst case is when computing something like Gamma[0.6]. In that case I need to:
- Multiply out the product P=(x+1)(x+2)... which has 23 terms (22 multiply plus 23 add) - Sum the terms of the series (12 terms: 12 add, 24 multiply and 12 divide) - compute Log[n] and Log[P] as always, plus the other lead terms (a few more primitive operations)
That's an awful lot of operations, which is why I still want to use Lanczos' method instead, or Gosper's latest work if that turns out to give me the digits I need.
The Lanczos method has the same requirement for two calls to Log[], but it appears I can get 34 digits of precision with about 17 terms, using 2 adds and one divide per term. For example, at [3] and [4] is a 15-digit precision version using 6 terms, which is 12 adds and 6 divides and Viktor Toth (see [5]) uses 80 terms to get about 165 digits of precision. I haven't multiplied out the matrices to find coefficients for 34-digit precision (not yet).
Maybe someone can set me straight on this...
- Robert
More references and more flexible source code at: http://mrob.com/pub/ries/lanczos-gamma.html
[1] http://rangevoting.org/WarrenSmithPages/homepage/gammprox.pdf
[2] http://en.wikipedia.org/wiki/Quadruple-precision_floating-point_format
[3] William Press et al., Numerical Recipes in C, section 6.1, function "gammln" (pp. 213-216 in the hardcover 2nd edition)
[4] http://paulbourke.net/miscellaneous/functions/ (which gives the same code as [3])
[5] http://www.vttoth.com/CMS/projects/41-the-lanczos-approximation
Thank you Mr Plouffe! I'll look it over more closely -- but right away I notice that all of the coefficients in both specfun.f90 and specfun.f have between 18 and 25 digits. There was a time when there were a variety of double-precision hardware implementations (most of them more precise than out present "double precision" but not nearly so precise as today's "quad precision"). For example, the Cray C90's double precision provided about 28 decimal digits. [1] For IEEE binary128, we would need to enhance the precision of the coefficients to 34 digits and perhaps generate new minimax models (for example). - Robert [1] I have a list of many machines, but unfortunately not many of the old mainframes, at http://mrob.com/pub/math/floatformats.html On Mon, Jan 9, 2012 at 12:05, Simon Plouffe <simon.plouffe@gmail.com> wrote:
Hello mr Munafo,
maybe this source code is of interest to you ?
It is about the computation of the gamma function to quadruple precision, if I recall this problem has been solved, I mean there are well established routines for doing that in Fortran or C, you can look at the code here : http://people.sc.fsu.edu/~**jburkardt/f77_src/specfun/**specfun.html<http://people.sc.fsu.edu/~jburkardt/f77_src/specfun/specfun.html>
best regards, simon plouffe
--
Robert Munafo -- mrob.com Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 - mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com
Do you want it to work for complex or only for real argument?
I've been looking for a numeric Gamma method that gives 34 digits of precision suitable for quad-precision floating-point [2] although I would settle for 31 digits suitable for "double-double".
By the way, on the theoretical front, one function that I did not study, but probably deserves study, is g(z) := Gamma(z+1/2)^(1/z). does this have a convergent continued fraction expansion? It has a divergent asymptotic series .36787944117144232159552377016146086744581113103177*z+ .33805859406623990237027945096151887426851375834020+ .13999922423542209328692750196611434611925705150275/z+ .33493044087723619416023793118971174508783220448207e-1/z^2+ .56720131162423401908596776888208182906669281170266e-2/z^3+ .11415650363382196952422045360200835828013297526599e-2/z^4+ .40066835176315930489725300869291932605805052785655e-4/z^5- .18459730484462595896459837456234342865776092168523e-3/z^6+ .12180790117916786446998042861431093723893158454545e-3/z^7+ .17568169709328468660924962236969118264708016900776e-3/z^8- .23036114978935780189878504071470392725192413835113e-3/z^9- .26501793497134862116045528462277497265836478867190e-3/z^10+ .59068400115606570603275441175404496195447459946218e-3/z^11+ .62050016149812141546059352772177092271212715711371e-3/z^12- .20943016044158043700802700877762421178065590234503e-2/z^13- .21036447470556705667818040016973832680344857168705e-2/z^14+ .99842712425733569346391610042789127327983738573824e-2/z^15+ .97772628167677602737106584758217900245031666550502e-2/z^16- .61986382824479867397792991698201564463192370984501e-1/z^17- .59747962332351720667423814714163223752965683269949e-1/z^18+... but perhaps it becomes convergent when converted to a continued fraction (aka "Pade summation"). This would probably be hard to prove, if true.
I only need it to work for real arguments. But if that were an issue I could stick with Lanczos. Any given set of Lanczos parameters work equally well on the entire half-plane Real[z]>0, and with the same guaranteed error bound. This article describes it really well: http://www.haoli.org/nr/bookcpdf/c6-1.pdf which is from: William H. Press et al., Numerical recipes in C: the art of scientific computing (2nd edition). Cambridge University Press, (1992) ISBN 0521431085. On Mon, Jan 9, 2012 at 12:06, Warren Smith <warren.wds@gmail.com> wrote:
Do you want it to work for complex or only for real argument?
I've been looking for a numeric Gamma method that gives 34 digits of precision suitable for quad-precision floating-point [2] although I would settle for 31 digits suitable for "double-double".
-- Robert Munafo -- mrob.com Follow me at: gplus.to/mrob - fb.com/mrob27 - twitter.com/mrob_27 - mrob27.wordpress.com - youtube.com/user/mrob143 - rilybot.blogspot.com
WS = Warren Smith <warren.wds@gmail.com> writes:
WS> By the way, on the theoretical front, one function that I did not WS> study, but probably deserves study, is WS> g(z) := Gamma(z+1/2)^(1/z). If I got it right – and chose a reasonable version of clgamma() – there is src to generate a complex plot of that function, as well as a couple of graphs covering real(z) ∈ [-16,16] at: http://jhcloos.com/wsmith-gamma/ A most interesting look. In the src you’ll see that I implemented it as cexp(clngam(z+.5)/z). I’ll add a README file after getting some food. -JimC -- James Cloos <cloos@jhcloos.com> OpenPGP: 1024D/ED7DAEA6
participants (5)
-
Bill Gosper -
James Cloos -
Robert Munafo -
Simon Plouffe -
Warren Smith