감마함수(gamma function)의 값을 (유효수자 15자리 까지 정밀하게) 계산하는 Common Lisp 언어 소스

Lanczos 알고리즘은 Stirlng 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다.  단지 exp 함수를 이용하는 부분에서는 exp 함수의 구현에 따라 오차가 더 있을 수 있다. 단, Common Lisp 언어의 float 타입은 유효숫자를 약 7개~8개 밖에 알아내지 못한다.

Comon Lisp 언어의 디퐁트 부동소수점수는 단정밀도 부동소수점수이다. 단정밀도 부동소수점수는 그 처리 가능한 양수의 범위가 최소 2⁻¹²⁶ ≈ 1.18 × 10⁻³⁸ 부터 최대  (2−2⁻²³) × 2¹²⁷ ≈ 3.4 × 10³⁸ 밖에 되지 않는다. 이 범위를 벗어나는 양수에 대해서는 overflow 에러가 나게 마련이다.

이런 이유로 100! 을 계산하기 위해 (gamma 101) 하면 floating point overflow 에러가 발생한다. 아래의 소스에서 (gamma 101) 을 구하는 부분을  ignore-errors 로 처리한 것은 floating point overflow 에러가 나서 실행이 멈추는 것을 방지하기 위함이다.

 

#!/usr/bin/env clisp

;; Filename: testLanczos-01.lsp
;;
;;           An approximation for the gamma function by using the Lanczos algorithm
;;
;; Execute: clisp testLanczos-01.lsp
;;         or
;; Execute: ./testLanczos-01.lsp
;;
;; See: http://en.wikipedia.org/wiki/Lanczos_approximation
;; See:http://www-history.mcs.st-and.ac.uk/Biographies/Lanczos.html
;; See: http://www.thinbasic.com/community/showthread.php?11279-Gamma-function

 
;; Coefficients used by the GNU Scientific Library
(setf *g* 7)
(setf *p* (list
            0.99999999999980993   676.5203681218851       -1259.1392167224028
          771.32342877765313     -176.61502916214059         12.507343278686905
           -0.13857109526572012     9.9843695780195716e-6    1.5056327351493116e-7) )
 
(defun gamma(z)
    (let ((tt (complex 0.0 0.0)))
        (if (realp z) (setf z (complex z 0.0)))
        ;; Reflection formula
        (if (< (realpart z) 0.5)
            (/ pi (sin (* (* pi z) (gamma (- 1 z)))))
            (progn
                (setf z (- z 1))
                (setf x (nth 0 *p*))
                (loop for i from 1 below (+ *g* 2) do
                    (setf x (+ x (/ (nth i *p*) (+ z i)))))
                (setf tt (+ z *g* 0.5))
                (* (sqrt (* 2 pi)) (expt tt (+ z 0.5)) (exp (- tt)) x) ))))

(defun factorial(n)
    (let ((k 0))
        (if (< n 2)
            1
            (progn
                (setf k 1)
                (if (evenp n)
                    (loop for i from 0 below (floor (/ n 2)) do
                        (setf k (* k (* (+ i 1) (- n i)))))
                    (progn
                        (loop for i from 0 below (floor (/ n 2)) do
                            (setf k (* k (* (+ i 1) (- n 1 i)))))
                        (setf k (* k n))))
                k))))

(defun facto(n)
    (let ((k 0))
        (if (< n 2)
            1
            (progn
                (setf k 1)
                (loop for i from 2 to n do
                    (setf k (* k i)) )
                k ))))

;; Begin here
(format t "gamma(10) = 9! = ~A asymptotically~%" (gamma 10))

(format t "gamma(21) = 20! = ~16A asymptotically~%" (gamma 21))
(format t "gamma(21) = 20! = ~16G asymptotically~%" (gamma 21))
(format t "gamma(21) = 20! = ~16G asymptotically~%" (realpart (gamma 21)))
(format t "gamma(21) = 20! = ~16F asymptotically~%" (gamma 21))
(format t "gamma(21) = 20! = ~16F asymptotically~%" (realpart (gamma 21)))
; (format t "gamma(101) = 100! = ~16A asymptotically~%" (gamma 101))
;; "{0.real:.15f}{0.imag:+.5f}j".format(gamma(101))

(ignore-errors
    (format t "gamma(101) = 100! = ~16A asymptotically~%" (gamma 101))
)

(setf i 10)
(format t "factorial(~D) = ~D! = ~D~%"  i i (factorial i))
(format t "facto(~D) = ~D! = ~D~%"  i i (facto i))

