fortran66のブログ

fortran について書きます。

中央二項係数 再び

中央二項係数

中央二項係数とは、二項係数の丁度真ん中の値になります。2nCn とも書き表せます。
以前も計算したのですが、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  を求めます。
中央二項係数の逆数和(その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 

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Modern Fortran Explained (Numerical Mathematics and Scientific Computation)

Numerical Computing with Modern Fortran (Applied Mathematics)

Numerical Computing with Modern Fortran (Applied Mathematics)

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77

Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77