2

I have derived the following alternating infinite series as an approximation to values of the Gamma Function which yields reasonably accurate results between Gamma(1) and Gamma(7)-ish. Beyond that, the approximations oscillate and increasingly diverge from the true Gamma values. Do you see any opportunities to rein in this series’ unruly behavior?

EDIT: Another path of inquiry has arisen. Since posting this, I’ve entered my formula into a spreadsheet which calculates the sum of the first 170 terms of the series (the iPad spreadsheet app will only calculate up to 170!). It now occurs to me that perhaps the deviations I’ve observed for large m might indeed be accounted for by either accumulated rounding error or some other deficiency in the app’s underlying calculation methods.

In my spreadsheet, sqrt(exp(e + pi)) is rounded (or truncated, perhaps) to 18.7264552066483. The intriguing thing is that, for large m, when I make a change of just 1 unit in the 13th and last decimal place, this radically changes the resulting sum and hence, the percent error of my approximation. Now this might be wishful thinking, but if the series is that utterly sensitive to minuscule changes, perhaps the mounting errors for large m are an artifact of inadequate precision in processing the calculations.

$$\displaystyle \Gamma(m) \approx \sum_{n=0}^{\infty}(-1)^n \frac{\left(\sqrt{e^{e+\pi}}\right)^{m+n}}{n!(m+n)}$$

enter image description here

Christopher Emery
  • 71
  • 2
    A naive range reduction would use the formula $\Gamma(x+1) = x\Gamma(x)$, valid for $x\gt 0$ even if $x$ is not a whole number. – hardmath Dec 21 '21 at 23:56
  • Sorry I must edit the domain. It should read m between 1 and 7. – Christopher Emery Dec 21 '21 at 23:58
  • The argument of Gamma is ‘m’. The index of summation is ‘n’. – Christopher Emery Dec 22 '21 at 00:06
  • 1
    It's hard to tell whether you have an approximation that you like well enough on $[1,7]$ and are mainly interested in extending one way or another to arguments somewhat above $7$ (but not well above), or you are interested in rethinking the approximation scheme from scratch to potentially get state-of-the-art numerical values for an extended range of arguments, perhaps up to several orders of magnitude more than $7$. – hardmath Dec 22 '21 at 01:06
  • @hardmath I suppose my initial curiosity surrounds the nature and cause of the alternating divergence for larger m, and if it could not be corrected by making the sqrt(exp(pi + e)) piece sensitive to n in some alternating fashion. My intuition is suggesting that such a sensitivity is still missing, one which would counteract the oscillating error that seems to orbit the true values of Gamma(m), and extend the useful range of the approximation. – Christopher Emery Dec 22 '21 at 01:48
  • @hardmath With respect to your question regarding the range of desired extension of the useful domain, for sufficiently large m, Stirling’s formula would be the obvious choice, and I’m not trying to reinvent that wheel. My curiosity is more surrounding the issues stated above. – Christopher Emery Dec 22 '21 at 01:55
  • I have now reviewed the results and am more convinced that the fluctuations reported earlier in Gamma(m) for large m values are due to rounding errors. I need assistance from mathematicians with superior calculating equipment. My spreadsheet seems to truncate decimal values beyond the 15th decimal place, which appears to be inadequate for evaluating Gamma(m) values beyond approx 7. – Christopher Emery Dec 22 '21 at 07:02
  • Your series doesn't actually converge to $\Gamma(m)$, and it takes quite a while to converge for large $m$, so I don't quite understand why it's interesting. ;) Using SageMath, it seems to be stable (once $n$ is large enough), even for large $m$. I'll post an answer shortly. – PM 2Ring Dec 22 '21 at 08:15
  • 1
    with $m=7$, Mathematica gives $719.5465729540$ and that is accurate to 10 places past the decimal. So, as PM2Ring says, this does not converge to $\Gamma(7)=720$. – robjohn Dec 22 '21 at 08:39
  • @PM 2Ring Yes, I am aware that it does not converge to Gamma (m), but the approximation is better with the precision with which you have treated it. Clearly, the 15 place accuracy of my spreadsheet is woefully inadequate for calculating stable values of this series. I think with a minor, still unknown, tweak to the piece you call g(x), this series can be made to follow the Gamma Function closely. – Christopher Emery Dec 22 '21 at 14:45
  • @robjohn Ah, true, but for an approximation to work, it needn’t be equal to the function being approximated; it only needs to differ by a predictable amount, and then the necessary “tweak” is discoverable. I still suspect a small modification of the sqrt(exp(pi + e)) piece, call it ‘x’, making it some function x(n) will do the trick. Clearly, I’ll need to upgrade my computational resources, as obtaining good numbers is pretty essential in such an endeavor. – Christopher Emery Dec 22 '21 at 15:51
  • 1
    For the values of $m$ up to $7$, the terms of the series seem to peak at around $n=17$ or so. The terms are huge there; even for $m=1$, they reach $1.25\times10^7$, which removes $7$ significant figures. For $m=7$, they reach $4.05\times10^{14}$, which removes $12$ significant figures. – robjohn Dec 22 '21 at 17:46
  • @robjohn Yes, the 17th to 18th terms seem to represent the peak per-term values, and are certainly contributing to the loss of precision. What is needed here is a near-infinity computational machine, which leaves the nettlesome question, how near infinity is “near”? Unfortunately, no matter how close to infinity we are, we are still infinitely far away from it. – Christopher Emery Dec 23 '21 at 02:52

