Family Day update

Task/Gamma-function/Stata/gamma-function-1.stata

@@ -0,0 +1,61 @@
+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_
+    y = 1
+    while (x<0.5) y = y/x++
+    while (x>1.5) y = --x*y
+    z = _gamma_coef[30]
+    x--
+    for (i=29; i>=1; i--) z = z*x+_gamma_coef[i]
+    return(y/z)
+}
+
+function map(f,a) {
+    n = rows(a)
+    p = cols(a)
+    b = J(n,p,.)
+    for (i=1; i<=n; i++) {
+        for (j=1; j<=p; j++) {
+            b[i,j] = (*f)(a[i,j])
+        }
+    }
+    return(b)
+}
+
+x=(1::1000)/100
+u=map(&gamma(),x)
+v=map(&gamma_(),x)
+max(abs((v-u):/u))
+end

Task/Gamma-function/Stata/gamma-function-2.stata

@@ -0,0 +1,31 @@
+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))

Task/Gamma-function/Stata/gamma-function.stata

@@ -1,42 +0,0 @@
-mata
-_gamma_coef = 1.0,
-5.7721566490153286061e-1,-6.5587807152025388108e-1,
--4.2002635034095235529e-2,1.665386113822914895e-1,
--4.2197734555544336748e-2,-9.6219715278769735621e-3,
-7.2189432466630995424e-3,-1.1651675918590651121e-3,
--2.1524167411495097282e-4,1.2805028238811618615e-4,
--2.0134854780788238656e-5,-1.2504934821426706573e-6,
-1.1330272319816958824e-6,-2.0563384169776071035e-7,
-6.1160951044814158179e-9,5.0020076444692229301e-9,
--1.1812745704870201446e-9,1.0434267116911005105e-10,
-7.782263439905071254e-12,-3.6968056186422057082e-12,
-5.100370287454475979e-13,-2.0583260535665067832e-14,
--5.3481225394230179824e-15
-
-function gamma_(x_) {
-	external _gamma_coef
-	x = x_
-	for (y=1; x>2;) y = --x*y
-	z = _gamma_coef[24]
-	x--
-	for (i=23; i>=1; i--) z = z*x+_gamma_coef[i]
-	return(y/z)
-}
-
-function map(f,a) {
-	n = rows(a)
-	p = cols(a)
-	b = J(n,p,.)
-	for (i=1; i<=n; i++) {
-		for (j=1; j<=p; j++) {
-			b[i,j] = (*f)(a[i,j])
-		}
-	}
-	return(b)
-}
-
-x=(1::10000)/100
-u=map(&gamma(),x)
-v=map(&gamma_(),x)
-max(abs((u-v):/u))
-end