XorShift256に使えるパラメタを列挙してみた

やったこと

• XorShift法で64ビット変数を4本使う場合のパラメタの列挙をするプログラム、作って動かした。
• 1ヶ月位かけて今朝未明にようやく終わったので記事化

ターゲット

• とりあえず高速で長い周期を持つ擬似乱数生成器が欲しい人
• XorShift法に関する研究がしたいが、パラメタの列挙から始めるのが面倒な人

テンプレート

std::uint64_t xorshift256(void) {
  static std::uint64_t x = 1, y = 0, z = 0, w = 0;  // 全ゼロ以外の値。種。
  std::uint64_t t = x ^ (x << A);
  x = y; y = z; z = w;
  return w = (w ^ (w >> C)) ^ (t ^ (t >> B));
}

パラメタ一覧

A	B	C
1	2	27
1	3	10
1	5	33
1	21	25
1	63	7
2	19	14
2	39	3
3	7	42
3	41	29
4	5	56
4	5	61
4	9	27
4	23	9
5	1	21
5	6	53
5	15	16
5	20	9
5	21	9
6	31	25
7	3	13
7	3	32
7	13	60
7	15	1
7	27	24
7	27	44
8	13	28
8	27	24
9	1	33
10	1	32
10	7	9
11	1	59
11	3	4
11	9	40
11	10	19
11	18	59
11	25	7
11	26	41
11	30	9
11	63	21
13	5	36
13	6	19
13	8	55
13	14	11
13	18	41
13	20	1
13	21	44
14	9	52
14	17	26
14	19	53
15	1	61
15	16	25
15	17	5
15	17	59
15	26	19
15	37	10
15	41	17
17	6	29
17	6	49
17	11	18
17	18	9
17	22	33
17	23	1
18	17	63
19	14	23
19	16	53
19	16	61
19	21	13
20	9	15
20	29	43
21	13	1
21	39	2
21	42	5
22	15	7
22	19	28
22	19	47
23	1	41
23	7	55
23	13	36
23	17	58
23	41	31
24	13	50
24	29	41
25	1	36
25	7	52
25	11	33
25	12	41
25	16	1
25	17	45
25	19	42
25	22	11
25	27	6
25	42	7
25	43	18
26	9	35
26	9	48
26	13	36
26	37	6
27	1	32
27	5	62
27	6	13
27	8	59
28	31	59
28	47	4
29	3	59
30	1	2
30	29	14
31	1	55
31	11	5
31	12	45
31	33	11
31	33	19
33	32	23
34	9	4
35	38	21
36	1	3
36	11	39
36	25	8
36	41	11
37	5	10
37	27	4
37	49	13
39	28	1
40	15	39
41	26	3
41	45	11
42	27	22
43	10	47
43	21	50
43	21	51
44	9	29
44	29	2
45	7	1
45	15	1
45	19	59
47	5	2
47	6	3
50	3	1
50	3	18
50	7	19
50	15	1
51	5	49
51	13	19
51	13	21
51	46	1
51	46	13
52	7	1
55	17	9
55	23	7
59	3	8
59	60	1
61	42	3
63	31	1

以上、152通り。

計算プログラム

#!/usr/bin/env ruby

require 'benchmark'
require 'matrix'
require 'parallel'

def unit_length
  64
end

def create_t(a, b, c)
  i = Matrix.identity(unit_length)
  o = Matrix.zero(unit_length)
  l = Matrix.build(unit_length) do |row, col|
    if row + 1 == col
      1
    else
      0
    end
  end
  r = Matrix.build(unit_length) do |row, col|
    if row == col + 1
      1
    else
      0
    end
  end
  g = Matrix[
    [o, o, o, (i + l**a) * (i + r**b)],
    [i, o, o, o],
    [o, i, o, o],
    [o, o, i, (i + r**c)]
  ]
  Matrix.build(unit_length * 4) do |row, col|
    g[row / unit_length, col / unit_length][row % unit_length, col % unit_length]
  end
end

@ds = [
  0x000000000000000000d3eafc3af14600ffffffffffffffffff2c1503c50eb9ff,
  0x000000000000013540775b48cc32ba00fffffffffffffecabf88a4b733cd45ff,
  0x0000000000042f00fffffffffffbd0ff0000000000042f00fffffffffffbd0ff,
  0x00000280fffffd7f00000280fffffd7f00000280fffffd7f00000280fffffd7f,
  0x00003d30f19cd100ffffc2cf0e632eff00003d30f19cd100ffffc2cf0e632eff,
  0x0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff,
  0x00663d80ff99c27f00663d80ff99c27f00663d80ff99c27f00663d80ff99c27f,
  0x00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff,
  0x0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f,
  0x3333333333333333333333333333333333333333333333333333333333333333,
  0x5555555555555555555555555555555555555555555555555555555555555555,
]

