Gamma function: Difference between revisions
→{{header|Stata}}: update
(→{{header|Stata}}: update) |
|||
Line 3,776:
=={{header|Stata}}==
This implementation uses the Taylor expansion of 1/gamma(1+x). The coefficients were computed with
The results are compared to Mata's '''[https://www.stata.com/help.cgi?mf_gamma gamma]''' function for each real between 1/100 and 100, by steps of 1/100.
<lang stata>mata
_gamma_coef = 1.0,
5.772156649015328606065121e-1,
-6.558780715202538810770195e-1,
-4.200263503409523552900393e-2,
1.665386113822914895017008e-1,
-4.219773455554433674820830e-2,
-9.621971527876973562114922e-3,
7.218943246663099542395010e-3,
-1.165167591859065112113971e-3,
-2.152416741149509728157300e-4,
1.280502823881161861531986e-4,
-2.013485478078823865568939e-5,
-1.250493482142670657345359e-6,
1.133027231981695882374130e-6,
-2.056338416977607103450154e-7,
6.116095104481415817862499e-9,
5.002007644469222930055665e-9,
-1.181274570487020144588127e-9,
1.04342671169110051049154e-10,
7.782263439905071254049937e-12,
-3.696805618642205708187816e-12,
5.100370287454475979015481e-13,
-2.05832605356650678322243e-14,
-5.348122539423017982370017e-15,
1.226778628238260790158894e-15,
-1.181259301697458769513765e-16,
1.186692254751600332579777e-18,
1.412380655318031781555804e-18,
-2.298745684435370206592479e-19,
1.714406321927337433383963e-20
function gamma_(x_) {
external _gamma_coef
x = x_
while (x<0.5) y = y/x++
z = _gamma_coef[24]▼
while (x>1.5) y = --x*y
x--
for (i=
return(y/z)
}
function map(f,a) {
n = rows(a)
Line 3,816 ⟶ 3,835:
}
x=(1::
u=map(&gamma(),x)
v=map(&gamma_(),x)
max(abs((
end</lang>
'''Output'''
<pre>
Here follows the Python program to compute coefficients.
<lang python>from mpmath import mp
mp.dps = 50
def gamma_coef(n):
a = [mp.mpf(1), mp.mpf(mp.euler)]
for k in range(3, n + 1):
s = sum((-1)**j * mp.zeta(j) * a[k - j - 1] for j in range(2, k))
a.append((s - a[1] * a[k - 2]) / (1 - k * a[0]))
return a
def horner(a, x):
y = 0
for s in reversed(a):
y = y * x + s
return y
gc = gamma_coef(30)
def gamma_approx(x):
y = mp.mpf(1)
while x < 0.5:
y /= x
x += 1
while x > 1.5:
x -= 1
y *= x
return y / horner(gc, x - 1)
for x in gc:
print(mp.nstr(x, 25))</lang>
=={{header|Tcl}}==
| |||