# Gamma and Zeta function for complex numbers # (C) Minato Nakazawa # rev. 1.02 on 9 December 2009. # rev. 1.03 on 19 December 2009. Bug fix of Zeta2() and increase precision. # # Note: Gamma function was translated from C source by # (C) Takuya OOURA , 1996. # http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf-j.html Gammac <- function(x) { xr <- Re(x) xi <- Im(x) if (xr < 0) { wr <- 1 - xr wi <- -xi } else { wr <- xr wi <- xi } ur <- wr + 6.00009857740312429 vr <- ur * (wr + 4.99999857982434025) - wi * wi vi <- wi * (wr + 4.99999857982434025) + ur * wi yr <- ur * 13.2280130755055088 + vr * 66.2756400966213521 + 0.293729529320536228 yi <- wi * 13.2280130755055088 + vi * 66.2756400966213521 ur <- vr * (wr + 4.00000003016801681) - vi * wi ui <- vi * (wr + 4.00000003016801681) + vr * wi vr <- ur * (wr + 2.99999999944915534) - ui * wi vi <- ui * (wr + 2.99999999944915534) + ur * wi yr <- yr + ur * 91.1395751189899762 + vr * 47.3821439163096063 yi <- yr + ui * 91.1395751189899762 + vi * 47.3821439163096063 ur <- vr * (wr + 2.00000000000603851) - vi * wi ui <- vi * (wr + 2.00000000000603851) + vr * wi vr <- ur * (wr + 0.999999999999975753) - ui * wi vi <- ui * (wr + 0.999999999999975753) + ur * wi yr <- yr + ur * 10.5400280458730808 + vr yi <- yi + ui * 10.5400280458730808 + vi ur <- vr * wr - vi * wi ui <- vi * wr + vr * wi t <- ur * ur + ui * ui vr <- yr * ur + yi * ui + t * 0.0327673720261526849 vi <- yi * ur - yr * ui yr <- wr + 7.31790632447016203 ur <- log(yr * yr + wi * wi) * 0.5 - 1 ui <- atan2(wi, yr) yr <- exp(ur * (wr - 0.5) - ui * wi - 3.48064577727581257) / t yi <- ui * (wr - 0.5) + ur * wi ur <- yr * cos(yi) ui <- yr * sin(yi) yr <- ur * vr - ui * vi yi <- ui * vr + ur * vi if (xr < 0) { wr <- xr * pi wi <- exp(xi * pi) vi <- 1 / wi ur <- (vi + wi) * sin(wr) ui <- (vi - wi) * cos(wr) vr <- ur * yr + ui * yi vi <- ui * yr - ur * yi ur <- 2*pi / (vr * vr + vi * vi) yr <- ur * vr yi <- ur * vi } (yr + yi*1i) } # Zeta <- function(s) { if (s==1) { ret <- Inf } else if (Re(s)>0) { a <- 1/1^s for (i in 2:10000) a <- a+(-1)^(i-1)/i^s ret <- 1/(1-2^(1-s))*a } else { t <- 1-s ret <- 2*(2*pi)^(-t)*cos(t*pi/2)*Gammac(t)*Zeta(t) } ret } # require(gsl) Zeta2 <- function(s) { if (s==1) { ret <- Inf } else if (Re(s)>0) { a <- 0 for (i in 100000:1) a <- a+(-1)^(i-1)/i^s ret <- 1/(1-2^(1-s))*a } else { t <- 1-s ret <- 2*(2*pi)^(-t)*cos(t*pi/2)*exp(lngamma_complex(t))*Zeta2(t) } ret } Zeta(-1) Zeta(1) Zeta(2) Zeta(-1+1i*5) Zeta(1/2+14.1347*1i) # Zeta2(-1) Zeta2(1) Zeta2(2) Zeta2(-1+1i*5) Zeta2(1/2+14.1347*1i)