1 Answers1

2

To simplify the notation, I define $$g = \exp((e+\pi) / 2)$$ $$g\approx$$ 18.72645520664827426714604454007043332367490879335810853397969311801774145367
to 256 bits of precision.

Your series can be written as $$\Gamma(x)\approx \frac{g^x}x -\frac{g^{x+1}}{1!(x+1)} + \frac{g^{x+2}}{2!(x+2)} - \frac{g^{x+3}}{3!(x+3)} + \ldots $$

Here's some Sage code that calculates that series, to any desired precision. (SageMath is a free open-source mathematics software system based on Python and built on top of many existing open-source packages: NumPy, SciPy, matplotlib, Sympy, Maxima, GAP, FLINT, R).

@interact
def cegamma(m=4, N=100, bits=128, auto_update=False):
    RF = RealField(prec=bits)
    g = RF(exp((e + pi) / 2))
    p = g^m
    q = 1
    tot = p / m
    print("n: total term")
    for i in range(1, N + 1):
        p *= -g
        q *= i
        t = p / (q * (m + i))
        tot += t
        print(i, ":", tot, t)

Here's the output for the default inputs of $m=4$, $N=100$ terms, with 128 bits of precision used for the calculations. (Standard floating point arithmetic on most modern systems has 53 bits of precision).


