中央二項係数
中央二項係数とは、二項係数の丁度真ん中の値になります。2nCn とも書き表せます。
以前も計算したのですが、32bit 整数の範囲だったのであまり多くは求められませんでした。
fortran66.hatenablog.com
今回 64bit 整数で計算してもいいのですが、ひとひねりして 128bit 実数の仮数部のビット数が 64bit より長いことを利用して、指数部を使わず仮数部の範囲で求めてみることにします。
また、せっかく求めたので中央二項係数を使って、特殊な数値を計算して検算とします。
ここのページを参考にして、 を求めます。
中央二項係数の逆数和(その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)
- 作者: Michael Metcalf,John Reid,Malcolm Cohen
- 出版社/メーカー: Oxford Univ Pr (Txt)
- 発売日: 2011/05/19
- メディア: ペーパーバック
- この商品を含むブログを見る
Numerical Computing with Modern Fortran (Applied Mathematics)
- 作者: Richard J. Hanson,Tim Hopkins
- 出版社/メーカー: Society for Industrial and Applied Mathematics
- 発売日: 2014/01/16
- メディア: ペーパーバック
- この商品を含むブログを見る
Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77
- 作者: Ian Chivers,Jane Sleightholme
- 出版社/メーカー: Springer
- 発売日: 2015/08/17
- メディア: ペーパーバック
- この商品を含むブログを見る