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

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

#!/usr/bin/env julia

# Filename: testLanczos-01.jl
#
#           An approximation for the gamma function by using the Lanczos algorithm
#
# Execute: julia testLanczos-01.jl
#         or
# Execute: ./testLanczos-01.jl
#
# 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
g = 7
p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,
     771.32342877765313, -176.61502916214059, 12.507343278686905,
     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7]
 
function gamma(z)
    z = complex(z)
    # Reflection formula
    # if z.real < 0.5
    if real(z) < 0.5
        return pi / (sin(pi*z)*gamma(1-z))
    else
        z -= 1
        x = p[1]
        for i =1:g+1
            x += p[i+1]/(z+i)
        end
        t = z + g + 0.5
        return sqrt(2*pi) * t^(z+0.5) * exp(-t) * x
    end
end

function factorial(n)
    if n < 2
        return 1
    else
        k = 1
        if n % 2 == 0
            for i = 0:div(n, 2)-1
                k *= (i + 1) * (n - i)
            end
        else
            for i = 0:div(n, 2)-1
                k *= (i + 1) * (n - 1 - i)
            end
            k *= n
        end
        return k
    end
end

function facto(n)
    if n < 2
        return BigInt(1)
    else
        k = BigInt(1)
        for i = 2:n
            k *= BigInt(i)
        end
        return k
    end
end

@printf("gamma(10) = 9! = %s asymptotically\n", real(gamma(10)))
@printf("gamma(101) = 100! = %s asymptotically\n", real(gamma(101)))
for i = 0:10
    @printf("%d! = %d\n", i, factorial(i))
end
  
i = 100
@printf("factorial(%d) = %d! = %s\n", i, i, factorial(BigInt(i)))
@printf("facto(%d) = %d! = %s\n", i, i, facto(BigInt(i)))

"""
Output:
gamma(10) = 9! = 362880.0000000015 asymptotically
gamma(101) = 100! = 9.332621544393794e157 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! = 93326215443944152681699238856266700490715968264381621468
59296389521759999322991560894146397615651828625369792082722375825118521091686400
0000000000000000000000
facto(100) = 100! = 933262154439441526816992388562667004907159682643816214685929
63895217599993229915608941463976156518286253697920827223758251185210916864000000
000000000000000000
"""