中央二項係数　再び - fortran66のブログ

中央二項係数

中央二項係数とは、二項係数の丁度真ん中の値になります。2nＣn とも書き表せます。
以前も計算したのですが、32bit 整数の範囲だったのであまり多くは求められませんでした。
fortran66.hatenablog.com

今回 64bit 整数で計算してもいいのですが、ひとひねりして 128bit 実数の仮数部のビット数が 64bit より長いことを利用して、指数部を使わず仮数部の範囲で求めてみることにします。

また、せっかく求めたので中央二項係数を使って、特殊な数値を計算して検算とします。

ここのページを参考にして、 $\sum_{i=1}^\infty{(-1)^{(i-1)} \over {}_{2i}C_{i}\,i^2}= 2\{ \mathrm{arcsinh}({1\over2})\}^2$ を求めます。
中央二項係数の逆数和（その２）

実行結果

 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.
  16                               601080390.
  17                              2333606220.
  18                              9075135300.
  19                             35345263800.
  20                            137846528820.
  21                            538257874440.
  22                           2104098963720.
  23                           8233430727600.
  24                          32247603683100.
  25                         126410606437752.
  26                         495918532948104.
  27                        1946939425648112.
  28                        7648690600760440.
  29                       30067266499541040.
  30                      118264581564861424.
  31                      465428353255261088.
  32                     1832624140942590534.
  33                     7219428434016265740.
  34                    28453041475240576740.
  35                   112186277816662845432.
  36                   442512540276836779204.
  37                  1746130564335626209832.
  38                  6892620648693261354600.
  39                 27217014869199032015600.
  40                107507208733336176461620.
  41                424784580848791721628840.
  42               1678910486211891090247320.
  43               6637553085023755473070800.
  44              26248505381684851188961800.
  45             103827421287553411369671120.
  46             410795449442059149332177040.
  47            1625701140345170250548615520.
  48            6435067013866298908421603100.
  49           25477612258980856902730428600.
  50          100891344545564193334812497256.
  51          399608854866744452032002440112.
  52         1583065848125949175357548128136.
  53         6272525058612251449529907677520.
  54        24857784491537440929618523018320.
  55        98527218530093856775578873054432.
  56       390590044887157789360330532465784.
  57      1548655265692941410446222812934512.
  58      6141219157058215937976400809912720.
  59     24356699707654619143838606602026720.
  60     96614908840363322603893139521372656.
 s = sum 1/((-1)**(i + 1) * 2iCi * i^2)
  0.463129641154388784993858142463065
   = 2 * {arcsinh(1/2)}^2
  0.463129641154388784993858142463066
続行するには何かキーを押してください . . .

最後の桁以外は一致しています。

ソース・プログラム

    program pascal 
      use, intrinsic :: iso_fortran_env
      implicit none
      integer, parameter :: kq = real128
      real(kq), parameter :: pi = 4.0_kq * atan(1.0_kq)
      real(kq), allocatable :: tri(:), xCy(:)
      integer(int64) :: i
      real(kq) :: s 
      xCy = [real(kq)::]
      tri = [1]
      do i = 1, 120
        tri = [tri, 0.0_kq] + [0.0_kq, tri]
        if (mod(i, 2) == 0.0) xCy = [xCy, tri(i / 2 + 1)] ! 2mCm
      end do  
      !
      print *, 'Central Binomial Coefficients 2mCm'
      print '(i4, x, f40.0)', (i, xCy(i), i = 1, size(xCy))
 
      s = 0.0_kq
      do i = 1, size(xCy)
        s = s + (-1)**(i + 1) * 1.0_kq / (xCy(i) * i * i)   
      end do    
      print *, 's = sum 1/((-1)**(i + 1) * 2iCi * i^2)', s
      print *, '  = 2 * {arcsinh(1/2)}^2              ', 2 * asinh(0.5_kq)**2
    end program pascal