正木です。

Γ 関数と Bessel 関数です。

#  DLogGamma[s]=Γ'(s)/ Γ(s) 
EvenPower=Function.new{"|x| Power[x]*Double "}
DMuTerm=Function.new{"|s| (Bernoulli*Double/Even/EvenPower[s]).Shift "}
#DMuTerm[s][n]=Bernoulli[2*n+2]/(2*n+2)/s**(2*n+2)
DMu=Function.new(){"|s| -DMuTerm[s].Sum "}
InverseSum=Function.new{"|s| (Nat+s).inverse.Sum0 "}
DLogGammaError=Function.new{"|s| 
  Sequence.new([],s){'|k,s| 
    Sequence.new([],k,s){\"|n,k,s| DMuTerm[s+n][k].abs1 \"} '} "}
DLogGammaM=Function.new{"|s| 
  Sequence.new([],s){'|k,s| 
    seq=Sequence.new([],k,s){\"|n,k,s| t=s+n
      DMu[t][k]+Log[t]-1/(2*t)-InverseSum[s][n] \"}
    RealComplex.new(seq,DLogGammaError[s][k]) '} "}
DLogGamma=Function.new{"|s| 
  if (s.kind_of?(Integer) || s.to_i==s) && s<=0
    nil
  elsif s.real<1
    DLogGamma[s+1]-1/s
  else
    seq=Sequence.new([],s){'|k,s| DLogGammaM[s][k] '}
    err=Sequence.new([],s){'|k,s| DLogGammaError[s][k][k] '}
    RealComplex.new(seq,err)
  end
  "}
EulerConstant=-DLogGamma[1]

#  LogGamma[s]=log(Γ(s))
MuTerm=Function.new{"|s| (Bernoulli*Double/Even).Shift/Odd/OddPower[s] "}
#MuTerm[s][n]=Bernoulli[2*n+2]/(2*n+2)/(2*n+1)/s**(2*n+1)
GammaMu=Function.new{"|s| MuTerm[s].Sum "}
LogProd=Function.new{"|s| Sequence.new([0],s){\"|n,s| Log[s+n-1] \"}.Sum "}
GammaMuError=Function.new{"|s| 
  Sequence.new([],s){'|k,s| 
    Sequence.new([],k,s){\"|n,k,s| MuTerm[s+n][k].abs1 \"}'}"}
LogStirling0=Function.new{"|s| (s-1/2)*Log[s]-s "}
LogGamma0J=Function.new{"|s| 
    Sequence.new([],s){\"|k,s| LogStirling0[s]+GammaMu[s][k]\"} "}
LogGamma0N=Function.new{"|s| 
  Sequence.new([],s){'|k,s| 
    Sequence.new([],k,s){\"|n,k,s| LogGamma0J[s+n][k]-LogProd[s][n]\"}'}"}
LogGamma0M=Function.new{"|s| 
  Sequence.new([],s){'|k,s| 
    RealComplex.new(LogGamma0N[s][k],GammaMuError[s][k+1]) '} "}
LogGamma0=Function.new{"|s| 
  err=Sequence.new([],s){'|k,s| MuTerm[s+k][k+1].abs1 '} 
  RealComplex.new(LogGamma0M[s],err)
  "}
LogSqrt2pi=-LogGamma0[1]

LogGamma=Function.new{"|s| 
  if s.real<1
    self[s+1]-Log[s]
  else
    LogGamma0[s]+LogSqrt2pi
  end
  "}
Gamma=Function.new{"|s| 
  if s.real<1
    self[s+1]/s
  elsif s.kind_of?(Integer)
    Fact[s-1]
  elsif s.kind_of?(Float) && s.to_i==s
    Fact[s.to_i-1]
  elsif s.Float?
    Exp[LogGamma0[s]]*sqrt(2*PI)
  else
    Exp[LogGamma[s]]
  end
  "}
#  InverseGamma[s]=1/Γ(s) 
InverseGamma=Function.new{"|s| 
  if s.real<1
    InverseGamma[s+1]*s
  elsif s.kind_of?(Integer)
    1/Fact[s-1]
  elsif s.kind_of?(Float) && s.to_i==s
    1/Fact[s.to_i-1]
  elsif s.Float?
    Exp[-LogGamma0[s]]/sqrt(2*PI)
  else
    Exp[-LogGamma[s]]
  end
  "}

def Math.Gamma(s,e=2**-53)
  Gamma[s].approx(e)
end

BesselT=Function.new{"|n| 
  Function.new({},n){\"|x,n| 
    Sequence.new([],x,n){'|i,x,n| (x/2)**(n+2*i)/Gamma[n+i+1]'}/Fact \"}
  "}
Bessel=Function.new{"|n| 
  Function.new({},n){\"|x,n| Real.new(BesselT[n][x],'A')\"}  "}

----