Gamma(x) Lanzcos Davide

8103c8d0-0e01-4f6f-b62e-e3020f899022 24 15 12 decimal

&[DATE] &[TIME] - &[FILENAME]

z Lanczos Davde Carpi 20200514 ε 11015 -^* : p 57.156235665862923517 59.597960355475491248 - 14.136097974741747174 0.49191381609762019978 - .33994649984811888699 104 -^* .46523628927048575665 104 -^* .98374475304879564677 - (104 -^* .15808870322491248884 103 -^* .21026444172410488319 - (103 -^* .21743961811521264320 103 -^* .16431810653676389022 - (103 -^* .84418223983852743293 104 -^* .26190838401581408670 - (104 -^* .36899182659531622704 105 -^* 141 mat : g 4.7421875 : z Re 0.5 < y π π z * sin 1 z - Lanczos * / : z z 1 - : x 0.99999999999999709182 : n 1 p length range x x p n el z n + / + : for t z g + 0.5 + : y 2 π * sqrt t z 0.5 +^* t - exp * x * : 51 line if y Im abs ε ≤ y Re y if 61 line :

g=4.7421875, n=15

x Lanczos

x 2.5 i * :

x Lanczos 0.016129054047794 0.026749108945949 i * -

x Lanczos abs 0.031235576093277

2.5 Lanczos 1.32934038817912

2.5 Gamma 1.32934227535731

From Collab NDTMA ... otherwise un-documented α β X Excel 14 mat

x Gamma x Lanczos :

x α β Γ 1 α Gamma β α^* / x α 1 - (^* x β / - exp * :

L H N xd U 0 : dx H L - N / : 12 mat i 1 N range U i el L dx - (dx i * + eval : for U 31 line :

ε 1015 -^:

α 0.7352 :

β 245.4756 :

x α β Γ

ε α β Γ 131.67

L ε : H 10 : N 1000 : 13 mat

U L H N xd :

data Σ Finite diff's Cumul Integration vx data 1 col : vy data 2 col : Y 1 el 0 : i 2 vx rows range Δx i el vx i el vx i 1 - el - : Y i el Y i 1 - el 12 / vy i 1 - el vy i el + (Δx i el * (* + : 21 line for vx Y augment 61 line :

data j 1 U rows range Y j el U j el α β Γ eval : for U Y augment 21 line :

sub data 2 data rows 12 submatrix :

X 30.4 63.9 387.1 330 187.3 495.9 18 163.3 128.8 91 mat :

sub Σ sub

α 0.7352 : β 245.4756 : ε 1012 -^: 31 sys

eval(,) drops timing ~ 0.75
acc acuuracy 50/100 plentty.
Smath Gamma(x) ... 4 decimals x α β CDF 1 α Gamma β α^* / t α 1 -^t β / - exp * t ε x int eval * :

CDF X α β CDF vectorize :

So damned wrong ! CDF 1.292 2.612 14.4916 12.4422 7.2545 18.3727 0.7877 6.3677 5.0819 91 mat

Excel 0.2233 0.3647 0.8684 0.8291 0.6644 0.9194 0.1551 0.6224 0.5508 91 mat :

NOT even close
acc 10000, 5 min 0.2339 0.387 1.0043 0.9449 0.7301 1.0935 0.1614 0.6796 0.5959 91 mat

Int(f,a,b,n) ... Simpson/Romberg Kernel

Simpson Integrator in the Romberg Transform Domain [0,1] <= [a,b] ...  
... integrates  Smath ainterp(X,Y,x), cinterp(X,Y,x), cubic(x) from  l_p_c_spline f a b n Int w s b a - (34 / 2 w * 1 - (* 14 / 2 w * 1 - (3^* - 12 / + (* a + : w G w s f 12 w * 1 - (2^- (* : α 0 : β 1 : h β α - n / : 14 mat Simpson 13 / h 2 / α G β G + (* h α k h * + G k 1 n 1 - sum * + (* 23 / h * α h 2 / - k h * + G k 1 n sum * + eval : 32 / b a - (* Simpson * 41 line :

t0 1 time :

x f x α β Γ :

1/1 ... The  4th decimals NOT to be argued n 100 :

f ε X n Int vectorize 0.2232 0.3646 0.8682 0.8289 0.6643 0.9192 0.1551 0.6223 0.5507 91 mat

Excel 0.2233 0.3647 0.8684 0.8291 0.6644 0.9194 0.1551 0.6224 0.5508 91 mat

1 time t0 - 9 s *

Romberg Integrand in the Transform Domain [0,1] <= [a,b] ... f a b x Romberg x s b a - (34 / 2 x * 1 - (* 14 / 2 x * 1 - (3^* - 12 / + (* a + : x G x s f 12 x * 1 - (2^- (* : 21 line :

x c 10 x ≤ (x 1 ≤ (& cases :

Fill area

x g f α β x Romberg :

data spec pts i 1 data rows range p i el a b c spec : for v p 1 el : Unest 'p' j 2 p rows range v v p j el stack : for data v augment 51 line :

sx 0 : N 1000 : 12 mat

xmin sx : xmax 1 sx - : 12 mat ymin 0 : ymax 0.005 : 12 mat 21 line

space_x xmax xmin - (: space_y ymax ymin - (: 21 line

A function that tests an arbitrary point 
for its membership in the critical region. x y Test y x g ≤ (x xmax ≤ (* :

A routine that finds an array of N random 
pairs in the critical region. Pts i 1 N range x 1 Random space_x * xmin + 1 el : y 1 Random space_y * ymin + 1 el : 1 x y Test - x 1 Random space_x * xmin + 1 el : y 1 Random space_y * ymin + 1 el : 21 line while P i 1 el x : P i 2 el y : 51 line for P 21 line :

a b c spec . 6 #808080 13 mat :

fill Pts a b c spec pts :

Romberg integrand => Simpson(n)

f α β x Romberg x c * 1 0 + 40 red 1 5 mat fill 3 1 sys

Romberg-Simpson . . .

The function ' f ' to integrate is first transformed

in the Romberg domain [0 .. 1]. Then passed

to Simpson decimated ' n ' steps.

Simpson is well suited for this type of integrand .

]]>