正木です。

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
正の整数,Rational     Realの近似値
Float                 Float
Complex(Float,Float)  Complex(Float,Float)

上記の Euler Constant の計算は Γ 関数の対数微分を使う方法と
実質的に同じ algorithm です。
(40桁計算するのに(Pentium 150MHzで)30秒ほどかかります。)