n: total term
1 : -429838.83278896014300861543965159631773 -460582.97023748601574373014926172270191
2 : 3.1639471510015803333943129379871407637e6 3.5937859837905404764029283776387370815e6
3 : -1.6064302062636693178837366376603219359e7 -1.9228249213638273512231679314590360123e7
4 : 6.2702530225184220154606066851246874073e7 7.8766832287820913333443433227850093432e7
5 : -1.9952387983836339644371199456901536544e8 -2.6222641006354761659831806142026223952e8
6 : 5.3706178846991760697930046975336777013e8 7.3658566830828100342301246432238313558e8
7 : -1.2543198379500963597446743218553873574e9 -1.7913816264200139667239747916087551275e9
8 : 2.5895187624336366330550062651247962383e9 3.8438386003837329927996805869801835957e9
9 : -4.7931962500454989518084789767542955116e9 -7.3827150124791355848634852418790917499e9
10 : 8.0444970770996706206881342310685345778e9 1.2837693327145169572496613207822830089e10
11 : -1.2353459569352453944704943162072305455e10 -2.0397956646452124565393077393140840033e10
12 : 1.7488838981189667921984214727394399419e10 2.9842298550542121866689157889466704873e10
13 : -2.2970199358730821951359045364201677816e10 -4.0459038339920489873343260091596077235e10
14 : 2.8141404911118247679385964197751670798e10 5.1111604269849069630745009561953348614e10
15 : -3.2309489903571938212505599708984098665e10 -6.0450894814690185891891563906735769464e10
16 : 3.4904849174656004699138438579359120775e10 6.7214339078227942911644038288343219441e10
17 : -3.5609790040632747637000007438029677474e10 -7.0514639215288752336138446017388798249e10
18 : 3.4416154116537558872627148582265409829e10 7.0025944157170306509627156020295087302e10
19 : -3.1600847129171983681838518307084261619e10 -6.6017001245709542554465666889349671448e10
20 : 2.7636822837943705157740568687226448857e10 5.9237669967115688839579086994310710477e10
21 : -2.3074563364817698140983228709125845620e10 -5.0711386202761403298723797396352294477e10
22 : 1.8430878164598942554200644196586507067e10 4.1505441529416640695183872905712352687e10
23 : -1.4110981066403108713858518183988657488e10 -3.2541859231002051268059162380575164555e10
24 : 1.0373586000890293632524325654542699095e10 2.4484567067293402346382843838531356583e10
25 : -7.3343535250623563553139088326158198152e9 -1.7707939525952649987838234487158518910e10
26 : 4.9946223130259442905834945012879017058e9 1.2328975838088300645897403333903721521e10
27 : -3.2805753133471535049059884590121326661e9 -8.2751976263730977954894829603000343719e9
28 : 2.0809410348331769319354430387277408418e9 5.3615163481803304368414314977398735079e9
29 : -1.2762901704804745829551774755211036605e9 -3.3572312053136515148906205142488445024e9
30 : 7.5770817549407482764698619817323349593e8 2.0339983459745494106021636736943371564e9
31 : -4.3588231589194238635005828848637122812e8 -1.1935904913860172139970444866596047240e9
32 : 2.4320636511565691701866693310800409250e8 6.7908868100759930336872522159437532061e8
33 : -1.3173979015404333901490397256517594779e8 -3.7494615526970025603357090567318004029e8
34 : 6.9337809008365433116973182237718570800e7 2.0107759916240877213187715480289451859e8
35 : -3.5488480254241547174434218342268647076e7 -1.0482628926260698029140740057998721788e8
36 : 1.7676775026216243345922401053015185351e7 5.3165255280457790520356619395283832427e7
37 : -8.5749526352079741427463085865512976093e6 -2.6251727661424217488668709639566482960e7
38 : 4.0539157157321214471807171590012216690e6 1.2628868350940095589927025745552519278e7
39 : -1.8690093733260692827102063205194826398e6 -5.9229250890581907298909234795207043088e6
40 : 840855.30225404567305320762756799592626 2.7098646755801149557634139480874785660e6
41 : -369351.21069256523972873242999209032980 -1.2102065129466109127819400575600862561e6
42 : 158510.85581726505688611867839767971877 527862.06650983029661485110838977004857
43 : -66481.388524207453414161252488088903761 -224992.24434147251030027993088576862253
44 : 27280.655850167965564546677350068619713 93762.044374375418978707929838157523474
45 : -10941.509571557402562477419495410709441 -38222.165421725368127024096845479329154
46 : 4307.4111957821591485495185442480168283 15248.920767339561711026938039658726269
47 : -1649.1643507199879751140989982501147648 -5956.5755465021471236636175424981315930
48 : 630.01145175839145589634871433145870693 2279.1758024783794310104477125815734717
49 : -224.59230095850356740344122079371073032 -854.60375271689502329978993512516943725
50 : 89.554381048304585544778889360915961100 314.14668200680815294822011015462669142
51 : -23.698418781126538432703232651073849551 -113.25279982943112397748212201198981065
52 : 16.358343432758760080037753193029554883 40.056762213885298512740985844103404434
53 : 2.4534162079940082391251375941237511826 -13.904927224764751840912615598905803701
54 : 7.1923146176673447229520362012211900768 4.7388984096733364838268986070974388942
55 : 5.6061572708996837794289209759232660508 -1.5861573467676609435231152252979240260
56 : 6.1277296416709635162891264946305925489 0.52157237077127973686020551870732649809
57 : 5.9591843157300883547699982440163650360 -0.16854532594087516151912825061422751294
58 : 6.0127248180348229651624209668492871192 0.053540502304734610392422722832922083275
59 : 5.9960009340624350703187837468793873896 -0.016723883972387894843637219969899729648
60 : 6.0011390280826633541177782518391173891 0.0051380940202282837989945049597299994886
61 : 5.9995859460142874278252527808576453053 -0.0015530820683759262925254709814720837853
62 : 6.0000479308422614938911503729529959594 0.00046198482797406606589759209535065411236
63 : 5.9999126576083463808496812369208468701 -0.00013527323391511304146913603214908932698
64 : 5.9999516565987928831111204125743468988 0.000038998990446502261439175653500028720287
65 : 5.9999405838510462469587343411194082222 -0.000011072747746636152386071454938676577467
66 : 5.9999436806862656714022487437920059774 3.0968352194244435144026725977551671362e-6
67 : 5.9999428273139256830517166984206840620 -8.5337233998835053204537132191540880087e-7
68 : 5.9999430590593018791908582675101932790 2.3174537619613914156908950921697851866e-7
69 : 5.9999429970256712892942121402218014410 -6.2033630589896646127288391838003541873e-8
70 : 5.9999430133966964499052153368511912929 1.6371025160611003196629389851968101246e-8
71 : 5.9999430091363633170554608250162157294 -4.2603331328497545118349755635128309536e-9
72 : 5.9999430102298520448201324167578584566 1.0934887277646715917416427271717553341e-9
73 : 5.9999430099529858744520564859001326373 -2.7686617036807593085772581933860562069e-10
74 : 5.9999430100221514304461347641795206497 6.9165555994078278279388012470218470110e-11
75 : 5.9999430100051003581219237372584625857 -1.7051072324211026921058064044656936620e-11
76 : 5.9999430100092492372709809695333203813 4.1488791490572322748577955894660585825e-12
77 : 5.9999430100082526838056275333439306802 -9.9655346535343618938970109048843438495e-13
78 : 5.9999430100084890213635117844530337075 2.3633755788425110910302736249415810407e-13
79 : 5.9999430100084336739935724552001836186 -5.5347369939329252850088924801725660612e-14
80 : 5.9999430100084464755089965467416642863 1.2801515424091541480667704605460067997e-14
81 : 5.9999430100084435507350773044398522116 -2.9247739192423018120747449311391490584e-15
82 : 5.9999430100084442109031235328072616165 6.6016804622836740940488878701243594042e-16
83 : 5.9999430100084440636680824299380556093 -1.4723504110286920600719908286659206992e-16
84 : 5.9999430100084440961187811555042218712 3.2450698725566166261953579766801133379e-17
85 : 5.9999430100084440890498562463739814729 -7.0689249091302403983730043904293207039e-18
86 : 5.9999430100084440905720081332769404616 1.5221518869029589887087497984043508652e-18
87 : 5.9999430100084440902479705153026690163 -3.2403761797427144528089936518594449819e-19
88 : 5.9999430100084440903161764083834548604 6.8205893080785844073251891880389663866e-20
89 : 5.9999430100084440903019795467973942664 -1.4196861586060594028937402987234883982e-20
90 : 5.9999430100084440903049020871057477514 2.9225403083534850822265471086800892166e-21
91 : 5.9999430100084440903043070021900327111 -5.9508491571504036161865649790148701791e-22
92 : 5.9999430100084440903044268690318281739 1.1986684179546284735575989282355608047e-22
93 : 5.9999430100084440903044029815049253990 -2.3887526902774914150341658456951760038e-23
94 : 5.9999430100084440903044076917615631123 4.7102566377132806714685796123251943412e-24
95 : 5.9999430100084440903044067726517094260 -9.1910985368629824302750048538538060978e-25
96 : 5.9999430100084440903044069501470511968 1.7749534177078185581174813598104523538e-25
97 : 5.9999430100084440903044069162197411318 -3.3927310064910409360666964909154428771e-26
98 : 5.9999430100084440903044069226392252728 6.4194841410037820600026173212230100362e-27
99 : 5.9999430100084440903044069214367297755 -1.2024954973092304107109356238914809193e-27
100 : 5.9999430100084440903044069216597493179 2.2301954238994490398931879464084125560e-28

