54 lines
1.4 KiB
R
54 lines
1.4 KiB
R
stirling <- function(z) sqrt(2*pi/z) * (exp(-1)*z)^z
|
|
|
|
nemes <- function(z) sqrt(2*pi/z) * (exp(-1)*(z + (12*z - (10*z)^-1)^-1))^z
|
|
|
|
lanczos <- function(z)
|
|
{
|
|
if(length(z) > 1)
|
|
{
|
|
sapply(z, lanczos)
|
|
} else
|
|
{
|
|
g <- 7
|
|
p <- c(0.99999999999980993, 676.5203681218851, -1259.1392167224028,
|
|
771.32342877765313, -176.61502916214059, 12.507343278686905,
|
|
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7)
|
|
z <- as.complex(z)
|
|
if(Re(z) < 0.5)
|
|
{
|
|
pi / (sin(pi*z) * lanczos(1-z))
|
|
} else
|
|
{
|
|
z <- z - 1
|
|
x <- p[1]
|
|
for (i in 1:8) {
|
|
x <- x+p[i+1]/(z+i)
|
|
}
|
|
tt <- z + g + 0.5
|
|
sqrt(2*pi) * tt^(z+0.5) * exp(-tt) * x
|
|
}
|
|
}
|
|
}
|
|
|
|
spouge <- function(z, a=49)
|
|
{
|
|
if(length(z) > 1)
|
|
{
|
|
sapply(z, spouge)
|
|
} else
|
|
{
|
|
z <- z-1
|
|
k <- seq.int(1, a-1)
|
|
ck <- rep(c(1, -1), len=a-1) / factorial(k-1) * (a-k)^(k-0.5) * exp(a-k)
|
|
(z + a)^(z+0.5) * exp(-z - a) * (sqrt(2*pi) + sum(ck/(z+k)))
|
|
}
|
|
}
|
|
|
|
# Checks
|
|
z <- (1:10)/3
|
|
all.equal(gamma(z), stirling(z)) # Mean relative difference: 0.07181942
|
|
all.equal(gamma(z), nemes(z)) # Mean relative difference: 0.003460549
|
|
all.equal(as.complex(gamma(z)), lanczos(z)) # TRUE
|
|
all.equal(gamma(z), spouge(z)) # TRUE
|
|
data.frame(z=z, stirling=stirling(z), nemes=nemes(z), lanczos=lanczos(z), spouge=spouge(z), builtin=gamma(z))
|