fortran66のブログ

fortran について書きます。

中央二項係数

Pascal の三角形によって二項係数が求まります。その中で _{2m}C_m*1 にあたる要素だけを取り出します。
これを中央二項係数と呼びます。Pascal の三角形の中央に並んだ列に当たります。

次にこの中央二項係数を用いて、級数 \sum_{k=1}^\infty {1\over _{2k}C_k k^2}={\pi^2\over18}http://oeis.org/A086463 を計算してみます。

似たような級数に、有名な \sum_{k=1}^\infty {1\over k^2}={\pi^2\over6} がありますが、こちらは収束が遅いことで知られていて、加速法の例題に取り上げられたりします。これに比べると今回用いる式は収束がとても早くなっています。

実行結果

32ビット整数の範囲で中央二項係数を求めました。そのため有効桁数も10桁程度になっています。15項までの和でここまで収束しています。

 Central Binomial Coefficients 2mCm
  1         2
  2         6
  3        20
  4        70
  5       252
  6       924
  7      3432
  8     12870
  9     48620
 10    184756
 11    705432
 12   2704156
 13  10400600
 14  40116600
 15 155117520

 s = sum 1/(2iCi * i^2) 0.548311355607652
 s / pi^2, 1/18         5.555555555470207E-002 5.555555555555555E-002
Press any key to continue . . .

ソース・プログラム

    program pascal 
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, allocatable :: tri(:), m2n(:)
      integer :: i
      real(real64) :: s 
      m2n = [integer::]
      tri = [1]
      do i = 1, 30
        tri = [tri, 0] + [0, tri]
        if (mod(i, 2) == 0) m2n = [m2n, tri(i / 2 + 1)] ! 2mCm
      !   print '(*(i10))', tri
      end do  
      !
      print *, 'Central Binomial Coefficients 2mCm'
      print '(i3, i10)', (i, m2n(i), i = 1, size(m2n))
      print *
      
      s = 0.0
      do i = 1, size(m2n)
        s = s + 1.0_real64 / (real(m2n(i), real64) * i * i)   
      end do    
      print *, 's = sum 1/(2iCi * i^2)', s
      print *, 's / pi^2, 1/18        ', s / (4 * atan(1.0_real64))**2, 1.0_real64 / 18.0_real64
      stop

    end program pascal 

式の証明

メモ帳代わりに式の証明の概略をダラダラと記しておきます。

\sum{1\over _{2k}C_k k^2}=\sum{1\over n^2}{n! n! \over (2n!)}=\sum{1\over n}{\Gamma(n)\Gamma(n+1)\over\Gamma(n+1)}=\sum{1\over n}B(n,n+1)=\sum\bigint_0^1{1\over n}u^n(1-u)^{n-1}du
結局-\bigint_0^1{\log(u^2-u+1)\over u}du と変形できます。
ここで u=1-u の変数変換をへて、-\log(1-x)=x+{x^2\over2}+{x^3\over3}+\... を用いました。

次に u^2-u+1=(u-e^{i{\pi\over3}})(u-e^{-{i\pi\over3}})}因数分解できるので、積分範囲に注意しつつ対数の二次式を一次式に分離してやれば、求める積分は二重対数関数\mathrm{Li}_2(x)=-\int_0^x{\log(1-u)\over u}du を用いて、\mathrm{Li}_2(e^{i{\pi\over3}})+\mathrm{Li}_2(e^{-i{\pi\over3}}) となります。ここで二重対数関数の別表記 \mathrm{Li}_2(x)=\sum{x^k\over k^2} を使って、e^{\pm i{\pi\over3}複素平面での対称性を考えると、実数の級数を考えればよいことがわかります。

結局 {1\over 1^2}+{-1\over 2^2}+{-2\over 3^2}+{-1\over 4^2}+{1\over 5^2}+{2\over 6^2}+\... を求めることになります。これは -\sum{(-1)^k\over k^2}+\sum{3(-1)^k\over (3k)^2} と変形でき、\sum{(-1)^k\over k^2}=(1-{1\over2})\sum{1\over k^2}={\pi^2\over12} から、求める値 {\pi^2\over 18} が得られます。

眠いので間違ってたらごめんwww

*1:括弧による表記がうまく出来ないので _nC_k 型の表記を用います。