Lanczos 알고리즘은 Stirlng 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다. 단지 exp 함수를 이용하는 부분에서는 exp 함수의 구현에 따라 오차가 더 있을 수 있다. 단, Common Lisp 언어의 float 타입은 유효숫자를 약 7개~8개 밖에 알아내지 못한다.
Comon Lisp 언어의 디퐁트 부동소수점수는 단정밀도 부동소수점수이다. 단정밀도 부동소수점수는 그 처리 가능한 양수의 범위가 최소 2−126 ≈ 1.18 × 10−38 부터 최대 (2−2−23) × 2127 ≈ 3.4 × 1038 밖에 되지 않는다. 이 범위를 벗어나는 양수에 대해서는 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.
;; 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
|#
'프로그래밍 > Common Lisp' 카테고리의 다른 글
| 스트링 벡터에서 스트링 찾기(find) with Common Lisp (0) | 2013.09.12 |
|---|---|
| Common Lisp 언어로 역삼각함수, 역쌍곡선함수 값을 구하는 예제 (0) | 2013.09.11 |
| Common Lisp 언어로 복소수 계산하기 (0) | 2013.09.07 |
| Common Lisp 를 이용한 간단한 수학식 계산하기 (0) | 2013.09.07 |
| 조립제법(Horner의 방법) 예제 2 for Common Lisp (0) | 2013.09.07 |