(setf i 20)
(format t "factorial(~D) = ~D! = ~D~%"  i i (factorial i))
(format t "facto(~D) = ~D! = ~D~%"  i i (facto i))

(setf i 100)
(format t "factorial(~D) = ~D! = ~D~%"  i i (factorial i))
(format t "facto(~D) = ~D! = ~D~%"  i i (facto i))

#|
if __name__ == "__main__":
    print "gamma(10) = 9! = %s asymptotically" % gamma(10)
    # print "gamma(101) = 100! = %16s asymptotically" % gamma(101)
    print "gamma(101) = 100! = %16s asymptotically" % "{0.real:.15f}{0.imag:+.5f}j".format(gamma(101))

    for i in range(11):
        print "%d! = %d" % (i, factorial(i))

    i = 100
    print "factorial(%d) = %d! = %d" % (i, i, factorial(i))
    print "facto(%d) = %d! = %d" % (i, i, facto(i))
|#

#|
Output:
gamma(10) = 9! = #C(362880.16 0.0) asymptotically
gamma(21) = 20! = #C(2.4329015E18 0.0) asymptotically
gamma(21) = 20! = #C(2.4329015E18 0.0) asymptotically
gamma(21) = 20! = 2432901500000000000.     asymptotically
gamma(21) = 20! = #C(2.4329015E18 0.0) asymptotically
gamma(21) = 20! = 2432901500000000000. asymptotically
factorial(10) = 10! = 3628800
facto(10) = 10! = 3628800
factorial(20) = 20! = 2432902008176640000
facto(20) = 20! = 2432902008176640000
factorial(100) = 100! = 93326215443944152681699238856266700490715968264381621468
59296389521759999322991560894146397615651828625369792082722375825118521091686400
0000000000000000000000
facto(100) = 100! = 933262154439441526816992388562667004907159682643816214685929
63895217599993229915608941463976156518286253697920827223758251185210916864000000
000000000000000000
|#

 

GNU CLisp 에서는 네가지 타입의 부동소수점수

             short-float, single-float, double-float, long-float

들을 지원한다. clisp 를 대화형 모드로 실행하여 간단한 테스트를 해보자.

 

