```正木です。

Riemann の ζ 関数です。

ζ(s)=Zeta[s]=ZetaA[s]+1/(s-1)

-----
N_s=Function.new(){"|s| Nat.Shift**s "}
BernoulliZ =((Bernoulli/Fact)*Double).Shift
Ms1s=Function.new{"|s|
Sequence.new([],s){\"|m,s| (s.zero?)? Log[m]:(m**s-1)/s \"}.Shift
"}
ZetaTerm0=Function.new{"|s| BernoulliZ*(Prods[s].Shift*Double) "}
ZetaA0=Function.new{"|s| N_s[-s].Sum0+(N_s[-s]/2)-Ms1s[1-s]"}
ZetaTerm=Function.new{"|s|
Sequence.new([],s){\"|m,s|
ZetaTerm0[s]/(Power[m]*Double).Shift*m**(1-s)\"}.Shift  "}
Zsec=Function.new{"|s| (s.kind_of?(Complex))? ((s.real>0)? s.imag.abs/s.real+1
: nil): 1"}
ZetaErrorFactor=Function.new{"|s| Sequence.new([],s){\"|k,s| Zsec[s+2*k+1] \"}
"}
ZetaErrorTerm=Function.new{"|s|
Sequence.new([],s){\"|m,s| ZetaErrorFactor[s]*ZetaTerm[s][m].abs \"}
"}
Zeta1=Function.new{"|s|
Sequence.new([],s){\"|m,s| ZetaTerm[s][m].Sum \"} "}
ZetaB=Function.new{"|s|
Sequence.new([],s){\"|k,s|
seq=Sequence.new([],k,s){'|m,k,s| Zeta1[s][m][k]+ZetaA0[s][m]'}
err=Sequence.new([],k,s){'|m,k,s| ZetaErrorTerm[s][m][k+1] '}
RealComplex.new(seq,err) \"} "}
ZetaA=Function.new{"|s|
seq=Sequence.new([],s){'|k,s| ZetaB[s][k] '}
err=Sequence.new([],s){'|k,s| ZetaErrorTerm[s][k][k] '}
if (s.integer? || s==s.to_i) && s<=0
RealComplex.new(seq,err).approx(0)
else
RealComplex.new(seq,err)
end
"}
Zeta=Function.new{"|s| ZetaA[s]+1/(s-1) "}
EulerConstant=ZetaA[1]
#0.5772156649015328606065120900824024310421

def Math.Zeta(s,e=2**-53)
ZetaA[s].approx(e)+1/(s-1)
end

class Sequence2
def approx(eps)
z=self[find_error{|e,i| e <= eps }[1]]
return z.approx(eps) if z.type==Sequence2 || z.type==Complex
z
end
end

class Rational
def to_ContFraction
a=@numerator
b=@denominator
s=Sequence.new([a,b]){"|n| self[n-2]%self[n-1] rescue nil "}
Sequence.new([],s){"|n,s| s[n].div(s[n+1]) rescue nil"}
end
def to_Real
Real.new(to_ContFraction,"CF")
end
end

Root=Function.new{"|x|
Sequence.new([],x){'|m,x| a=1
s=Sequence.new([[a,1,1]],m,x){\"|n,m,x|
a,e,f=self[n-1]
a1=a.to_Real.approx(e**2/f)
e1=x/a1**(m-1)-a1
a2=a1+e1/m
e2=e1**2
f*=m if e<=e2
[a2,e2,f]
\"}.Shift
RealComplex.new(s,\"W\")
'}"}

class Integer,Rational
def **(y)
x=self
case y
when Integer
if y >= 0
Power[self][y]
elsif y < 0
1/Power[self][-y]
end
when Rational
return Sqrt[x] if y==1/2
return Root[x**y.numerator][y.denominator] if 0<y && y<1
n=y.floor
x**(y-n)*x**n
when Complex
return Float(self) ** y if y.real.type==Float
Exp[Log[x]*y]
when Float
Float(self) ** y
else
a , b = y.coerce(self)
a ** b
end
end

-----

Math.Zeta の引数と戻り値は
s                     Math.Zeta(s)
0 又は負の整数        Rational 又は 0

Float                 Float
Complex(Float,Float)  Complex(Float,Float)

(40桁計算するのに(Pentium 150MHzで)30秒ほどかかります。)

```