fortran66のブログ

fortran について書きます。

LOG(2) のバイナリー表示 昨日の続き

二進法表示で小数点以下 n+1 桁目から始まる LOG(2) のビット列表示です。

計算する式
\{2^d\log2\}=\{\{\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
  REAL(8) :: x, y, v, w
  DO n = 0, 1000
! 1st term
   x = 0.0d0
   DO k = 1, n
    v = 1.0d0
    w = 2.0d0
    DO i = 1, 31
     IF ( IBITS(n - k, i - 1, 1) == 1 ) v = v * w
     w = MOD(w * w, REAL(k, 8))
    END DO
    x = x + v / k
   END DO
! 2nd term
   y = 0.0d0
   v = 1.0d0
   DO k = 1, 100
    v = v * 2.0d0
    y = y + 1.0d0 / ( v * (k + n) )
   END DO
!
   PRINT '(i5, a, b30.30 )', n, ' : ', INT( 2**30 * MOD(x + y, 1.0d0) )
  END DO

  STOP
END PROGRAM test

出力結果

一ビットづつずれていっているのが分かります。
このプログラムでは倍精度で計算すると、n=700 あたりで有効数字が失われる感じです。これがもっともな数字のか吟味してません・・・