(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")))