Intel Fortran 組み込み関数の検定

大体において、組み込み関数の乱数は腐っているもので、真面目な計算に使うと怖い人達に叱られてしまいます。Intel Visual Fortran コンパイラー XE 12.1.4 で軽く検証してみます。
Fortran90 で導入されたサブルーチン RANDOM_NUMBER(x) は、elemental なサブルーチンなので引数に配列を渡すと、全要素に[0..1)の乱数を計算して返してくれます。サブルーチン RANDOM_NUMBER の乱数は、平均 0.5 の完全な一様分散になっているようです。この場合標準偏差 1/12 を持つ*1ことが知られています。

計算してみたところ、単精度の時にはデータ数が 10^6 を超えたあたりから、分散が理論値からずれ始めてきます。特に debug mode の時に影響が強く出ているようです。データを倍精度にすると理論値通りになります。

乱数や統計の知識に乏しいため原因は不明ですが、HELP FILE を見る限り乱数の周期は 10^6 よりもはるかに長いようです。生成種を変えても結果が大きく変わりません。

追記[H24(2012)-05-28]
問題は乱数ではなく、組み込み関数 SUM の方にあるようです。これを Kahan-sum に置き換えてやると、ちゃんと理論値通り 12 の周りで揺らぐようになりました。

実行結果

32bitの結果。64bit版でも結果に大差なし。

ソース・プログラム

program test
  implicit none
  integer, parameter :: kd = kind(1.0d0)
  integer, parameter :: nmax = 10**7
  integer :: i
  real(kd) :: x(nmax), av, sd
  call random_seed()
  do i = 1, 100
    call random_number(x)
    av = sum(x) / nmax
    sd = sum((x - 0.5_kd)**2) / nmax 
    print *, i, av, 1.0_kd / sd
  end do
  stop
end program test

*1:[tex:\int_0^1(x-1/2)^2 dx]