DivineJK’s diary

Z cjikf c evf'

群論と離散対数のお時間です - 2

前回の続きです。一般化されたBSGSはまとまりをよくするために次回に回し、今回は、

  • Pohlig-Hellman algorithm
  • BSGS, Pohlig-Hellman algorithmの実装

について書いていきたいと思います。

Pohlig-Hellman algorithm

 \varphi(M)素因数分解したときの指数が大きい場合に有効なアルゴリズムであるPohlig-Hellman algorithmを紹介します。ようやくここで群論が登場します。

ここでも M Xが互いに素である必要があります。ところで、これは X^{\varphi(M)}\equiv 1\bmod Mが成り立つことを意味しますが、 \varphi(M) X^{k}\equiv 1\bmod Mを満たすものの中で最小ではない可能性があります。しかし、この k \varphi(M)の約数であることが示せます。そこで、 \varphi(M)の約数のうちで等式を満たすような最小なものを探しておきます。そして、これが巡回群の位数となります。

以下、この巡回群の位数を nとします。そして、 n = p_1^{e_1}\times p_2^{e_2}\times\cdots\times p_m^{e_m}のように素因数分解できるとします。ここで注意するべきなのは、今考えているのはあくまで \bmod Mでの話題なので、 nではなく Mでmodをとらなくてはなりません(私が実装でハマった原因はこれでした)。

まず、Pohlig-Hellman algorithmの特別な場合の部分問題として、 n = p^{e}のように、素因数分解してひとつの素数のみのべき乗の形になっている場合を考えます。解く問題としては、 X^{K}\equiv Y\bmod Mなる 0\le K\lt p^{e}を求める、という問題です。位数が p^{e}なので、解があるとすればこの範囲に必ずおさまることが保証されます。この Kは、 K=a_0 + a_1 p +\cdots +a_ {e-1} p^{e-1} \, (0\le a_0 ,a_1 ,\cdots , a _ {e-1} \lt p)と表せるので、この a_0, a_1, \cdots, a _ {e-1}を求めたいです。

 a_0, a_1, \cdots, a _ {k-1}までが求まっているとして、 R = a_0 +a_1p+ \cdots + a _ {k-1}p^{k-1}とします(何も求まっていないときは R = 0とします)。このとき、方程式を変形して X^{a _ k p^{k} + \cdots + a _ {e-1} p^{e-1}} \equiv Y(X^{-1})^{R}\bmod Mを得ることができます。この両辺を p^{e-k-1}乗します。今考えている巡回群の位数は p^{e}ですから、 X^{cp^{e}} \equiv 1\bmod M\quad (c=0,1,2,\cdots)となります。よって、左辺は X^{a _ k p^{e-1}} = (X^{p^{e-1}})^{a _ k}のみ残ります。

ここで、 \gamma = X^{p^{e-1}}とおくことで、 \gamma^{a _ k} \equiv (YX^{-R})^{p^{e-k-1}}\bmod Mです。 \gammaの生成する巡回群の位数は pであるので、これを満たす a _ kをBSGSで求めることができます。もし解が存在しなければその旨を報告し、終了です。そうでないとき、 R a _ kp^{k}を足し、 a _ {k+1}を求める計算に進みます。

この操作を k=0から順に行い、 K=Rが求める答えになります。 0\le k\lt eで、各 kあたりの計算量は \mathrm{O}(\sqrt{p})なので、全体的な計算量は \mathrm{O}(e\sqrt{p})です。

もとの問題に戻ります。位数が n = p_1^{e_1}\times p_2^{e_2}\times\cdots\times p_m^{e_m}と表され、 X _ i = X^{n/p _ i^{e _ i}} p _ i^{e _ i}巡回群を生成します。同時に、 Y _ i = Y^{n/p _ i^{e _ i}}とします。問題の方程式の両辺を n/p _ i^{e _ i}乗することで、方程式 X _ i^{K} \equiv Y _ i\bmod Mを得ます。この方程式は先程の部分問題の解き方で解くことができます。また、この方程式の解を K _ iとすると、これは K \bmod p _ i^{e _ i}での値になります。

そうして得られた位数 p _ i^{e _ i}での方程式の解 K _ iをもとに、中国剰余定理から Kを求めることができます。

全体の計算量は \mathrm{O}(\sum_i e _ i(\sqrt{p _ i} + \log n))です。

