serial-correlation-test.rb

require "mathn"

# 線形合同法の周期全体にわたるシリアル相関係数 (1つ次の値との相関係数)を求める
# TAOCP参照

# 「良い乗数」とされているものがとくにいい結果になるわけでもないようだ

def sigma(h, k, c)
  if k == 1
    0
  else
    - sigma(k % h, h, c % h) + f(h, k, c)
  end
end

def f(h, k, c)
  h / k + k / h + 1 / (h*k) + 6 * c * c / (h * k) - 6 * (c / h).floor - 3 * e(h, c)
end

def e(h, c)
  (c == 0 ? 1 : 0) + (c % h != 0 ? 1 : 0)
end

srand 0
a = [69069, 1664525, 39894229, 48828125, 1566083941, 1812433253, 2100005341, 1]
10.times { a << rand(2**30) * 4 + 1 }
a.each do |x|
  puts "%.8x: %s" % [x, sigma(x, 2**32, 1).to_f]
end

~result

00010dcd: 62174.535943921655
0019660d: 2349.9696628712118
0260bcd5: 100.54606653749943
02e90edd: 207.35500532761216
5d588b65: 47.27574636787176
6c078965: -2.4704124704003334
7d2b89dd: 41.992662850767374
00000001: 4294967290.0
31fc2ab1: 125.99428170360625
5f12a8bd: 468.62229411676526
dc5a99d5: -150.3242875188589
60873301: 0.3630951661616564
693acd0d: 5.522999849170446
6e894bed: -53.10637180879712
2df5db0d: 3.014538612216711
6395f59d: 31.868721161037683
b1d29025: 8.611387364566326
7e877b4d: 19.70915138348937