fortran66のブログ

fortran について書きます。

ζ(3) アペリーの定数

中央二項係数でアペリーの定数

中央二項係数の逆数和(その6)
より、
ζ(3)=Σ1/n^3=5/2Σ(-1)^(n-1)/n^3(2n,n)
を用いて、ζ(3)アペリーの定数を計算してみます。

実行結果

 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.
 zeta(3) Apery's constant
 s = sum 1/((-1)**(i + 1) * 2iCi * i^2)
   1.20205690315959428539973816151145
続行するには何かキーを押してください . . .


最後の桁まであっています。
ζ(3) = 1.202056903159594285399738161511449990764986292
Apéry's constant - Wikipedia



ソース・プログラム

    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_int64) == 0) xCy = [xCy, tri(i / 2 + 1)] ! 2mCm
      end do  
      !
      print *, 'Central Binomial Coefficients 2mCm'
      print '(i4, 1x, 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**3)   
      end do    
      s = s * 5 / 2.0_kq
      print *, 'zeta(3) Apery''s constant'
      print *, 's = sum 1/((-1)**(i + 1) * 2iCi * i^2)', s
    end program pascal