63 lines
2.6 KiB
Scheme
63 lines
2.6 KiB
Scheme
(import (scheme base)
|
|
(scheme inexact)
|
|
(scheme write))
|
|
|
|
(define PI 3.14159265358979323846264338327950)
|
|
(define e 2.7182818284590452353602875)
|
|
|
|
(define gamma-lanczos
|
|
(let ((p '(676.5203681218851 -1259.1392167224028 771.32342877765313
|
|
-176.61502916214059 12.507343278686905 -0.13857109526572012
|
|
9.9843695780195716e-6 1.5056327351493116e-7)))
|
|
(lambda (x)
|
|
(if (< x 0.5)
|
|
(/ PI (* (sin (* PI x)) (gamma-lanczos (- 1 x))))
|
|
(let* ((x2 (- x 1))
|
|
(t (+ x2 7 0.5))
|
|
(a (do ((ps p (cdr ps))
|
|
(idx 0 (+ 1 idx))
|
|
(res 0.99999999999980993 (+ res
|
|
(/ (car ps)
|
|
(+ x2 idx 1)))))
|
|
((null? ps) res))))
|
|
(* (sqrt (* 2 PI)) (expt t (+ x2 0.5)) (exp (- t)) a))))))
|
|
|
|
(define (gamma-stirling x)
|
|
(* (sqrt (* 2 (/ PI x))) (expt (/ x e) x)))
|
|
|
|
(define gamma-taylor
|
|
(let ((a (reverse
|
|
'(1.00000000000000000000 0.57721566490153286061
|
|
-0.65587807152025388108 -0.04200263503409523553
|
|
0.16653861138229148950 -0.04219773455554433675
|
|
-0.00962197152787697356 0.00721894324666309954
|
|
-0.00116516759185906511 -0.00021524167411495097
|
|
0.00012805028238811619 -0.00002013485478078824
|
|
-0.00000125049348214267 0.00000113302723198170
|
|
-0.00000020563384169776 0.00000000611609510448
|
|
0.00000000500200764447 -0.00000000118127457049
|
|
0.00000000010434267117 0.00000000000778226344
|
|
-0.00000000000369680562 0.00000000000051003703
|
|
-0.00000000000002058326 -0.00000000000000534812
|
|
0.00000000000000122678 -0.00000000000000011813
|
|
0.00000000000000000119 0.00000000000000000141
|
|
-0.00000000000000000023 0.00000000000000000002))))
|
|
(lambda (x)
|
|
(let ((y (- x 1)))
|
|
(do ((as a (cdr as))
|
|
(res 0 (+ (car as) (* res y))))
|
|
((null? as) (/ 1 res)))))))
|
|
|
|
(do ((i 0.1 (+ i 0.1)))
|
|
((> i 2.01) )
|
|
(display (string-append "Gamma ("
|
|
(number->string i)
|
|
"): "
|
|
"\n --- Lanczos : "
|
|
(number->string (gamma-lanczos i))
|
|
"\n --- Stirling: "
|
|
(number->string (gamma-stirling i))
|
|
"\n --- Taylor : "
|
|
(number->string (gamma-taylor i))
|
|
"\n")))
|