```正木です。

Γ 関数と 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')\"}  "}

----

```