実装例

BSGSの実装例

def divisors(n): # nの約数
    L = []
    p = 1
    while p * p <= n:
        if n % p == 0:
            L.append(p)
        p += 1
    for i in range(len(L), 0, -1):
        if L[i-1] * L[i-1] != n:
            L.append(n // L[i-1])
    return L # nの約数が昇順に並んでいる。
def totient(n): # nと互いに素なn以下の正整数の個数
    a = n
    p = 2
    res = n
    if a % p == 0:
        res >>= 1
        while 1 - a % 2:
            a >>= 1
    p += 1
    while a != 1:
        if a % p == 0:
            res //= p
            res *= p - 1
            while a % p == 0:
                a //= p
        p += 2
        if p * p > n and a != 1:
            res //= a
            res *= a - 1
            break
    return res
def babystep_giantstep(n, a, b): # a^k = b (mod n)なる最小の非負整数kを求める(存在しなければ-1)。
    baby_step = {}
    candidate = divisors(totient(n)) # 巡回群の位数の候補
    pnt = 0
    while pow(a, candidate[pnt], n) != 1: # a^k = 1 (mod n)となる最小のkを探す。それが巡回群の位数となる。
        pnt += 1
        if pnt == len(candidate):
            return -1 # nとaが互いに素でない(一般化した時に扱う)。
    p = candidate[pnt] # 巡回群の位数
    l, r = 0, p
    m = (l + r) // 2 # m = ceil(sqrt(p))
    while r - l > 1:
        if m * m <= p:
            l = m
        else:
            r = m
        m = (l + r) // 2
    if m * m < p:
        m += 1
    x, y, u, v, k, l = 1, 0, 0, 1, a, n
    while l:
        x, y, u, v, k, l = u, v, x - u * (k // l), y - v * (k // l), l, k % l
    bas = pow(x, m, n) # bas = a^(-m) (mod n)
    e = 1
    for i in range(m): # baby_stepに(a^i, i)を記録していく。
        baby_step[e] = i
        e *= a
        e %= n
    y = b
    for i in range(m):
        if y in baby_step:
            return i * m + baby_step[y] # これが最小の解
        y *= bas # giantstep
        y %= n
    return -1 # 見つからなければ-1を返す。

Pohlig-Hellman algorithmの実装例

def divisors(n): # nの約数
    L = []
    p = 1
    while p * p <= n:
        if n % p == 0:
            L.append(p)
        p += 1
    for i in range(len(L), 0, -1):
        if L[i-1] * L[i-1] != n:
            L.append(n // L[i-1])
    return L # nの約数が昇順に並んでいる。
def totient(n): # nと互いに素なn以下の正整数の個数
    a = n
    p = 2
    res = n
    if a % p == 0:
        res >>= 1
        while 1 - a % 2:
            a >>= 1
    p += 1
    while a != 1:
        if a % p == 0:
            res //= p
            res *= p - 1
            while a % p == 0:
                a //= p
        p += 2
        if p * p > n and a != 1:
            res //= a
            res *= a - 1
            break
    return res
def babystep_giantstep(n, a, b): # a^k = b (mod n)なる最小の非負整数kを求める(存在しなければ-1)。
    baby_step = {}
    candidate = divisors(n-1) # 巡回群の位数の候補
    pnt = 0
    while pow(a, candidate[pnt], n) != 1: # a^k = 1 (mod n)となる最小のkを探す。それが巡回群の位数となる。
        pnt += 1
        if pnt == len(candidate):
            return -1 # nとaが互いに素でない(一般化した時に扱う)。
    p = candidate[pnt] # 巡回群の位数
    l, r = 0, p
    m = (l + r) // 2 # m = ceil(sqrt(p))
    while r - l > 1:
        if m * m <= p:
            l = m
        else:
            r = m
        m = (l + r) // 2
    if m * m < p:
        m += 1
    x, y, u, v, k, l = 1, 0, 0, 1, a, n
    while l:
        x, y, u, v, k, l = u, v, x - u * (k // l), y - v * (k // l), l, k % l
    bas = pow(x, m, n) # bas = a^(-m) (mod n)
    e = 1
    for i in range(m): # baby_stepに(a^i, i)を記録していく。同じkeyは小さい方をとる。
        if e not in baby_step:
            baby_step[e] = i
        e *= a
        e %= n
    y = b
    for i in range(m):
        if y in baby_step:
            return i * m + baby_step[y] # これが最小の解
        y *= bas # giantstep
        y %= n
    return -1 # 見つからなければ-1を返す。
def extgcd(a, b, c): # ax + by = cを満たす整数の組(x, y)を探す。
    if b == 0:
        if a == 0:
            if c == 0:
                return 0, 0
            return None
        if c % a == 0:
            return c // a, 0
        return None
    if b < 0:
        a, b, c = -a, -b, -c
    tk, tl = a, b
    while tl:
        tk, tl = tl, tk % tl
    if c % tk:
        return None
    a //= tk
    b //= tk
    c //= tk
    x, y, u, v, k, l = 1, 0, 0, 1, a, b
    while l:
        x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
        k, l = l, k % l
    x *= c
    k = x // b
    x -= k * b
    y *= c
    y += k * a
    return x, y
def CRT(num, a_list, m_list): # 複数modの条件を満たす整数を見つける。
    r = 0
    bas = 1
    x, y = 0, 0
    for i in range(num):
        x, y = extgcd(bas, -m_list[i], a_list[i]-r)
        r += bas * x
        bas *= m_list[i]
    return r % bas
def pohlig_hellman(mod, a, b):
    candidate = divisors(totient(mod)) # 巡回群の位数の候補
    pnt = 0
    while pow(a, candidate[pnt], mod) != 1: # a^k = 1 (mod)となる最小のkを探す。それが巡回群の位数となる。
        pnt += 1
        if pnt == len(candidate):
            return -1 # nとaが互いに素でない(一般化した時に扱う)。
    n = candidate[pnt] # 巡回群の位数
    factors = [] # nの素因数分解
    res_seed = [] # 素数のべき乗での答えを格納
    mod_list = []
    r = 0
    cnt = 0
    tmp = n
    # 素因数分解パート
    while tmp % 2 == 0:
        cnt += 1
        tmp >>= 1
    if cnt:
        r += 1
        res_seed.append(0)
        mod_list.append(0)
        factors.append((2, cnt))
    pr = 3
    while tmp != 1:
        cnt = 0
        while tmp % pr == 0:
            tmp //= pr
            cnt += 1
        if cnt:
            r += 1
            res_seed.append(0)
            mod_list.append(0)
            factors.append((pr, cnt))
        pr += 2
        if pr * pr > n and tmp != 1:
            res_seed.append(0)
            mod_list.append(0)
            r += 1
            factors.append((tmp, 1))
            break
    # 素数のべき乗についての部分問題を解く。
    for i in range(r):
        p = factors[i][0]
        e = factors[i][1]
        p_i = pow(p, e)
        mod_list[i] = p_i
        g = pow(a, n//p_i, mod) # gは位数p^eの群をなす。
        h = pow(b, n//p_i, mod)
        # g^res_seed[i] = h (mod) となるres_seed[i]を求める。
        # res_seed[i] = x_0 + x_1 * p + ... + x_{e-1} * p^{e-1} (0 <= x_j < p)の形で表せる。
        bas = 1
        inv_g = extgcd(g, mod, 1)[0] # gの逆元
        gamma = pow(g, pow(p, e-1), mod) # これは位数pの巡回群をなす。
        for k in range(e):
            # (x_0, x_1, .... x_{k-1})までわかっているとする。
            # 式変形により、g^(x_k*p^k+...+x_{e-1}*p^{e-1}) = h * g^(-(x_0+x_1*p+...+x_{k-1}*p^{k-1}))
            # 両辺p^{e-k-1}することで、gamma^(x_k) = ...の形になるので、x_kをbsgsで求める。
            hk = pow(h*pow(inv_g, res_seed[i], mod)%mod, pow(p, e-k-1), mod)
            xk = babystep_giantstep(mod, gamma, hk)
            if xk < 0:
                return -1 # x_kが見つからなければ-1
            res_seed[i] += xk * bas
            bas *= p
    # CRTでそれぞれの部分問題の答えをまとめる。
    res = CRT(r, res_seed, mod_list)
    return res

verify

atcoder.jp  X Pが固定なので、辞書や位数を一度計算してメモしておくように実装すると、計算量を減らすことができます。