Using the SageMath RealField class for the arithmetic, the series appears to be quite stable, once $N$ is large enough, even at relatively low precision.

Here's a live link to the script, running on the SageMathCell server, so you can run it in your browser, even on an iPad or phone.

FWIW, Sage gives this expression for the sum of your series:

$$\frac{{\left({\left(x + 1\right)} e^{\left(\frac{1}{2} \, \pi x + \frac{1}{2} \, x e\right)} \gamma\left(x + 2, -e^{\left(\frac{1}{2} \, \pi + \frac{1}{2} \, e\right)}\right) - {\left(x^{2} + x\right)} e^{\left(\frac{1}{2} \, \pi x + \frac{1}{2} \, x e\right)} \gamma\left(x + 1, -e^{\left(\frac{1}{2} \, \pi + \frac{1}{2} \, e\right)}\right)\right)} e^{\left(-\frac{1}{2} \, \pi - \frac{1}{2} \, e\right)}}{x^{2} \left(-e^{\left(\frac{1}{2} \, \pi + \frac{1}{2} \, e\right)}\right)^{x}}$$

where $$\gamma(a,z)=\int_0^z t^{a-1}e^{-t}\,\mathrm{d}t$$ is the lower incomplete gamma function

It turns out that the expression for the sum of your series can be simplified to $$\gamma(x, g)$$

