Visual Studio 2019 와 MSYS2 MinGW64 에서 테스트 된 소스입니다.

 혹시 MinGW 에서 컴파일되지 않으면

$ packman -S mpfr

명령으로 mpfr 라이브러리를 설치하고 컴파일하면 된다.

 

// Filename: calcGammaFn.c
//
// Compile: gcc -o calcGammaFn calcGammaFn.c -lmpfr -lgmp
// Execute: ./calcGammaFn
//     Or
// Compile: cl calcGammaFn.c /I. mpfr.lib
// Execute: calcGammaFn
//
// Date: 2021.01.28


#include <stdio.h>
#include <math.h>    // for log(10)
#include <mpfr.h>

int main()
{
  mpfr_t x; 
  int i, inex;

  mpfr_set_emin (-41);
  mpfr_init2 (x, 42);

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_gamma (x, x, MPFR_RNDZ);      
      mpfr_subnormalize (x, inex, MPFR_RNDZ);
      mpfr_dump (x);
  }

  printf("\n");

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_gamma (x, x, MPFR_RNDZ);
      mpfr_printf("Gamma(%2d) = %2d! = %Rf\n", i, i - 1, x);
  }

  printf("\n");

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_lngamma (x, x, MPFR_RNDZ);
      mpfr_printf("LogGamma(%2d) = Log(%2d!) = %Rf\n", i, i - 1, x);
  }

  printf("\n");

 double t10 = log(10.0);
 printf("log(10) = %f\n", t10);
 
  printf("\n");

  for (i = 1; i <= 17; i++)
  {
      mpfr_set_ui (x, i, MPFR_RNDN);
      inex = mpfr_lngamma (x, x, MPFR_RNDZ);      
      inex = mpfr_div_d(x, x, t10,  MPFR_RNDZ);
      mpfr_printf("Log10Gamma(%2d) = Log10(%2d!) = %Rf\n", i, i - 1, x);
  }

  mpfr_clear (x);
  
  return 0;
}

/*
Output:
LogGamma( 1) = Log( 0!) = 0
LogGamma( 2) = Log( 1!) = 0
LogGamma( 3) = Log( 2!) = 0.693147
LogGamma( 4) = Log( 3!) = 1.791759
LogGamma( 5) = Log( 4!) = 3.178054
LogGamma( 6) = Log( 5!) = 4.787492
LogGamma( 7) = Log( 6!) = 6.579251
LogGamma( 8) = Log( 7!) = 8.525161
LogGamma( 9) = Log( 8!) = 10.604603
LogGamma(10) = Log( 9!) = 12.801827
LogGamma(11) = Log(10!) = 15.104413
LogGamma(12) = Log(11!) = 17.502308
LogGamma(13) = Log(12!) = 19.987214
LogGamma(14) = Log(13!) = 22.552164
LogGamma(15) = Log(14!) = 25.191221
LogGamma(16) = Log(15!) = 27.899271
LogGamma(17) = Log(16!) = 30.671860

log(10) = 2.302585

Log10Gamma( 1) = Log10( 0!) = 0
Log10Gamma( 2) = Log10( 1!) = 0
Log10Gamma( 3) = Log10( 2!) = 0.301030
Log10Gamma( 4) = Log10( 3!) = 0.778151
Log10Gamma( 5) = Log10( 4!) = 1.380211
Log10Gamma( 6) = Log10( 5!) = 2.079181
Log10Gamma( 7) = Log10( 6!) = 2.857332
Log10Gamma( 8) = Log10( 7!) = 3.702431
Log10Gamma( 9) = Log10( 8!) = 4.605521
Log10Gamma(10) = Log10( 9!) = 5.559763
Log10Gamma(11) = Log10(10!) = 6.559763
Log10Gamma(12) = Log10(11!) = 7.601156
Log10Gamma(13) = Log10(12!) = 8.680337
Log10Gamma(14) = Log10(13!) = 9.794280
Log10Gamma(15) = Log10(14!) = 10.940408
Log10Gamma(16) = Log10(15!) = 12.116500
Log10Gamma(17) = Log10(16!) = 13.320620
*/

 

 

 

 

 

Posted by Scripter

댓글을 달아 주세요

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

댓글을 달아 주세요