swh:1:snp:81eadaa089e8253d8469bcef66aa332632c6c669
Tip revision: 113add4c7a9bfd37466a9b41cf2dc060b7b87091 authored by HwB on 06 December 2012, 00:00:00 UTC
version 1.3.1
version 1.3.1
Tip revision: 113add4
zeta.R
##
## z e t a . R eta and zeta Functions
##
eta <- function(z) {
if (!is.numeric(z) && !is.complex(z))
stop("Argument 'z' must be a real or comples vector.")
# find special cases
rez <- Re(z); imz <- Im(z)
i0 <- which(z == 0)
i2 <- which(z == (round(z/2)*2) & rez < 0 & imz == 0)
# i1 <- which(z == (round((z-1)/2)*2+1) & rez < 0 & imz == 0)
# reflection point
r <- 0.5
L <- which(rez < r)
if (length(L) > 0) {
zL <- z[L]
z[L] <- 1 - zL
}
# series coefficients are precalculated using the binomial distribution
cc <- c(
.99999999999999999997, -.99999999999999999821, .99999999999999994183, -.99999999999999875788,
.99999999999998040668, -.99999999999975652196, .99999999999751767484, -.99999999997864739190,
.99999999984183784058, -.99999999897537734890, .99999999412319859549, -.99999996986230482845,
.99999986068828287678, -.99999941559419338151, .99999776238757525623, -.99999214148507363026,
.99997457616475604912, -.99992394671207596228, .99978893483826239739, -.99945495809777621055,
.99868681159465798081, -.99704078337369034566, .99374872693175507536, -.98759401271422391785,
.97682326283354439220, -.95915923302922997013, .93198380256105393618, -.89273040299591077603,
.83945793215750220154, -.77148960729470505477, .68992761745934847866, -.59784149990330073143,
.50000000000000000000, -.40215850009669926857, .31007238254065152134, -.22851039270529494523,
.16054206784249779846, -.10726959700408922397,
.68016197438946063823e-1, -.40840766970770029873e-1, .23176737166455607805e-1, -.12405987285776082154e-1,
.62512730682449246388e-2, -.29592166263096543401e-2, .13131884053420191908e-2, -.54504190222378945440e-3,
.21106516173760261250e-3, -.76053287924037718971e-4, .25423835243950883896e-4, -.78585149263697370338e-5,
.22376124247437700378e-5, -.58440580661848562719e-6, .13931171712321674741e-6, -.30137695171547022183e-7,
.58768014045093054654e-8, -.10246226511017621219e-8, .15816215942184366772e-9, -.21352608103961806529e-10,
.24823251635643084345e-11, -.24347803504257137241e-12, .19593322190397666205e-13, -.12421162189080181548e-14,
.58167446553847312884e-16, -.17889335846010823161e-17, .27105054312137610850e-19)
nz <- length(z)
cc <- rev(cc)
ncc <- length(cc)
Z <- matrix(rep(z, each = ncc), nrow = ncc, ncol = nz)
N <- matrix(rep(seq(ncc, 1, by = -1), times = nz), nrow = ncc, ncol = nz)
# now compute the 'infinite' series
f <- drop(cc %*% N^-Z)
# and handle the special cases
if (length(L) > 0) {
zz <- z[L]
t <- (2-2^(zz+1)) / (2.^zz-2) / pi^zz
f[L] <- t * cos(pi/2*zz) * gammaz(zz) * f[L]
if (length(i0) > 0) f[i0] <- 0.5
if (length(i2) > 0) f[i2] <- 0.0
}
return(f)
}
zeta <- function(z) {
zz <- 2^z
k <- zz / (zz - 2)
f <- k * eta(z)
i1 <- which(z == 1)
if (length(i1) > 0) f[i1] <- Inf
return(f)
}