線形合同法のパラメータを乱数列から求めてみる。

線形合同法とは
とりあえず実装
c だけがわからないとき
a と c がわからないとき
全部わからないとき

線形合同法とは

一応文字の置き方の準備なども含め線形合同法について簡単に解説します。
線形合同法(Linear Congruential Generator, LCG)とは、以下の漸化式で定義されている疑似乱数生成アルゴリズムの一つです。

$\left\{ \begin{array}{l} x_{0} = seed\\ x_{i+1} = ax_{i} + c \mod M\\ \end{array} \right. \ \ \ \ (a,c,M,seed \in \mathbb{N})$

パラメータ $a,c,M$ の値によって周期や出てくる数値が代わり、最大周期は $M$ で $[0,\ M-1$ ] の自然数が出てきます(mod M なので)。また、 $x_{0}$ は初期値で $seed$ と呼ばれます。一般的にその時の時間やCPUの温度などを使って設定されます。
LCGは現在も広く使われ、その理由に実装が簡単なのと分かりやすい漸化式で定義されているので理解しやすいといったところがあると思います。しかし実は割と限界で、最大周期は $2^{32}$ と比較的短いです。その理由もあってか最近ではより高速でより周期も大きいメルセンヌツイスタアルゴリズムなどにとって代わられています。
後々使うので、ここでLCGの周期が最大になるようなパラメータの一つを紹介しておきます。

$\left\{ \begin{array}{l} a = 1664525\\ b = 1013904223\\ c = 2^{32}\\ \end{array} \right.$

他にもLCGの特徴は多くありますがそれらについてはググればたくさん出てくるのでこのあたりにしておきます。今回は最終的に、与えられた乱数列から各パラメータの値を求めてみます。わかりやすい日本語の記事がなかったのと自分の勉強がてら書こうと思います。

とりあえず実装

普通に線形合同法を実装してみます。漸化式の通り実装します。インライン引数に各パラメータと生成する乱数列の大きさをとります。

# LCG.py
import sys
import time

def LCG(seed, a, c, M, size):
    res = []
    x = seed

    for i in range(size):
        x = (a * x + c) % M
        res.append(x)

    return res

if __name__ == "__main__":
    a = int(sys.argv[1])
    c = int(sys.argv[2])
    M = int(sys.argv[3])
    size = int(sys.argv[4])

    res = LCG(time.time_ns(), a, c, M, size)

    print(res)

漸化式そのままですね。seedは適当に時間使ってます。気持ちはC言語のsrand(time(NULL))です。確かに、ランダムに数字が動いてるように見えます。
f:id:mi__tsu:20190310020837p:plain

ここで線形合同法で生成した既知の乱数列があり、その続きの値を予測することを考えます。ここではseedを適当な時間の値にしましたが、仮に $a,c,M$ すべての値がわかってるとしたらseedの値を既知の乱数列の最後の数値にしてあげそのまま生成するとその次の値、さらにその次の値と予測できます。漸化式を見ればわかりますが、一周期のうち取りうる値は一回ずつしか出てこないからですね。

c だけがわからないとき

連続した乱数列の大きさが2以上なら解けます。この場合は簡単で以下のような式変形で解くことができます。

$\begin{eqnarray} x_{1} &\equiv& ax_{0} + c \mod M\\ c &\equiv& x_{1} - ax_{0} \mod M \end{eqnarray}$

以下実装。引数に $a,M$ と標準入力で乱数列を渡します。

# getc.py
import sys

if __name__ == "__main__":
    a = int(sys.argv[1])
    M = int(sys.argv[2])

    x = []
    buf = input("number1 > ")
    x.append(int(buf))
    buf = input("number2 > ")
    x.append(int(buf))

    c = (x[1] - a * x[0]) % M

    print("c = {}".format(c))

確かに $c$ の値が復元できてますね。
f:id:mi__tsu:20190311171110p:plain

a と c がわからないとき

連続した乱数列の大きさが3以上なら解けます。値が二つわからないので一見複雑そうですが意外とそうでもありません。以下のような連立合同式を立ててあげます。

$\left\{ \begin{array}{l} x_{1} \equiv ax_{0} + c \mod M\\ x_{2} \equiv ax_{1} + c \mod M \end{array} \right.$

乱数の周期が2以上のとき次のように変形すると解け、 $a$ の値を求めることができます。

$\begin{eqnarray} x_{2} - x_{1} &\equiv& ax_{1} - ax_{0} \mod M\\ a &\equiv& \frac{x_{2} - x_{1}}{x_{1} - x_{0}} \mod M \end{eqnarray}$

のこりの $c$ については、先程の $c$ だけが分からないときの問題を解けば求めることができます。
また、周期が1のときは $a = 0$ となり $c$ の値も一般的に
$c = Mm + x_{0}\ \ \ (m \in \mathbb{N})$ と解けます。
実装で注意すべき点は有限体上の除算で、これを行うためにモジュラ逆数を用います。拡張ユークリッドの互助法(egcd)を使うことでmodの逆元(modinv)が実装できます。
以下実装。引数に $M$ と標準入力で乱数列をとります。周期が1の時は $c$ の最小の値を表示します。

# getac.py
import sys

def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    g, y, x = egcd(b % a, a)
    return (g, x - (b // a) * y, y)

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception("No modinv")
    return x % m

if __name__ == "__main__":
    M = int(sys.argv[1])

    x = []
    buf = input("number1 > ")
    x.append(int(buf))
    buf = input("number2 > ")
    x.append(int(buf))
    buf = input("number3 > ")
    x.append(int(buf))

    print(x)
    if x[0] == x[1]:
        a = 0
        c = x[0]
    else:
        a = (x[2] - x[1]) * modinv(x[1] - x[0], M) % M
        c = (x[2] - a * x[1]) % M

    print("a = {}".format(a))
    print("c = {}".format(c))

はい。
f:id:mi__tsu:20190311183354p:plain

全部わからないとき

乱数列が6つあればだいたい解けます。
今まで通り乱数列を増やして連立させてみることを考えます。すると以下のような式が立っていきます。

$\left\{ \begin{array}{l} x_{1} \equiv ax_{0} + c \mod M\\ x_{2} \equiv ax_{1} + c \mod M\\ x_{3} \equiv ax_{2} + c \mod M\\ x_{4} \equiv ax_{3} + c \mod M\\ \vdots \end{array} \right.$

上の式は次のように言い換えることもできます。

$\left\{ \begin{array}{l} x_{1} - (ax_{0} + c) = k_{1}M\\ x_{2} - (ax_{1} + c) = k_{2}M\\ x_{3} - (ax_{2} + c) = k_{3}M\\ x_{4} - (ax_{3} + c) = k_{4}M\\ \vdots \end{array} \right. \ \ \ \ (k_{n} \in \mathbb{N})$

式が一つ増えるごとに未知数 $k_{n}$ が増えます。
これを解決するには、ある自然数 $m$ に対し、ランダムないくつかの $m$ の倍数のgcdを取ると高い確率で $m$ になるという事実を用います。これを利用することで $x \equiv 0 \mod M$ となるような $M$ がgcdを取るだけで簡単に求まります。
ここで適当な配列 $T_{n}$ を以下のように定義しておきます。

$\begin{eqnarray} T_{0} &=& x_{1} - x_{0}\\ T_{1} &=& x_{2} - x_{1} \equiv (ax_{1} + c) - (ax_{0} + c) \equiv a(x_{1} - x_{0}) \equiv aT_{0} \mod M\\ T_{2} &=& x_{3} - x_{2} \equiv (ax_{2} + c) - (ax_{1} + c) \equiv a(x_{2} - x_{1}) \equiv aT_{1} \mod M\\ T_{3} &=& x_{4} - x_{3} \equiv (ax_{3} + c) - (ax_{2} + c) \equiv a(x_{3} - x_{2}) \equiv aT_{2} \mod M\\ \end{eqnarray}$

この $T_{n}$ を用いるとこのような変形が可能になります。

$T_{0}T_{2} - T_{1}^{2} \equiv a^{2}T_{0}^{2} - a^{2}T_{0}^{2} \equiv 0 \mod M$

これらを利用すれば $M$ の値を求めることができます。他のパラメータについても今までの方法を利用すれば問題なく求まるようになります。
以下実装。標準入力に乱数列をとります。

import sys
from functools import reduce
import math

def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    g, y, x = egcd(b % a, a)
    return (g, x - (b // a) * y, y)

def modinv(a, m):
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception("No modinv")
    return x % m

if __name__ == "__main__":

    # inputs
    x = []
    buf = input("number1 > ")
    x.append(int(buf))
    buf = input("number2 > ")
    x.append(int(buf))
    buf = input("number3 > ")
    x.append(int(buf))
    buf = input("number4 > ")
    x.append(int(buf))
    buf = input("number5 > ")
    x.append(int(buf))
    buf = input("number6 > ")
    x.append(int(buf))

    # solve M
    T = []
    for x0, x1 in zip(x, x[1:]):
        T.append(x1 - x0)

    zeros = []
    for t0, t1, t2 in zip(T, T[1:], T[2:]):
        zeros.append(t2 * t0 - t1 * t1)

    M = abs(reduce(math.gcd, zeros))

    # solve a, c
    if x[0] == x[1]:
        a = 0
        c = x[0]
    else:
        a = (x[2] - x[1]) * modinv(x[1] - x[0], M) % M
        c = (x[2] - a * x[1]) % M


    # answer
    print("a = {}".format(a))
    print("c = {}".format(c))
    print("M = {}".format(M))