I found that by telling Sage to use a variable for $g$ rather than its symbolic expression in terms of $\pi$ and $e$.

var('x,g,n')
f = (-1)^n * g^(x + n) / (factorial(n) * (x + n))
ceg = sum(f, n, 0, oo)
print(ceg)
show(latex(ceg))

live link

And here's an updated version of the main script which prints the value calculated using the lower incomplete gamma function ("gil") before it starts summing the series. This value doesn't suffer from the loss of precision caused by the large terms in your series, so its value differs from the sum of the series in the last 10 digits when calculating the series for $m=4$.

Here's a plot of $\Gamma(x)-\gamma (x,g)$ over $[0,4]$

Plot of Gamma(x)-gamma(x,g)

created using this plotting script, which can plot the difference over any specified domain.

The Wikipedia article on the incomplete gamma functions gives generalised continued fractions for both the upper and lower incomplete gamma functions. This script calculates both functions using those continued fractions, and sums them to give an estimate of gamma.

PM 2Ring
  • 5,544
  • 1
    I appreciate the legwork. I will have to try out this SageMath gizmo, as it appears to have the computational power I need to probe these matters. A casual examination of the series will reveal that it is the Maclaurin series for 1/e^x, multiplied through by x^m, then integrated with respect to x, term by term. This leaves only the question of what is an appropriate value (or function, x(n)) of x. – Christopher Emery Dec 22 '21 at 15:32
  • 1
    @ChristopherEmery: Since your Question is about that series, adding a paragraph to the Question's body about how you derived it, perhaps not much more detailed than your Comment above, would substantially improve your post. – hardmath Dec 23 '21 at 16:15
  • 1
    @hardmath Yes, thank you. I will add a little motivation/derivation paragraph when I get a spare moment. – Christopher Emery Dec 23 '21 at 17:40
  • @Christopher I've just updated my answer. – PM 2Ring Dec 24 '21 at 01:18
  • @PM 2Ring The difference graph looks surprisingly smooth and suspiciously exponential. That’s mildly encouraging because it implies that the series expression is perhaps missing some small exponential term which might correct its divergence from Gamma. – Christopher Emery Dec 30 '21 at 15:13
  • @Christopher Well, the difference between them equals the upper incomplete gamma function, which appears to be better behaved than the lower, despite the infinity in the upper limit of its integral. I'll add a script to my answer that computes it via a continued fraction expansion. – PM 2Ring Dec 30 '21 at 18:55