底とmodが互いに素でない離散対数問題

離散対数問題とは

整数 p \ge 1, \ 0 \le a, b \lt p に対して, a^{x} \equiv b \mod p を満たす最小の非負整数 x を求める問題です.

例えば p=40, \ a = 2, \ b = 24 のとき,

  • a^{0} \equiv 1
  • a^{1} \equiv 2
  • a^{2} \equiv 4
  • a^{3} \equiv 8
  • a^{4} \equiv 16
  • a^{5} \equiv 32
  • a^{6} \equiv 24

より x=6 が解となります.

Baby-step Giant-stepアルゴリズム

底 a と mod p が互いに素なときに使える離散対数問題の O(\sqrt{p}) 解法です.

h := \lceil \sqrt{p} \rceil とします. x = ih - j \ (0 \lt i, j \le h) とおくと, a^{x} \equiv b \iff a^{ih} \equiv a^{j}b となります.

したがって a^{1}b, \dots, a^{h}b を hash table に保存しておいて, a^{h}, a^{2h}, \dots, a^{h^{2}} の順に存在判定をすれば良いです.

底とmodが互いに素でない場合

p の素因数分解を p = p_1^{d_1} \cdots p_m^{d_m} \ q_1^{e_1} \cdots q_n^{e_n} とおきます. ただし p_1, \dots, p_m は a の素因数であり, q_1, \dots, q_n は a の素因数ではないとします.

このとき, x \ge \lfloor \log_2 p \rfloor に対して a^{x} は g := p_1^{d_1} \dots p_m^{d_m} の倍数になります. したがって,

(1). b が g の倍数でないとき

x \ge \lfloor \log_2 p \rfloor の範囲で a^{x} \equiv b \mod p の解は存在しません.

(2). b が g の倍数であるとき

x \ge \lfloor \log_2 p \rfloor に対して

a^{x} \equiv b \mod p \iff a^{x} / g \equiv b / g \mod p / g

が成立します. a と p/g = q_1^{e_1} \cdots q_n^{e_n} は互いに素なので, これは Baby-step Giant-step アルゴリズムで解くことができます.

全体の流れ

  1. 0 \le x \lt \lfloor \log_2 p \rfloor については愚直に a^{x} \equiv b \mod p かどうかを判定する. 解が見つかれば終了.
  2. b が g = p_1^{d_1} \dots p_m^{d_m} の倍数でないなら, 解なし.
  3. a^{x} / g \equiv b / g \mod p / g を Baby-step Giant-step アルゴリズムで解く.

実装解説

離散対数問題(Mod-Log) | Luzhiled’s memo の実装を, 上記の説明に照らし合わせて解説します.

(※ 元のコードを一部変更しております.)

int64_t mod_log(int64_t a, int64_t b, int64_t p) {
  // g := p_1^{d_1} ... p_m^{d_m} を求める.
  int64_t g = 1;
  for(int64_t i = p; i; i /= 2) (g *= a) %= p;
  g = __gcd(g, p);

  // a^{x} が g の倍数になる手前までは愚直に計算する.
  int64_t t = 1, c = 0;
  for(; t % g; c++) {
    if(t == b) return c;
    (t *= a) %= p;
  }
  // a^{c} = t

  // b が g の倍数でないなら解なし.
  if(b % g) return -1;

  // a^{y} (t / g) = b / g mod p / g を解いて, y + c を解とする.
  t /= g;
  b /= g;
  p /= g;
  // a^{y} t = b mod p を解く.

  int64_t h = 0, gs = 1;

  // h := ceil(sqrt(p))
  // gs := a^{h}
  for(; h * h < p; h++) (gs *= a) %= p;

  // a^{ih-j} t = b mod p (0 < i, j <= h)
  // (a^{h})^{i} t = b a^{j} mod p
  // (gs)^{i} t = b a^{j} mod p

  // b a^{j} (0 < j <= h) を bs に登録.
  unordered_map< int64_t, int64_t > bs;
  for(int64_t s = 0, e = b; s < h; bs[e] = ++s) {
    (e *= a) %= p;
  }

  // (gs)^{i} t (0 < i <= h) で bs を検索. y = s - bs[e]
  for(int64_t s = 0, e = t; s < p;) {
    (e *= gs) %= p;
    s += h;
    if(bs.count(e)) return c + s - bs[e];
  }

  // 見つからなければ解なし.
  return -1;
}

参考