群論と離散対数のお時間です - 2
前回の続きです。一般化されたBSGSはまとまりをよくするために次回に回し、今回は、
- Pohlig-Hellman algorithm
- BSGS, Pohlig-Hellman algorithmの実装
について書いていきたいと思います。
Pohlig-Hellman algorithm
を素因数分解したときの指数が大きい場合に有効なアルゴリズムであるPohlig-Hellman algorithmを紹介します。ようやくここで群論が登場します。
ここでもとが互いに素である必要があります。ところで、これはが成り立つことを意味しますが、はを満たすものの中で最小ではない可能性があります。しかし、このはの約数であることが示せます。そこで、の約数のうちで等式を満たすような最小なものを探しておきます。そして、これが巡回群の位数となります。
以下、この巡回群の位数をとします。そして、のように素因数分解できるとします。ここで注意するべきなのは、今考えているのはあくまででの話題なので、ではなくでmodをとらなくてはなりません(私が実装でハマった原因はこれでした)。
まず、Pohlig-Hellman algorithmの特別な場合の部分問題として、のように、素因数分解してひとつの素数のみのべき乗の形になっている場合を考えます。解く問題としては、なるを求める、という問題です。位数がなので、解があるとすればこの範囲に必ずおさまることが保証されます。このは、と表せるので、このを求めたいです。
までが求まっているとして、とします(何も求まっていないときはとします)。このとき、方程式を変形してを得ることができます。この両辺を乗します。今考えている巡回群の位数はですから、となります。よって、左辺はのみ残ります。
ここで、とおくことで、です。の生成する巡回群の位数はであるので、これを満たすをBSGSで求めることができます。もし解が存在しなければその旨を報告し、終了です。そうでないとき、にを足し、を求める計算に進みます。
この操作をから順に行い、が求める答えになります。で、各あたりの計算量はなので、全体的な計算量はです。
もとの問題に戻ります。位数がと表され、はの巡回群を生成します。同時に、とします。問題の方程式の両辺を乗することで、方程式を得ます。この方程式は先程の部分問題の解き方で解くことができます。また、この方程式の解をとすると、これはのでの値になります。
そうして得られた位数での方程式の解をもとに、中国剰余定理からを求めることができます。
全体の計算量はです。
実装例
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 とが固定なので、辞書や位数を一度計算してメモしておくように実装すると、計算量を減らすことができます。