def pop32(n)
  n = (n & 0x55555555) + ((n >> 1) & 0x55555555)
  n = (n & 0x33333333) + ((n >> 2) & 0x33333333)
  n = (n & 0x0f0f0f0f) + ((n >> 4) & 0x0f0f0f0f)
  n = (n & 0x00ff00ff) + ((n >> 8) & 0x00ff00ff)
  (n & 0x0000ffff) + (n >> 16)
end

def pop(n)
  s = 0
  until n == 0
    s += pop32(n & 0xffffffff)
    n >>= 32
  end
  s
end

def multiply_bin(a,b)
  as = a.to_a.map do |item|
    item.join.to_i 2
  end
  bs = b.t.to_a.map do |item|
    item.join.to_i 2
  end
  Matrix.build(a.row_size, b.column_size) do |l, n|
    pop(as[l] & bs[n]) & 1
  end
end

def regular_bin_p(a)
  t = a.to_a
  len = a.row_size
  for y in 0...len
    if t[y][y] == 0
      processed = false
      for z in y+1...len
        unless t[z][y] == 0
          s = t[z]
          t[z] = t[y]
          t[y] = s
          processed = true
          break
        end
      end
      return false unless processed
    end

    for z in 0...len
      next if z == y
      unless t[z][y] == 0
        for x in 0...len
          t[z][x] ^= t[y][x]
        end
      end
    end
  end
  return true
end

def clear_cache
  @power_cache = {}
end

def highpower_bin(a,n)
  t = a
  until n == 0
    if n == 1
      return t
    else
      t = multiply_bin(t,t)
      n >>= 1
    end
  end
end

def power_bin(a,n)
  return highpower_bin(a,n) if (n & (n - 1)) == 0

  s = Matrix.identity(a.row_size)
  i = 1
  j = 0
  @power_cache[i] = a unless @power_cache[i]
  t = @power_cache[i]
  until n == 0
    if n & 1 == 1
      j += i
      @power_cache[j] = multiply_bin(s,t) unless @power_cache[j]
      s = @power_cache[j]
    end
    i *= 2
    @power_cache[i] = multiply_bin(t,t) unless @power_cache[i]
    t = @power_cache[i]
    n >>= 1
  end
  return s
end

def precheck(a,b,c)
  t = create_t(a,b,c)

  return nil unless regular_bin_p t
  return nil unless power_bin(t, (2**(unit_length * 4))) == t

  return [a,b,c]
end

def postcheck(a,b,c)
  t = create_t(a,b,c)
  clear_cache

  for d in @ds
    return nil if power_bin(t, d + 1) == t
  end

  return [a,b,c]
end

works = []
for a in 1...unit_length
  for b in 1...unit_length
    for c in 1...unit_length
      works << [a,b,c]
    end
  end
end

temp_results = []
results = []
temp_results = Parallel.map(works, progress:   ' Pre-check  ...') { |work|
  precheck work[0], work[1], work[2]
}.compact
results = Parallel.map(temp_results, progress: ' Post-check ...') { |work|
  postcheck work[0], work[1], work[2]
}.compact

puts "#{works.size} -> #{temp_results.size} -> #{results.size}"

s = ''
for result in results
  s << "#{result[0]}/#{result[1]}/#{result[2]}\n"
end

IO.write("xorshift-#{unit_length * 4}.log", s)

多倍長整数をネイティブに書きたかったのでrubyで組んでみたり、
そのくせ普通に書いたら遅すぎて使い物にならず高速化のためにぐちゃぐちゃ弄っていたり、
そもそもコメントがなかったりしますがご愛嬌。

計算時間とか

マシンスペック

CPU: i7-5820K (3.3 GHz, 12 thread)

rubyの自己診断

前段: 600時間23分2秒使い、250047パターンから有効な432パターンを抽出。
後段: 7時間25分29秒使い、432パターンから有効な152パターンを抽出。

timeによる診断

real: 36468分32.311秒
user: 435482分24.948秒
sys: 14分35.224秒

シングルスレッドで実行してたら10ヶ月掛かってたな・・・

まとめ

以下のことがわかった。

• A, B, C の最大値が最も小さいものは(1, 3, 10)の組み合わせ。
• A+B+C が最も小さくなるのも(1, 3, 10)の組み合わせ。