(defun gamma (x)
  (if (< x 0.5)
    (/ %pi%
       (* (sin (* %pi% x))
          (gamma (- 1 x))))
    (let* ((cof #(676.5203681218851 -1259.1392167224028
                  771.32342877765313 -176.61502916214059
                  12.507343278686905 -0.13857109526572012
                  9.9843695780195716e-6 1.5056327351493116e-7))
           (ser0 0.99999999999980993)
           (z (pred x))
           (tmp (+ z (len cof) -0.5))
           (ser (+ ser0 (sum [mapcar / cof (succ z)]))))
      (* (sqrt (* 2 %pi%))
         (expt tmp (+ z 0.5))
         (exp (- tmp))
         ser))))

(each ((i (rlist 0.1..1.0..0.1 2..10)))
  (put-line (pic "##.# ######.######" i (gamma i))))