Lanczos 알고리즘은 Stirlng 공식에 의한 알고리즘 보다 정밀하며, 십진수로 유효숫자 약 15자리 까지는 정확하게 계산해 준다.  단지 exp 함수를 이용하는 부분에서는 exp 함수의 구현에 따라 오차가 더 있을 수 있다.


#!/usr/bin/env python
# -*- encoding:utf-8 -*-

# Filename: testLanczos-01.py
#
#           An approximation for the gamma function by using the Lanczos algorithm
#
# Execute: python testLanczos-01.py
#         or
# Execute: ./testLanczos-01.py
#
# 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

from cmath import *
 
# Coefficients used by the GNU Scientific Library
g = 7
p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
     771.32342877765313, -176.61502916214059, 12.507343278686905,
     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
 
def gamma(z):
    z = complex(z)
    # Reflection formula
    if z.real < 0.5:
        return pi / (sin(pi*z)*gamma(1-z))
    else:
        z -= 1
        x = p[0]
        for i in range(1, g+2):
            x += p[i]/(z+i)
        t = z + g + 0.5
        return sqrt(2*pi) * t**(z+0.5) * exp(-t) * x

def factorial(n):
    if n < 2:
        return 1
    else:
        k = 1
        if n % 2 == 0:
            for i in xrange(n/2):
                k *= (i + 1) * (n - i)
        else:
            for i in xrange(n/2):
                k *= (i + 1) * (n - 1 - i)
            k *= n
        return k

def facto(n):
    if n < 2:
        return 1
    else:
        k = 1
        for i in xrange(2, n + 1):
                k *= i
        return k

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! = (362880+0j) asymptotically
gamma(101) = 100! = 9.33262154439379e+157+0.00000j asymptotically
0! = 1
1! = 1
2! = 2
3! = 6
4! = 24
5! = 120
6! = 720
7! = 5040
8! = 40320
9! = 362880
10! = 3628800
factorial(100) = 100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
facto(100) = 100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
"""



Posted by Scripter

댓글을 달아 주세요


감마함수 Γ(x)의 정의

(i) \textrm{$\alpha > 0$일 때는}

         \Gamma (\alpha) = \int_0^\infty e^{-t} t^{\alpha - 1} \ dt

(ii) \textrm{$\alpha < 0$ 이고 $\alpha \ne$정수 일 때는}

   \textrm{$\alpha + k> 0$ 이 되는 최소의 양의 정수 $k$를 찾아서}

         \Gamma (\alpha) = \dfrac{\Gamma(\alpha + k)}{\alpha (\alpha + 1) \cdots (\alpha +k - 1)}


[감마함수 Γ(x)의 특징]

(1) \Gamma (\alpha + 1) = \alpha \cdot \Gamma(\alpha)

(2) \textrm{특히 $n$이 양의 정수일 때는 \ } \Gamma (n) = (n - 1)!



아래에서는 Ubuntu에서 gnuplot을 이용하여 감마함수 Γ(x)의 그래프를 단계적으로 완성해 나가는 과정을 보여준다. 과정은 모두 5단계로 나뉘어져 있다. 일단 그려주고 나서 옵션을 조금씩 변경 또는 추가 하면서 그림을 다시 그려주는 방식을 취한다. (그려진 그림을 다시 그리는 gnuplot 명령은 replot이다)


* 단계 1: 감마함수의 그래프를 일단 그린다. (옵션은 차후에 조절한다.)

$ gnuplot

    G N U P L O T
    Version 4.2 patchlevel 6
    last modified Sep 2009
    System: Linux 2.6.32-33-generic

    Copyright (C) 1986 - 1993, 1998, 2004, 2007 - 2009
    Thomas Williams, Colin Kelley and many others

    Type `help` to access the on-line reference manual.
    The gnuplot FAQ is available from http://www.gnuplot.info/faq/

    Send bug reports and suggestions to <http://sourceforge.net/projects/gnuplot>


Terminal type set to 'wxt'
gnuplot> plot [-5:5] [-5:5] gamma(x)




* 단계 2: x좌표축, y좌표축에 레이블을 붙여주고, 그래프의 제목을 표시한다.

gnuplot> set xlabel "x values"
gnuplot> set ylabel "y values"
gnuplot> set title "Gamma function"
gnuplot> replot




* 단계 3: 좌표를 읽기 쉽게 정수 간격으로 그리드를 만든다. (set grid 명령은 그리드를 ON 하는 명령)

gnuplot> set xtics 1
gnuplot> set ytics 1
gnuplot> set grid
gnuplot> replot




* 단계 4: x축과 y축을 (0.8 두께의 청색 실선으로) 그려준다.

gnuplot> set xzeroaxis linetype 3 linewidth 0.8
gnuplot> set yzeroaxis linetype 3 linewidth 0.8
gnuplot> replot





* 단계 5: 수직 방향의 점근선 부근에 잘못된 부분 처리하기  ('set samples 숫자' 명령 사용)

gnuplot> set samples 1600
gnuplot> replot





Posted by Scripter

댓글을 달아 주세요


감마함수 Γ(x)의 정의

(i) \textrm{$\alpha > 0$일 때는}

         \Gamma (\alpha) = \int_0^\infty e^{-t} t^{\alpha - 1} \ dt

(ii) \textrm{$\alpha < 0$ 이고 $\alpha \ne$정수 일 때는}

   \textrm{$\alpha + k> 0$ 이 되는 최소의 양의 정수 $k$를 찾아서}

         \Gamma (\alpha) = \dfrac{\Gamma(\alpha + k)}{\alpha (\alpha + 1) \cdots (\alpha +k - 1)}


[감마함수 Γ(x)의 특징]

(1) \Gamma (\alpha + 1) = \alpha \cdot \Gamma(\alpha)

(2) \textrm{특히 $n$이 양의 정수일 때는 \ } \Gamma (n) = (n - 1)!



* Ubuntu의 KmPlot을 이용하여 감마함수 Γ(x)의 그래프를 그린 화면

    (함수식에 f(x) = gamma(x) 라고 적어주고, 'Create" 버튼을 클릭하면 그려진다.)



* KmPlot 의 메뉴에서 "File" -> "Export..." 를 선택하여 출력한 감마함수 Γ(x)의 그래프




Posted by Scripter

댓글을 달아 주세요