複素平面で有効なガンマ関数を数値計算します。
目指すは大浦様が公開してくださっている複素ガンマ関数cdgamma,cqgammaの計算速度を上回ることです。
結果として、本稿の複素ガンマ関数の計算コードは大浦様のコードの0.85~0.9倍ほどの時間で計算できるコードになりました。
大浦様のコードは、Lanczos近似という方法に基づいています。
Lanczos近似はガンマ関数を
\(
\begin{align}
\Gamma(z)=\frac{(z+g-0.5)^{z-0.5}}{e^{z+g-0.5}}L_g(z) \\
L_g(z)=C_0+\sum_{k=1}^{N-1}\frac{C_N}{z+k-1}
\end{align}
\)
の形で展開します[6]。係数\(C_k\)は定数\(g\)の値に依存します。
定数\(g\)は魔法の数のようなもので、値によって収束具合が異なります。このあたりの理論に関しては[5-6]が詳しいのでそちらをご参考ください。
おおよそ10-20項程度で収束し、非常にシンプルなコードになっているのが特徴です。
同じ考えでは超えられるはずがないので、ガンマ関数の特性をフルに使って計算していきます。
実数のみのガンマ関数が欲しい場合、gfortranでは組み込み関数dgamma, gammaを最適化オプション -Oまたは-O3 とともにコンパイルすると500分の1くらいの時間で計算できるため、圧倒的に早いです。
[adsense1]
方針
方針としては
- 漸化式
- 相反公式
- 逆ガンマ関数のTaylor展開
- 対数ガンマ関数の漸近展開
を用いて計算していきます。
詳細は以下に折りたたんでおきます。
まず、
漸化式
\(
\Gamma(z+1)=z\Gamma(z)
\)
を用いて、実部の範囲を\(0\lt Re(z)\lt 1\)の空間に移します。
具体的に、\(z=3.4\)とかでしたら、
\(
\Gamma(3.4)=2.4\times 1.4\times 0.4 \times \Gamma(0.4)
\)
として移すのです。もしもzが負でしたら、
相反公式
\(
\displaystyle \Gamma(1-z)\Gamma(z)=\frac{\pi}{\sin(\pi z)}
\)
を用いて、\(\Gamma(1-z)\)が\([0,1]\)になるようにします。
最後に
相反公式で\(\Gamma(z)\)を計算します。
\(0\lt Re(z)\lt 1\)の範囲では、逆ガンマ関数のTaylor展開を用います。
ここで、逆ガンマ関数\(f(z)\)とは単にガンマ関数の逆数をとったもので、
\(
\displaystyle f(z)=\frac{1}{\Gamma(z)}
\)
と定義され[2,9],
\(
\displaystyle f(z)=\sum_{k=1}^{\infty} a_k z^k
\)
と展開されます。収束半径は無限(\(|z|\lt \infty\))ですが、数値計算では\(|z|\lt 1\)の範囲辺りまでにするべきでしょう。
この逆ガンマ関数の利点は、z=0で発散せず、大きい値を取らない事です。
逆ガンマ関数のz=0回りのTaylor展開の係数は
\(
\begin{align}
a_1&=1 \\
a_2&=\gamma \\
a_k&=k a_1 a_k-a_2 a_{k-1}+\sum_{j=2}^k (-1)^j \zeta(j)a_{k-j}, k\gg 2
\end{align}
\)
で与えられます[1,3,9]。ここで、\(\gamma\)はオイラーの定数、\(\zeta(s)\)はリーマンのゼータ関数です。
この係数を求めるために、ゼータ関数の値をwolframから持ってきました。javaによる任意精度演算を行うと、係数は
a(1)=1q0
a(2)=0.5772156649015328606065120900824024310421593359399235988057672q0
a(3)=-0.655878071520253881077019515145390481279766380478584347292363q0
a(4)=-0.42002635034095235529003934875429818711394500401106093522067q-1
a(5)=0.166538611382291489501700795102105235717781502247174340570468q0
a(6)=-0.42197734555544336748208301289187391301652684189822486376919q-1
a(7)=-0.9621971527876973562114921672348198975362942252113002105139q-2
a(8)=0.7218943246663099542395010340446572709904800880238318001094q-2
a(9)=-0.1165167591859065112113971084018388666809333795384057443408q-2
a(10)=-0.215241674114950972815729963053647806478241923378338750350q-3
a(11)=0.128050282388116186153198626328164323394892099693677214901q-3
a(12)=-0.20134854780788238655689391421021818382294833297979115262q-4
a(13)=-0.1250493482142670657345359473833092242322655621153959815q-5
a(14)=0.1133027231981695882374129620330744943324004838621075654q-5
a(15)=-2.05633841697760710345015413002057283651257902629337946q-7
a(16)=6.116095104481415817862498682855342867275865719712321q-9
a(17)=5.002007644469222930055665048059991303044612742494482q-9
a(18)=-1.181274570487020144588126565436505577738759504932588q-9
a(19)=1.04342671169110051049154033231225019140070982312582q-10
a(20)=7.782263439905071254049937311360777226068086181393q-12
a(21)=-3.696805618642205708187815878085766236570963451362q-12
a(22)=5.10037028745447597901548132286323180272688606971q-13
a(23)=-2.0583260535665067832224295448552374197460910808q-14
a(24)=-5.348122539423017982370017318727939948989715478q-15
a(25)=1.226778628238260790158893846622422428165455750q-15
a(26)=-1.18125930169745876951376458684229783121155729q-16
a(27)=1.186692254751600332579777242928674071088494q-18
a(28)=1.412380655318031781555803947566709037086351q-18
a(29)=-2.29874568443537020659247858063369926028451q-19
a(30)=1.7144063219273374333839633702672570668127q-20
a(31)=1.33735173049369311486478139512226802287q-22
a(32)=-2.05423355176667278932502535135573379668q-22
a(33)=2.7360300486079998448315099043309820149q-23
a(34)=-1.732356445910516639057428451564779799q-24
a(35)=-2.3606190244992872873434507354275310q-26
a(36)=1.8649829417172944307184131618786669q-26
a(37)=-2.218095624207197204399716913626860q-27
a(38)=1.29778197494799366882441448633059q-28
a(39)=1.180697474966528406222745415510q-30
a(40)=-1.124584349277088090293654674262q-30
a(41)=1.27708517514086620399020667775q-31
a(42)=-7.391451169615140823461289330q-33
a(43)=1.1347502575542157609541653q-35
a(44)=4.6391346410587220299448049q-35
a(45)=-5.347336818439198875077418q-36
a(46)=3.20799592361335262286124q-37
a(47)=-4.445829736550756882102q-39
a(48)=-1.311174518881988712901q-39
a(49)=1.64703335254381388682q-40
a(50)=-1.0562331785035812186q-41
となります。qは指数を表します。
javaプログラムは↓
import java.math.BigDecimal;
public class ak {
public static void main(String[] args) {
int k,j;
BigDecimal ans = new BigDecimal("0");
BigDecimal s=new BigDecimal("0");
BigDecimal m1=new BigDecimal("-1");
BigDecimal dk=new BigDecimal("0");
BigDecimal a[]=new BigDecimal[51];
BigDecimal zeta[]=new BigDecimal[51];
ans = ans.add(new BigDecimal("0.1"));
ans = ans.add(new BigDecimal("0.11"));
ans = ans.add(new BigDecimal("0.222"));
ans = ans.add(new BigDecimal("0.3333"));
ans = ans.add(new BigDecimal("0.44444"));
ans = ans.add(new BigDecimal("0.555555"));
zeta[ 2] = new BigDecimal("1.64493406684822643647241516664602518921894990120679843773556");
zeta[ 3] = new BigDecimal("1.20205690315959428539973816151144999076498629234049888179227");
zeta[ 4] = new BigDecimal("1.08232323371113819151600369654116790277475095191872690768298");
zeta[ 5] = new BigDecimal("1.03692775514336992633136548645703416805708091950191281197419");
zeta[ 6] = new BigDecimal("1.01734306198444913971451792979092052790181749003285356184241");
zeta[ 7] = new BigDecimal("1.00834927738192282683979754984979675959986356056523870641728");
zeta[ 8] = new BigDecimal("1.00407735619794433937868523850865246525896079064985002032911");
zeta[ 9] = new BigDecimal("1.00200839282608221441785276923241206048560585139488875654860");
zeta[10] = new BigDecimal("1.00099457512781808533714595890031901700601953156447751725779");
zeta[11] = new BigDecimal("1.00049418860411946455870228252646993646860643575820861711914");
zeta[12] = new BigDecimal("1.00024608655330804829863799804773967096041608845800340453304");
zeta[13] = new BigDecimal("1.00012271334757848914675183652635739571427510589550984513670");
zeta[14] = new BigDecimal("1.00006124813505870482925854510513533374748169616915454948276");
zeta[15] = new BigDecimal("1.00003058823630702049355172851064506258762794870685817750657");
zeta[16] = new BigDecimal("1.00001528225940865187173257148763672202323738899047153115311");
zeta[17] = new BigDecimal("1.00000763719763789976227360029356302921308824909026267909538");
zeta[18] = new BigDecimal("1.00000381729326499983985646164462193973045469721895333114317");
zeta[19] = new BigDecimal("1.00000190821271655393892565695779510135325857114483863023593");
zeta[20] = new BigDecimal("1.00000095396203387279611315203868344934594379418741059575006");
zeta[21] = new BigDecimal("1.00000047693298678780646311671960437304596644669478493760021");
zeta[22] = new BigDecimal("1.00000023845050272773299000364818675299493504182177965826985");
zeta[23] = new BigDecimal("1.00000011921992596531107306778871888232638725499778451985860");
zeta[24] = new BigDecimal("1.00000005960818905125947961244020793580122750391883730279586");
zeta[25] = new BigDecimal("1.00000002980350351465228018606370506936601184473091954331240");
zeta[26] = new BigDecimal("1.00000001490155482836504123465850663069862886478816788591055");
zeta[27] = new BigDecimal("1.00000000745071178983542949198100417060411945471903188256583");
zeta[28] = new BigDecimal("1.00000000372533402478845705481920401840242323289305929581152");
zeta[29] = new BigDecimal("1.00000000186265972351304900640390994541694806166533046920067");
zeta[30] = new BigDecimal("1.00000000093132743241966818287176473502121981356795513681619");
zeta[31] = new BigDecimal("1.00000000046566290650337840729892332512200710626918533694731");
zeta[32] = new BigDecimal("1.00000000023283118336765054920014559759404950248298228453031");
zeta[33] = new BigDecimal("1.00000000011641550172700519775929738354563095165224717276359");
zeta[34] = new BigDecimal("1.00000000005820772087902700889243685989106305417312260461722");
zeta[35] = new BigDecimal("1.00000000002910385044497099686929425227884046410698198743303");
zeta[36] = new BigDecimal("1.00000000001455192189104198423592963224531842098380889412404");
zeta[37] = new BigDecimal("1.00000000000727595983505748101452086901233805926485092555547");
zeta[38] = new BigDecimal("1.00000000000363797954737865119023723635587327351264602838490");
zeta[39] = new BigDecimal("1.00000000000181898965030706594758483210073008503058930961866");
zeta[40] = new BigDecimal("1.00000000000090949478402638892825331183869490875386000099088");
zeta[41] = new BigDecimal("1.00000000000045474737830421540267991120294885703390452991144");
zeta[42] = new BigDecimal("1.00000000000022737368458246525152268215779786912138298219892");
zeta[43] = new BigDecimal("1.00000000000011368684076802278493491048380259064374359028425");
zeta[44] = new BigDecimal("1.00000000000005684341987627585609277182967524068553057158899");
zeta[45] = new BigDecimal("1.00000000000002842170976889301855455073704942662074368826531");
zeta[46] = new BigDecimal("1.00000000000001421085482803160676983430714173953767869860563");
zeta[47] = new BigDecimal("1.00000000000000710542739521085271287735447995680002274204359");
zeta[48] = new BigDecimal("1.00000000000000355271369133711367329846953405934299214565550");
zeta[49] = new BigDecimal("1.00000000000000177635684357912032747334901440027957015550858");
zeta[50] = new BigDecimal("1.00000000000000088817842109308159030960913863913863256088715");
a[0]=new BigDecimal("0");
a[1]=new BigDecimal("1");
a[2]=new BigDecimal("0.5772156649015328606065120900824024310421593359399235988057672348848677");
for(k=3;k<=50;k++){
s = new BigDecimal("0");
dk = BigDecimal.valueOf(k);
for(j=2;j<=k;j++){
s=s.add(((m1.pow(j)).multiply(zeta[j])).multiply(a[k-j]));
}
a[k]=(m1.multiply(a[2]).multiply(a[k-1])).add(s);
a[k]=a[k].divide(a[1].subtract(dk), 60, BigDecimal.ROUND_HALF_UP);
System.out.println(a[k].toString());
}
}
}
上記方法では\(z\)の実部の値はいくらでも\([0,1]\)の範囲に持っていくことが出来ますが、虚部が大きいと使うことはできません。
そこで、対数ガンマ関数の漸近展開を用います。
対数ガンマ関数\(\ln{\Gamma(z)}\)はガンマ関数の対数をとったものです。漸近領域では[8]より、
\(
\begin{align}
\ln\Gamma(z)\sim \left(z-\frac{1}{2}\right)\ln{z} -z+\frac{1}{2} \ln(2\pi)+\sum_{m=1}^{\infty}\frac{B_{2m}}{2m(2m-1)z^{2m-1}}
\end{align}
\)
と表されます。これにexpをとったものが答えとなります。
係数
\(
\displaystyle \frac{B_{2m}}{2m(2m-1)}
\)
の値は、wolframより求められますtable[(bernoulli b(2*i))/(2*i*(2*i-1)), {i,0,50}] with 36 digits。
\(z\)が程々の値、すなわち逆ガンマ関数のTaylor展開が使ないくらい遠く、漸近展開が使えないくらい近い場合を考えます。
この場合、漸化式
\(
\displaystyle \Gamma(z)=\frac{1}{z}\Gamma(z+1)
\)
を用いて\(|z|\)を大きな\(|z|\)へ移します。
その後、漸近展開を用います。
補足
数値計算で
\(
\displaystyle s=\sum_{k=1}^{N} a(k)z^{mk}
\)
を計算する際の効率の良い計算方法は
w=z**m
s=a(N)
do i=N-1,1,-1
s=s*w+a(i)
enddo
s=s*w
とすることです。doループの中で乗数が入らないため、高速です。
しかし、\(N\)を逐次的に計算し増やす事はできません。
———–折りたたみここまで———–
領域ごとに使っている関係式をまとめると、このような感じです。
※範囲変更(cgamma,円の大きさ8→9)
考え方は以上です。
実際にFortran90でプログラムしてみます。
fortran90によるプログラム
用意してあるのは、
複素倍精度用のガンマ関数cgamma
(本当のガンマ関数値との相対誤差 < 2e-14 [-10,10]の範囲, 8e-14 [-40,40]の範囲)
と
4倍精度複素ガンマ関数zgamma
(本当のガンマ関数値との相対誤差 < 3e-32 [-10,10]の範囲, 8e-32 [-40,40]の範囲)
です。
2017/01/10 逆ガンマ関数の項数の最適化
2017/10/27 正の整数の時、1ずれることの修正
function cgamma(z)
! sikinote (http://slpr.sakura.ne.jp/qp)
! author : sikino
! date : 2017/01/09
! 2017/01/10 optimize the number of a(M)
! 2017/10/27 fix for positive integer
implicit none
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::cgamma
integer::i,n,M
double precision,parameter::pi=3.14159265358979324d0
double precision,parameter::e1=0.3678794411714423216d0
double precision,parameter::ln2pi2=0.91893853320467274d0
double precision,parameter::sq2pi=2.5066282746310005d0
double precision::a(1:30),AB(1:14),dz,iz
complex(kind(0d0))::q,s,w,r,q1,q2
dz=dble(z)
iz=imag(z)
if(dz.eq.nint(dz).and.iz.eq.0d0)then
if(dz.gt.0d0)then
s=dcmplx(1d0,0d0)
n=nint(dz)
do i=2,n-1
s=s*i
enddo
else
s=1d+300
endif
cgamma=s
else
q=z
if(dz.lt.0d0)q=1d0-z
if(abs(iz).lt.1.45d0)then
! For x=arb, |y|\lt 1.45
n=int(dble(q))
s=1d0
do i=1,n
s=s*(q-i)
enddo
w=q-dble(n)
a(1) = 1d0
a(2) = 0.57721566490153286d0
a(3) =-0.65587807152025388d0
a(4) =-0.42002635034095236d-1
a(5) = 0.16653861138229149d0
a(6) =-0.42197734555544337d-1
a(7) =-0.96219715278769736d-2
a(8) = 0.72189432466630995d-2
a(9) =-0.11651675918590651d-2
a(10)=-0.21524167411495097d-3
a(11)= 0.12805028238811619d-3
a(12)=-0.20134854780788239d-4
a(13)=-0.12504934821426707d-5
a(14)= 0.11330272319816959d-5
a(15)=-0.20563384169776071d-6
a(16)= 0.61160951044814158d-8
a(17)= 0.50020076444692229d-8
a(18)=-0.11812745704870201d-8
a(19)= 0.10434267116911005d-9
a(20)= 0.77822634399050713d-11
a(21)=-0.36968056186422057d-11
a(22)= 0.51003702874544760d-12
a(23)=-0.20583260535665068d-13
a(24)=-0.53481225394230180d-14
a(25)= 0.12267786282382608d-14
a(26)=-0.11812593016974588d-15
a(27)= 0.11866922547516003d-17
a(28)= 0.14123806553180318d-17
a(29)=-0.22987456844353702d-18
a(30)= 0.17144063219273374d-19
M=int(11.3*abs(w)+13)
if(M.gt.30)M=30
r=a(M)
do i=M-1,1,-1
r=r*w+a(i)
enddo
cgamma=s/(r*w)
else
! For x=arb, |y|\gt 1.45
s=1d0
if(abs(q).lt.9d0)then
do i=0,7
s=s*(q+i)
enddo
s=1d0/s
q=q+8d0
endif
AB(1) = 0.83333333333333333d-1
AB(2) =-0.27777777777777778d-2
AB(3) = 0.79365079365079365d-3
AB(4) =-0.59523809523809524d-3
AB(5) = 0.84175084175084175d-3
AB(6) =-0.19175269175269175d-2
AB(7) = 0.64102564102564103d-2
AB(8) =-0.29550653594771242d-1
q1=1d0/q
q2=q1*q1
r=AB(8)
do i=7,1,-1
r=r*q2+AB(i)
enddo
cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
endif
if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
endif
return
end function cgamma
function zgamma(z)
! sikinote (http://slpr.sakura.ne.jp/qp)
! author : sikino
! date : 2017/01/09
! 2017/01/10
! 2017/10/27 fix for positive integer
implicit none
complex*32,intent(in)::z
complex*32::zgamma
integer::i,n,M
real*16,parameter:: pi=3.14159265358979323846264338327950288q0
real*16,parameter:: e1=0.36787944117144232159552377016146087q0
real*16,parameter::ln2pi2=0.91893853320467274178032973640561764q0
real*16,parameter:: sq2pi=2.50662827463100050241576528481104525q0
real*16::a(1:50),AB(1:30),dz,iz,aw
complex*32::q,s,w,r,q1,q2
dz=real(z); iz=imag(z)
if(dz.eq.nint(dz).and.iz.eq.0q0)then
if(dz.gt.0q0)then
s=dcmplx(1q0,0q0)
n=nint(dz)
do i=2,n-1
s=s*i
enddo
else
s=1q+3000
endif
zgamma=s
else
q=z
if(dz.lt.0q0)q=1q0-z
if(abs(imag(q)).lt.1.35q0)then
! For x=arb, |y|\lt 1.35
n=int(real(q))
s=1q0
do i=1,n
s=s*(q-i)
enddo
w=q-real(n)
a(1) = 1q0
a(2) = 0.577215664901532860606512090082402431q0
a(3) =-0.655878071520253881077019515145390481q0
a(4) =-0.420026350340952355290039348754298187q-1
a(5) = 0.166538611382291489501700795102105236q0
a(6) =-0.421977345555443367482083012891873913q-1
a(7) =-0.962197152787697356211492167234819898q-2
a(8) = 0.721894324666309954239501034044657271q-2
a(9) =-0.116516759185906511211397108401838867q-2
a(10)=-0.215241674114950972815729963053647806q-3
a(11)= 0.128050282388116186153198626328164323q-3
a(12)=-0.201348547807882386556893914210218184q-4
a(13)=-0.125049348214267065734535947383309224q-5
a(14)= 0.113302723198169588237412962033074494q-5
a(15)=-2.056338416977607103450154130020572837q-7
a(16)= 6.116095104481415817862498682855342867q-9
a(17)= 5.002007644469222930055665048059991303q-9
a(18)=-1.181274570487020144588126565436505578q-9
a(19)= 1.043426711691100510491540332312250191q-10
a(20)= 7.782263439905071254049937311360777226q-12
a(21)=-3.696805618642205708187815878085766237q-12
a(22)= 5.100370287454475979015481322863231803q-13
a(23)=-2.058326053566506783222429544855237420q-14
a(24)=-5.348122539423017982370017318727939949q-15
a(25)= 1.226778628238260790158893846622422428q-15
a(26)=-1.181259301697458769513764586842297831q-16
a(27)= 1.186692254751600332579777242928674071q-18
a(28)= 1.412380655318031781555803947566709037q-18
a(29)=-2.298745684435370206592478580633699260q-19
a(30)= 1.714406321927337433383963370267257067q-20
a(31)= 1.337351730493693114864781395122268023q-22
a(32)=-2.054233551766672789325025351355733797q-22
a(33)= 2.736030048607999844831509904330982015q-23
a(34)=-1.732356445910516639057428451564779799q-24
a(35)=-2.360619024499287287343450735427531031q-26
a(36)= 1.864982941717294430718413161878666894q-26
a(37)=-2.218095624207197204399716913626860384q-27
a(38)= 1.297781974947993668824414486330595180q-28
a(39)= 1.180697474966528406222745415509988346q-30
a(40)=-1.124584349277088090293654674261653627q-30
a(41)= 1.277085175140866203990206677750563204q-31
a(42)=-7.391451169615140823461289330001078362q-33
a(43)= 1.134750257554215760954165252401478595q-35
a(44)= 4.639134641058722029944804894433875948q-35
a(45)=-5.347336818439198875077418068570325956q-36
a(46)= 3.207995923613352622861238213439296964q-37
a(47)=-4.445829736550756882101576639628564322q-39
a(48)=-1.311174518881988712901096702293607192q-39
a(49)= 1.647033352543813886818056435068597372q-40
a(50)=-1.056233178503581218582902184722885523q-41
aw=abs(w)
if(aw.le.0.42q0)then
M=int(35q0*aw+20q0)+1
else
M=int(15.5q0*aw+28.2q0)+1
endif
if(M.gt.50)M=50
r=a(M)
do i=M-1,1,-1
r=r*w+a(i)
enddo
zgamma=s/(r*w)
else
! For x=arb, |y|\ge 1.35
s=1q0
if(abs(q).lt.18q0)then
do i=0,16
s=s*(q+i)
enddo
s=1q0/s
q=q+17q0
endif
AB(1) = 0.83333333333333333333333333333333333q-1
AB(2) =-0.27777777777777777777777777777777778q-2
AB(3) = 0.79365079365079365079365079365079365q-3
AB(4) =-0.59523809523809523809523809523809524q-3
AB(5) = 0.84175084175084175084175084175084175q-3
AB(6) =-0.19175269175269175269175269175269175q-2
AB(7) = 0.64102564102564102564102564102564103q-2
AB(8) =-0.29550653594771241830065359477124183q-1
AB(9) = 0.17964437236883057316493849001588940q0
AB(10)=-1.3924322169059011164274322169059011q0
AB(11)= 13.402864044168391994478951000690131q0
AB(12)=-156.84828462600201730636513245208897q0
AB(13)= 2193.1033333333333333333333333333333q0
AB(14)=-36108.771253724989357173265219242231q0
AB(15)= 691472.268851313067108395250775673468q0
AB(16)=-1.52382215394074161922833649588867805q7
q1=1q0/q; q2=q1*q1
r=AB(16)
do i=15,1,-1
r=r*q2+AB(i)
enddo
zgamma=s*exp((q-0.5q0)*log(q)-q+ln2pi2+r*q1)
endif
if(dz.lt.0q0)zgamma=pi/(zgamma*sin(pi*z))
end if
return
end function zgamma
[adsense2]
計算速度の比較
大浦様のコードと計算速度を比較します。
比較方法は、
乱数を[-20,20]の範囲で40000点振り、それを10回繰り返した時の平均cpu時間を計測します。
上の計算を20回繰り返し、その平均を取ることによって計算時間を比較します。
その結果は以下のようになりました。
大浦様のコードによる結果は黒、私のコードによる結果は赤です。
大浦様のコードと比べて、私のプログラムの方が計算時間が0.85~0.9倍になっていることが分かります。
嬉しい結果となりました。満足です。
コードの利用に関しては連絡先とサイトについてをご覧下さい。
ガンマ関数の概形
最後に、ガンマ関数の概形です。
参考文献
[1] Gamma Function -mathworld
[2]Reciprocal gamma function -wikipedia
[3]Gamma function – OeisWiki
[4] Wrench, J.W. (1968). Concerning two series for the gamma function. Mathematics of Computation, 22, 617–626
[5]Wrench, J.W. (1973). Erratum: Concerning two series for the gamma function. Mathematics of Computation, 27, 681–682
※[4]の係数表には間違いがあり、同著者によって[4]で訂正されています。
[6]The Lanczos Approximation -Boost
[7]Lanczos_approximation -wikipedia
[8] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 6.1.40, Page 257)
[8] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 6.1.34, Page 256)
