LOG(2) のバイナリー表示 - fortran66のブログ

九月二日の続き。
Pi Directory by David H. Bailey下のDavid H. Bailey, "The BBP Algorithm for Pi," manuscript, Sep 2006. に基づいて、LOG(2) の任意位置での二進表示の計算というものをやってみました。
アルゴリズムを理解できたかの試し計算なので、あんまり考えてません。

式は、以下の通りです。
$\{2^d\log2\}=\{\{\sum_{k=1}^n{2^{n-k} \over k}\}+\sum^\infty_{k=n+1}{2^{n-k} \over k}\}$
$=\{\{\sum_{k=1}^n{2^{n-k} \mathrm{mod} k \over k}\}+\sum^\infty_{k=n+1}{2^{n-k} \over k}\}$
ここで {} は、小数部を意味します。

ソースコード

PROGRAM test ! LOG(2) 
  IMPLICIT NONE
  INTEGER :: i, k, n
  LOGICAL :: qq(31)
  REAL(8) :: x, y, v, w
  CHARACTER(64) :: tmp
  DO n = 1, 200
! 1st term
   x = 0.0d0
   DO k = 1, n
    DO i = 1, 31
     qq(i) = (IBITS(n - k, i - 1, 1) == 1) 
    END DO
    v = 1.0d0
    w = 2.0d0
    DO i = 1, 31
     IF (qq(i)) v = v * w 
     w = MOD(w * w, REAL(k, 8))
    END DO
    x = x + v / REAL(k, 8)
   END DO
! 2nd term
   y = 0.0d0
   DO k = n + 1, n + 100
    y = y + 2.0d0**REAL(n - k, 8) / REAL(k, 8)
   END DO
!
   WRITE(tmp, '(b64.64)') MOD(x + y, 1.0d0) !, x, y
   PRINT '(i5, x, 3a)', n, tmp(1:12), ':', tmp(13:64)
  END DO

  STOP
END PROGRAM test

出力結果

出力は LOG(2) の二進表示の n+1 桁目から始まる数値の倍精度浮動小数のビット列表示です。指数部と仮数部の区切りに「:」を書いておきました。確かに、１ビットづつずれた結果が得られているのがわかります。（指数部が一定になるように調整しなかったので、ところどころガクガクになっていますが・・）