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

댓글을 달아 주세요

아래의 소스는 윈도우에서

            Luna MinGW & GNU C 4.5.0 (gcc), 

로 테스트되었다. long 타입으로는 13! 까지만 정확하계 계산되지만 GMP 를 이용한 계산은 아무리 큰 수의 곱셈이라도 정확히 계산해준다.

  • 윈도우에 Luna MinGW (with GCC 4.5.0) 설치하기:
     1) Luna MinGW 홈페이지(MinGW GCC C/C++ Compiler package with installer) 
     2) Luna MinGW 설치파일 다운로드
  •  영문  위키피디아에서 소개하는 MinGW
  • MinGW 의 공식 웹사이트에서 MinGW 를 설치하면 gcc 버전이 3.4.5 밖에 되지 않고,
    gmp 라이브러리도 수동으로 설치헤야 하는 번거로움이 있으므로,
    여기서는 Luna MinGW 를 설치하기로 한다. Luna MinGW 를 설치하면,
    현재 최신 버전인 gcc4.5.0 을 쓸 수 있다. 그 대신 g++ 로 컴파일 후에 쓰일 런타임 다이나믹  라이브러리 libgmpxx-4.dll 은 별도로 구하여 설치해야 한다.
  • 윈도우에 Luna MinGW (with GCC 4.5.0) 설치하기:
        * Luna MinGW 홈페이지(MinGW GCC C/C++ Compiler package with installer) 
        * Luna MinGW 설치파일 다운로드
    1) 설치는 C:\MinGW 에 한것으로 간주한다.
    2) gcc 에서 gmp 라이브러리는 추가 설치 없이 바로 쓸 수 있는데,
        g++ 로는 컴파일은 되지만 컴파일된 실행파일을 실행하면 
        limgmpxx-4.dll 파일이 없다고 에러가 날 것이다. 이때는 아래 처럼 
        limgmpxx-4.dll 파일을 컴파일해서 만들던가 이미 만들어진 것을 PATH 가 걸려 있는 곳에 넣어두어야 한다. 아래는 limgmpxx-4.dll 파일을 구할 수 있는 곳이다. 
  • g++에서 gmp를 사용하기 위해 libgmpxx-4.dll 구하여 설치하기:
       * 소스포지에서 libgmpxx-5.0.1-1-mingw32-dll-4.tar.lzma 다운로드하기
    1)  확장명이 lzma 인 파일은 7zip 으로 풀어야 한다. (7zip 은 http://www.7-zip.org 에서 구한다.)
    2) 위의 파일을 받아서 7zip으로 압축 해제하여 libgmpxx-4.dll 을  C:\MinGW\bin 에 복사힌다.  (C:\MinGW 은 Luna MinGW 가 설치된 디렉토리이다)
  • gcc 버전 확인하기
    프롬프트> gcc -v
    Using built-in specs.
    COLLECT_GCC=gcc
    COLLECT_LTO_WRAPPER=c:/mingw/bin/../libexec/gcc/mingw32/4.5.0/lto-wrapper.exe
    Target: mingw32
    Configured with: ../gcc-4.5.0/configure --enable-languages=c,c++,ada,fortran,obj
    c,obj-c++ --disable-sjlj-exceptions --with-dwarf2 --enable-shared --enable-libgo
    mp --disable-win32-registry --enable-libstdcxx-debug --enable-version-specific-r
    untime-libs --disable-werror --build=mingw32 --prefix=/mingw
    Thread model: win32
    gcc version 4.5.0 (GCC)
  • g++ 버전 확인하기
    프롬프트> g++ -v
    Using built-in specs.
    COLLECT_GCC=g++
    COLLECT_LTO_WRAPPER=c:/mingw/bin/../libexec/gcc/mingw32/4.5.0/lto-wrapper.exe
    Target: mingw32
    Configured with: ../gcc-4.5.0/configure --enable-languages=c,c++,ada,fortran,obj
    c,obj-c++ --disable-sjlj-exceptions --with-dwarf2 --enable-shared --enable-libgo
    mp --disable-win32-registry --enable-libstdcxx-debug --enable-version-specific-r
    untime-libs --disable-werror --build=mingw32 --prefix=/mingw
    Thread model: win32
    gcc version 4.5.0 (GCC)



    * 소스 파일명: recFactGMPcpp01.cc

/*
 *  Filename: recFactGMPcpp01.cc
 *
 *  Compile: g++ -o recFactGMPcpp01 recFactGMPcpp01.cc -lgmpxx -lgmp
 */

#include <iostream>
#include <gmpxx.h>

using namespace std;

mpz_class factorial(int n) {
    static mpz_class v;
    v = 1;
    for (int i = 1; i <= n; i++) {
     v = v*i;
    }
    return v;
}

int main(int argc, char *argv[]) {
    int n1 = 9;
    int n2 = 30;
    for (int i = n1; i <= n2; i++) {
        cout << i << "! = " << factorial(i) << endl;
        if (i == 13)
            cout << "-- below calc are regular ----" << endl;
    }
    return 0;
}




