線形合同法であるseedが0からいくつ進めたものかを得る その3 (法が2のべきとは限らない一般の場合)

• 線形合同法であるseedが0からいくつ進めたものかを得る - oupoの日記
• 線形合同法であるseedが0からいくつ進めたものかを得る その2 - oupoの日記

法が2のべきとは限らない一般の場合で求める方法ができた。

法が素数pのべきの場合、{x_n}が法p^eで最大周期のLCG数列なら{x_{pn} / p}が法p^{e-1}で最大周期のLCGの数列になることを使う。
一般の法mに対しては、mを素因数分解したのち、各素数について計算し、中国剰余定理で結果を得る。

require "prime"

# プログラムの結果が正しいかテストする
def test()
  m = 2 * (3 ** 4) # テストする法mの値

  factors = Prime.prime_division(m)
  product = factors.map{|p,e| p }.inject(1, :*)
  if m % 4 == 0 then product *= 2 end

  # 最大周期になる全てのaとbを試す
  (0..m/product-1).each do |k|
    a = product * k + 1
    m.times do |b|
      next if gcd(b, m) != 1
      list = gen_list(a, b, m)
      list.each_with_index do |s, i|
        index = calc_index(s, a, b, m)
        if i != index
          raise "wrong result (#{s},#{a},#{b},#{m}) expected: #{i}, got: #{index}"
        end
      end
      puts "ok a=#{a}, b=#{b}"
    end
  end
end

# x_0 = 0, x_{n+1} = (ax_n + b) mod m で定義される数列の最初のm項を得る
def gen_list(a, b, m)
  list = []
  s = 0
  m.times do
    list << s
    s = (a * s + b) % m
  end
  list
end

# x_0 = 0, x_{n+1} = (ax_n + b) mod m で定義される数列でx_n = sとなるnを求める
# ただし数列の周期はmでなければならない
def calc_index(s, a, b, m)
  factors = Prime.prime_division(m)
  res = 0
  product = 1
  factors.each do |p, e|
    i = calc_index_pe(s, a, b, p, e)
    res = chinese(res, product, i, p**e)
    product *= p**e
  end
  res
end

# calc_indexでm = p^e (pは素数)の場合
def calc_index_pe(s, a, b, p, e)
  if e == 0 then
    0
  else
    l = (inverse(b, p) * s) % p
    pe = p ** e
    ss = step(a, b, s, p - l, pe)
    (aa, bb) = step_param(a, b, p, pe)
    (calc_index_pe(ss / p, aa, bb / p, p, e - 1) * p - p + l) % pe
  end
end

# f(x) = (ax + b) mod mに対してf^n(s)を求める
def step(a, b, s, n, m)
  if n == 0 then
    s
  elsif n % 2 == 0 then
    step((a * a) % m, ((a + 1) * b) % m, s, n / 2, m)
  else
    (step((a * a) % m, ((a + 1) * b) % m, s, n / 2, m) * a + b) % m
  end
end

# 関数 f(x) = (ax + b) mod m のn個の合成f^nは
# f^n(x) = (a'x + b') mod mと書ける。そのa', b'を得る
def step_param(a, b, n, m)
  if n == 0
    [1, 0]
  elsif n % 2 == 0
    step_param((a * a) % m, (a + 1) * b % m, n / 2, m)
  else
    (aa, bb) = step_param((a * a) % m, (a + 1) * b % m, n / 2, m)
    [(aa * a) % m, (aa * b + bb) % m]
  end
end

# 中国剰余定理
# x ≡ a (mod p), x ≡ b (mod q) となるxを求める (p, qは互いに素)
def chinese(a, p, b, q)
  (a * q * inverse(q, p) + b * p * inverse(p, q)) % (p*q)
end

# ax + by = k の解を得る (ただしkはaとbの最大公約数)
def extgcd(a, b)
  if b == 0 then
    return [1, 0]
  end
  (q, r) = a.divmod(b)
  (x, y) = extgcd(b, r)

  # 今 bx + ry = k が成り立っている。r = a - bqを使えば
  # ay + b(x - qy) = kを得る
  [y, x - q * y]
end

def gcd(a, b)
  if b == 0 then
    a
  else
    gcd(b, a % b)
  end
end

# a x ≡ 1 (mod m) となるxを得る (aとmは互いに素)
def inverse(a, m)
  (x, y) = extgcd(a, m)
  x
end