ガンマ関数の数値計算

複素平面で有効なガンマ関数を数値計算します。
目指すは大浦様が公開してくださっている複素ガンマ関数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展開
  • 対数ガンマ関数の漸近展開

を用いて計算していきます。

詳細は以下に折りたたんでおきます。

領域ごとに使っている関係式をまとめると、このような感じです。

※範囲変更(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)

Javaによる任意精度演算

javaによる任意精度演算のプログラムです。
fortranでやろうとしましたがサポートされてないようです。

私はjavaに精通していません。とりあえず何のプログラムでもいいので任意精度演算を使おうとした際、javaが一番手軽だっただけです。

ここでは、ベルヌーイ数を計算することを考えます。
擬コードとして、wikipediaより、

for m in 0, 1, …, n do
     A[m] ← 1 / (m + 1)
     for j in m, m − 1, …, 1 do
         A[j − 1] ← j × (A[j − 1] - A[j])
     end for
end for
return A[0]    # which is Bn

でベルヌーイ数が得られます。
倍精度演算(10進数で16桁)では桁落ちのため\(B_{14}\)で6桁しか一致しません。

計算過程で100桁まで用いて(割り算で100桁目で四捨五入)、\(B_{50}\)を計算してみましょう。

import java.math.BigDecimal;
 
public class bernoulli {

  final static int divideN=100;
 
  public static void main(String[] args) {
    int m,j;
    BigDecimal b  = new BigDecimal("0");
    BigDecimal bm = new BigDecimal("0");
    BigDecimal a[]= new BigDecimal[51];

    a[0]=new BigDecimal("0");

    for(m=0;m<=50;m++){
       bm = BigDecimal.valueOf(m);
         a[m]=(new BigDecimal("1")).divide(bm.add(new BigDecimal("1")), divideN, BigDecimal.ROUND_HALF_UP);
       for(j=m;j>=1;j--){
             a[j-1]=(BigDecimal.valueOf(j)).multiply(a[j-1].subtract(a[j]));
         }
       b=a[0];
       String text=String.format("%.36e",b);
       System.out.println(text);
      }
    }
}

コンパイルと実行は

javac bernoulli.java
java bernoulli

です。

出力は36桁としています。上記プログラムによる結果は、
7.500866746076964366855720075757575758e+24
となり、wolframによる値Table[bernoulli b(i), {i,0,50}] with 36 digits
7.50086674607696436685572007575757576e+24
です。全桁一致しました。
[adsense1]


以下、javaの任意精度演算に関するメモです。

メモ

変数の定義方法
  • 整数型
    int i,j;
  • 倍精度実数型
    double x;
  • 任意精度型
    BigDecimal a,b,c;
配列の定義方法
BigDecimal r[]= new BigDecimal[51];

↑で配列r[0]~r[50]が確保される。

配列の初期化

for文でまわしましょう。

for(i=0;i<=N;i++){
    r[i]=new BigDecimal("0");
}
変数へ値の代入

aに任意精度型で\(1.2000000000\ldots\)を代入する。

a=new BigDecimal("1.2");
足し算
\(
a=b+1.2
\)
a=b.add(new BigDecimal("1.2"))

(整数型の型変換)
\(
a=b+i
\)
\(i\)は整数型

a=b.add(BigDecimal.valueOf(i));

引き算
\(
a=b-1.2
\)
a=b.subtract(new BigDecimal("1.2"));

掛け算
\(
a=b\times c
\)
a=b.multiply(c);

割り算

\(
a=b/c
\)
丸めの桁を100桁にする時

a=b.divide(c, 100, BigDecimal.ROUND_HALF_UP);

足し算と割り算

\(
a=b/(c+9.92)
\)
丸めの桁を100桁にする時

a=b.divide(c.add(new BigDecimal("9.92")), 100, BigDecimal.ROUND_HALF_UP);

累乗
\(
a=b^c
\)
a=b.pow(c);

ループ

\(
\displaystyle s=\sum_{i=0}^{N-1}a_i
\)

b=new BigDecimal("0");
for(i=0;i<=N-1;i++){
    b=b.add(r[i]);
}

端末への指数表記,有効桁数36桁で出力
String text=String.format("%.36e",a);
System.out.println(text);

[adsense2]

参考文献


【Java】任意精度整数(多倍長整数)で演算する -ヽ|∵|ゝ(Fantom) の 開発blog?
Java BigDecimalの四則演算とフォーマット処理 -ITのおもちゃ箱
intとBigDecimalの型変換 -Javaプログラミング講座

逆行列の数値計算

Lapackを用いて数値計算で逆行列を求めます。
基本的に数値的にあらわに逆行列を求める意味はそんなにありません。
連立一次方程式を解きたければ、行列のLU分解を通じて解けば良いです。この場合、逆行列を求めることはしません。

そうはいっても使ってみたい時があります。そんな状況を考えましょう。

Lapack, もしくはインテルのマス・カーネル・ライブラリ(MKL)を用います。

用いるルーチンは

?getrf, ?getri

の2つです。

ここに用意したのは配列の大きさNと配列Aを入れると逆行列を出力するプログラムです。
倍精度、複素倍精度のルーチンを書きました。

ifort -mkl main.f90

でコンパイルしてください。

program main
  implicit none
  integer::i,j,N
  complex(kind(0d0)),allocatable::A(:,:),B(:,:)
 
  N=3
  allocate(A(1:N,1:N),B(1:N,1:N))
  A(1,1:3)=(/(2d0,0d0),(3d0,0d0),(0d0,0d0)/)
  A(2,1:3)=(/(2d0,0d0),(9d0,0d0),(1d0,3d0)/)
  A(3,1:3)=(/(3d0,4d0),(5d0,0d0),(3d0,1d0)/)

  B=A
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  call zinvmat(size(A,1),A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  A=matmul(B,A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  stop
end program main
   
subroutine invmat(N,A)
  implicit none
  integer,intent(in)::N
  double precision,intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  double precision::tw(1:1)
  double precision,allocatable::work(:)
 
  info=0
  call dgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call dgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(tw(1))
  allocate(work(1:lwork+1))
  work=0d0
 
  call dgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine invmat


subroutine zinvmat(N,A)
  implicit none
  integer,intent(in)::N
  complex(kind(0d0)),intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  complex(kind(0d0))::tw(1:1)
  complex(kind(0d0)),allocatable::work(:)
 
  info=0
  call zgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call zgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(dble(tw(1)))
  allocate(work(1:lwork+1))
  work=dcmplx(0d0,0d0)
 
  call zgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine zinvmat