명령프롬프트> clisp
  i i i i i i i       ooooo    o        ooooooo   ooooo   ooooo
  I I I I I I I      8     8   8           8     8     o  8    8
  I  \ `+' /  I      8         8           8     8        8    8
   \  `-+-'  /       8         8           8      ooooo   8oooo
    `-__|__-'        8         8           8           8  8
        |            8     o   8           8     o     8  8
  ------+------       ooooo    8oooooo  ooo8ooo   ooooo   8

Welcome to GNU CLISP 2.49 (2010-07-07) <http://clisp.cons.org/>

Copyright (c) Bruno Haible, Michael Stoll 1992, 1993
Copyright (c) Bruno Haible, Marcus Daniels 1994-1997
Copyright (c) Bruno Haible, Pierpaolo Bernardi, Sam Steingold 1998
Copyright (c) Bruno Haible, Sam Steingold 1999-2000
Copyright (c) Sam Steingold, Bruno Haible 2001-2010

Type :h and hit Enter for context help.

[1]> pi
3.1415926535897932385L0
[2]> (/ 1 3)
1/3
[3]> (+ (/ 1 3) 0.0)
0.33333334
[4]> (coerce (/ 1 3) 'double-float)
0.3333333333333333d0
[5]> (coerce (/ 1 3) 'long-float)
0.33333333333333333334L0
[6]> (defun double-float(x)
    (coerce x 'double-float))
DOUBLE-FLOAT
[7]> (double-float (/ 1 3))
0.3333333333333333d0
[8]> (defun long-float(x)
    (coerce x 'long-float))
LONG-FLOAT
[9]> (long-float (/ 1 3))
0.33333333333333333334L0
[10]> (quit)
Bye.

 

 

#!/usr/bin/env clisp

;; Filename: testLanczos-02.lsp
;;
;;           An approximation for the gamma function by using the Lanczos algorithm
;;
;; Execute: clisp testLanczos-02.lsp
;;         or
;; Execute: ./testLanczos-02.lsp
;;
;; See: http://en.wikipedia.org/wiki/Lanczos_approximation
;; See:http://www-history.mcs.st-and.ac.uk/Biographies/Lanczos.html
;; See: http://www.thinbasic.com/community/showthread.php?11279-Gamma-function

 
;; Coefficients used by the GNU Scientific Library
(setf *g* 7)
(setf *p* (list
            0.99999999999980993L0   676.5203681218851L0       -1259.1392167224028L0
          771.32342877765313L0     -176.61502916214059L0         12.507343278686905L0
           -0.13857109526572012L0     9.9843695780195716L-6    1.5056327351493116L-7) )
 
(defun gamma(z)
    (let ((z (coerce z 'long-float))
          (tt (complex 0.0L0 0.0L0)))
        (if (realp z) (setf z (complex z 0.0L0)))
        ;; Reflection formula
        (if (< (realpart z) 0.5L0)
            (/ pi (sin (* (* pi z) (gamma (- 1L0 z)))))
            (progn
                (setf z (- z 1))
                (setf x (nth 0 *p*))
                (loop for i from 1 below (+ *g* 2) do
                    (setf x (+ x (/ (nth i *p*) (+ z i)))))
                (setf tt (+ z *g* 0.5L0))
                (* (sqrt (* 2L0 pi)) (expt tt (+ z 0.5L0)) (exp (- tt)) x) ))))

(defun factorial(n)
    (let ((k 0))
        (if (< n 2)
            1
            (progn
                (setf k 1)
                (if (evenp n)
                    (loop for i from 0 below (floor (/ n 2)) do
                        (setf k (* k (* (+ i 1) (- n i)))))
                    (progn
                        (loop for i from 0 below (floor (/ n 2)) do
                            (setf k (* k (* (+ i 1) (- n 1 i)))))
                        (setf k (* k n))))
                k))))

(defun facto(n)
    (let ((k 0))
        (if (< n 2)
            1
            (progn
                (setf k 1)
                (loop for i from 2 to n do
                    (setf k (* k i)) )
                k ))))

;; Begin here
(format t "gamma(10) = 9! = ~A asymptotically~%" (gamma 10))

(format t "gamma(21) = 20! = ~16A asymptotically~%" (gamma 21))
(format t "gamma(21) = 20! = ~16G asymptotically~%" (gamma 21))
(format t "gamma(21) = 20! = ~16G asymptotically~%" (realpart (gamma 21)))
(format t "gamma(21) = 20! = ~16F asymptotically~%" (gamma 21))
(format t "gamma(21) = 20! = ~16F asymptotically~%" (realpart (gamma 21)))
; (format t "gamma(101) = 100! = ~16A asymptotically~%" (gamma 101))
;; "{0.real:.15f}{0.imag:+.5f}j".format(gamma(101))

(ignore-errors
    (format t "gamma(101) = 100! = ~16A asymptotically~%" (gamma 101))
)

(format t "tgamma(~D) = ~D! = ~D~%"  101 100 (tgamma 101))

(setf i 10)
(format t "factorial(~D) = ~D! = ~D~%"  i i (factorial i))
(format t "facto(~D) = ~D! = ~D~%"  i i (facto i))

(setf i 20)
(format t "factorial(~D) = ~D! = ~D~%"  i i (factorial i))
(format t "facto(~D) = ~D! = ~D~%"  i i (facto i))

(setf i 100)
(format t "factorial(~D) = ~D! = ~D~%"  i i (factorial i))
(format t "facto(~D) = ~D! = ~D~%"  i i (facto i))

 

#|
Output:
gamma(10) = 9! = #C(362880.0000000007251L0 0.0L0) asymptotically
gamma(21) = 20! = #C(2.4329020081766424835L18 0.0L0) asymptotically
gamma(21) = 20! = #C(2.4329020081766424835L18 0.0L0) asymptotically
gamma(21) = 20! = 2432902008176642483.5     asymptotically
gamma(21) = 20! = #C(2.4329020081766424835L18 0.0L0) asymptotically
gamma(21) = 20! = 2432902008176642484. asymptotically
gamma(101) = 100! = #C(9.332621544393798272L157 0.0L0) asymptotically
tgamma(101) = 100! = 9.332621544394416d157
factorial(10) = 10! = 3628800
facto(10) = 10! = 3628800
factorial(20) = 20! = 2432902008176640000
facto(20) = 20! = 2432902008176640000
factorial(100) = 100! = 93326215443944152681699238856266700490715968264381621468
59296389521759999322991560894146397615651828625369792082722375825118521091686400
0000000000000000000000
facto(100) = 100! = 933262154439441526816992388562667004907159682643816214685929
63895217599993229915608941463976156518286253697920827223758251185210916864000000
000000000000000000
|#