Rで Hurwitz's zeta function , Hurwitz zeta function フルビッツゼータ関数の近似

RでHurwitz's zeta関数 フルビッツゼータ関数の近似
よいライブラリがないので、とりあえず、つくってみた。

library("Rmpfr")

Hurwit_zeta_approx<-function(beta4,x){
	#N<-10^7
	#k<-0:N
	#s<-beta4
	ans2<-NULL
	if(beta4==beta4){
		for(j in 1:length(x)){
			a3<-x[j]
			#sum(1/(k+a)^s)
			#(a3)^(1-beta4)/(beta4-1)+0.5*(a3)^(-beta4)
			n<-30
			n2<-30
			s<-beta4
			k<-1:n
			kk<-0:n
			B1<-(as.numeric(Bernoulli(2*(1:n2))))
			B2<-lgamma(2*(1:n2)+1)
			#B2<-gamma(2*(1:n2)+1)

			B3<-NULL
			for(i in 1:n2){
				B3[i]<-log(s)+sum(log(s+1:(2*i)))
				#B3[i]<-(s)*prod((s+1:(2*i)))

			}	
			B4<-(s+2*(1:n2)+1)*log(a3+n)
			#B4<-(a3+n)^(s+2*(1:n2)+1)
			ans2[j]<-sum((a3+kk)^-s)+(a3+n)^(1-s)/(s-1)-1/2/(a3+n)^s+sum(B1*exp(-B2+B3-B4))[[1]]
			#sum((a3+kk)^-s)+(a3+n)^(1-s)/(s-1)-1/2/(a3+n)^s+sum(B1/B2*B3/B4)
		}
	
	}
	as.vector(ans2)
}