실행 결과:
프롬프트> recFactGMPcpp01
9! = 362880
10! = 3628800
11! = 39916800
12! = 479001600
13! = 6227020800
-- below calc are regular ----
14! = 87178291200
15! = 1307674368000
16! = 20922789888000
17! = 355687428096000
18! = 6402373705728000
19! = 121645100408832000
20! = 2432902008176640000
21! = 51090942171709440000
22! = 1124000727777607680000
23! = 25852016738884976640000
24! = 620448401733239439360000
25! = 15511210043330985984000000
26! = 403291461126605635584000000
27! = 10888869450418352160768000000
28! = 304888344611713860501504000000
29! = 8841761993739701954543616000000
30! = 265252859812191058636308480000000



Posted by Scripter

댓글을 달아 주세요

아래의 소스는 윈도우에서

            Luna MinGW & GNU C 4.5.0 (gcc), 

로 테스트되었다. long 타입으로는 13! 까지만 정확하계 계산되지만 GMP 를 이용한 계산은 아무리 큰 수의 곱셈이라도 정확히 계산해준다.


* 소스 파일명: recFactGMP01.c

/*
 *  Filename: recFactGMP01.c
 *
 *  Compile: gcc -o recFactGMP01 recFactGMP01.c -lgmp
 */

#include <stdio.h>
#include <gmp.h>

void factorial(mpz_t v, int n) {
 int i;
    mpz_t t;
    mpz_init (t);
    mpz_init_set_str (v, "1", 0);
    for (i = 1; i <= n; i++) {
        mpz_init_set_si (t, i);
        mpz_mul(v, v, t);
    }
}

int main(int argc, char *argv[]) {
    int n1 = 9;
    int n2 = 30;
    int i;
    mpz_t v;
    mpz_init (v);
    for (i = n1; i <= n2; i++) {
        mpz_init_set_si (v, 1);
        factorial(v, i);
        gmp_printf("%d! = %Zd\n", i, v);
        if (i == 13)
        printf("-- below calc are regular ----\n");
    }
    return 0;
}



실행 결과:
프롬프트> recFactGMP01
9! = 362880
10! = 3628800
11! = 39916800
12! = 479001600
13! = 6227020800
-- below calc are regular ----
14! = 87178291200
15! = 1307674368000
16! = 20922789888000
17! = 355687428096000
18! = 6402373705728000
19! = 121645100408832000
20! = 2432902008176640000
21! = 51090942171709440000
22! = 1124000727777607680000
23! = 25852016738884976640000
24! = 620448401733239439360000
25! = 15511210043330985984000000
26! = 403291461126605635584000000
27! = 10888869450418352160768000000
28! = 304888344611713860501504000000
29! = 8841761993739701954543616000000
30! = 265252859812191058636308480000000



Posted by Scripter

댓글을 달아 주세요

  1. 진우 2012.09.21 16:47  댓글주소  수정/삭제  댓글쓰기

    글 잘봤습니다
    저도 gmp 라이브러리에 관심이 많은데요
    혹시 윈도우 visual studio 2005나 2008 2010 등에서
    gmp 라이브러리를 사용하는 방법 아시나요?

  2. Scripter 2013.01.23 01:14 신고  댓글주소  수정/삭제  댓글쓰기

    Visual Studio 2008의 경우에는 http://scripting.tistory.com/65 를 참고하십시오.
    Visual Studio 용으로는 gmp 라이브러리 보다 mpir 라이브러리 추천합니다. http://scripting.tistory.com/626 을 참고하십시오.


* 꼬리 재귀호출과 패턴 매칭을 이용하여 구현한 팩토리얼과 피보나치 수열 계산

{-
    Filename: fact.hs
         Rapid factorial and fibonacci seq implementations
         by pattern matching and tail recursive call

    Compile: ghc fact.hs
    Execute: main

    Date: 2010/07/20
    Author: phkim  pkim (AT) scripts.pe.kr
-}

module Main where

  factorial :: Integer -> Integer
  factorial n =  recfact 1 n

  recfact :: Integer -> Integer -> Integer
  recfact acc k = case k of
        0 -> acc
        k-> recfact (acc*k) (k - 1)

  fib :: Integer -> Integer
  fib n = fibGen 0 1 n

  fibGen :: Integer -> Integer -> Integer -> Integer
  fibGen a b n = case n of
        0 -> a
        n -> fibGen b (a + b) (n - 1)

  main :: IO ()
  main = do
        print ("Factorial(30000) has " ++ show (length (show (factorial 30000))) ++ " digits")
        print ("Fibonacci(200000) has " ++ show (length (show (fib 200000))) ++ " digits")

{-
    Expected result:
    "Factorial(30000) has 121288 digits"
    "Fibonacci(200000) has 41798 digits"
-}



크리에이티브 커먼즈 라이선스
Creative Commons License

 

Posted by Scripter

댓글을 달아 주세요