離散対数問題 (Discrete Logarithm Problem) を任意 mod で解く

# 離散対数問題 (Discrete Logarithm Problem) を任意 mod で解く

離散対数問題 (Discrete Logarithm Problem) を解くアルゴリズムとして，Baby-step Giant-step が知られています．素数 mod の場合に関しては多くの記事が存在しますが，任意 mod を取り扱っている記事があまり見当たらなかったので，その辺りを書いていきます．

まず，離散対数問題とは以下のような問題を指します．


正整数 $M$ と $0 \leq X, Y \lt M$ なる整数 $X, Y$ が与えられる．
$$ X^K\equiv Y \pmod M \tag{1} $$ を満たす (最小の) 非負整数 $K$ を求めよ．


## M が素数の場合

まず，$X=0$ かどうかで場合分けをします．

$X=0$ のとき，$Y=1$ なら $K=0$，$Y=0$ なら $K=1$，$Y\geq 2$ なら $K$ は存在しません．ただし，ここでは $0^0=1$ と定義しています．

続いて，$X\neq 0$ のケースを考えます．$p=\left\lceil\sqrt{M}\right\rceil$ として，$K=p\cdot i+q\;(0\leq q\lt p)$ のように $i$, $q$ の組によって $K$ を表現します．Fermat の小定理より，調べるべき $i$ の範囲は $0\leq i\lt p$ となります．

いま，$M$ が素数かつ $0\lt X\lt M$ なので $\gcd(X,M)=1$ が成り立ち，
$$
X\cdot X'\equiv 1\pmod M,\quad 0\lt X'\lt M\tag{2}
$$
を満たす整数 $X'$ がただ一つ存在します．この $X'$ を $X^{-1}$ と書くことにします．$X^{-1}$ を用いると，$(1)$ 式は次のように変形できます．

```math
\begin{align}
X^K&\equiv Y\pmod M\\
X^{p\cdot i+q}&\equiv Y\pmod M\\
X^q&\equiv Y\cdot A^i\pmod M\tag{3}\\
\end{align}
```

ここで，$A=(X^{-1})^p$ としました．

$(3)$ の左辺がとりうる値は $X^0 \bmod M,X^1\bmod M,\dots,X^{p-1}\bmod M$ の高々 $p=\left\lceil\sqrt{M}\right\rceil$ 通りです．そこで，$X^j\bmod M$ を $j$ に写すような連想配列を前計算しておきます．最小の $K$ を求めたい場合は，$j$ の昇順に連想配列へと挿入し，上書きはしないようにします．

以上より，$(3)$ の右辺で $i=0,1,\dots,p-1$ を全て調べれば，$(1)$ を満たす $K$ が存在するならば目的の $i$, $q$ を得ることができます．このとき，$i$ の昇順に連想配列を検索すれば，初めて見つかった $i$, $q$ は $(1)$ を満たす最小の $K$ を構成します．どの $i$ に対しても対応する $j$ が存在しない場合は，$(1)$ を満たす $K$ は存在しません．

計算量については，連想配列の前計算と $i$ の全探索がボトルネックとなり，全体で $O(\sqrt{M})$ です．

## M が任意の場合

この記事の本題です．$M$ が合成数になると上のアルゴリズムをそのまま適用することは出来ません．というのは，$X^{-1}$ が存在しない可能性があるからです．しかし，$M$ が合成数の場合も ほとんど同じアルゴリズムで問題を解くことができます．

まず，$M=1$ のケースを考えます．このとき，任意の $X$, $Y$ に対して $K=0$ は $(1)$ を満たし，これは最小です．

以下，$M\gt 1$ とします．

$g=\gcd{(X,M)}$ とします．$g=1$ であれば，$X^{-1}$ が存在するので先に述べたアルゴリズムをそのまま適用することが出来ますが，$g\neq 1$ だとそうはいきません．

$X=g\cdot X'$，$M=g\cdot M'$ のように $X$, $M$ を分解します．このとき，$X'$ と $M'$ は互いに素であることに注意します．

これらを $(1)$ に代入すると，$(g\cdot X')^K\equiv Y \pmod{g\cdot M'}$ となります．従って，ある整数 $n$ が存在して，次の式 $(4)$ が成立します．
$$(g\cdot X')^K-Y=n\cdot g\cdot M'\tag{4}$$
$Y=1$ のとき，$K=0$ が求める答えです．また，$Y\neq 1$ のとき，$K=0$ は式 $(4)$ を満たしません．

$Y\neq 1$ のときは $K\neq 0$ なので，$Y$ は $g$ の倍数である必要があります．そこで，$Y$ も $Y=g\cdot Y'$ のように分解します．ここで，$Y$ が $g$ の倍数でなければ，$(1)$ を満たす $K$ は存在しないということになります．

$Y=g\cdot Y'$ を用いて $(4)$ 式の計算を進めることで，式 $(5)$ を得ます．

```math
\begin{align}
gX'\cdot(gX')^{K-1}-gY'&=n\cdot gM'\\
X'\cdot (gX')^{K-1}-Y'&=nM'\\
X'\cdot (gX')^{K-1}&\equiv Y'\pmod{M'}\tag{5}
\end{align}
```

ここで，$\gcd{(X',M')}=1$ より $X'\cdot{X''}\equiv 1 \pmod{M'}$，$0\lt X''\lt M'$ なる $X''$ がただ一つ存在します．これを $X'^{-1}$ とおき，$(5)$ 式に代入することで $(6)$ 式を得ます．

```math
\begin{align}
X'\cdot(gX')^{K-1}&\equiv Y'\pmod{M'}\\
(gX')^{K-1}&\equiv Y'\cdot{X'}^{-1}\pmod{M'}\\
X^L&\equiv Z\pmod{M'}\tag{5}\\
\end{align}
```

2 行目から 3 行目の変形では $L=K-1,\ Z=Y'\cdot {X'}^{-1}$ と置きました．

$(5)$ は元の問題よりもサイズの小さい離散対数問題になっているので，再帰的にこの問題を解くことが出来ます．$(5)$ を満たす最小の非負整数 $L$ が存在すれば $K=L+1$ が $(1)$ の最小解であり，$L$ が存在しなければ $K$ も存在しません．

ここで，一つ注意すべきなのは $\gcd{(X,M')}=1$ とは限らないということです．

例えば $(X,M)=(20,24)$ とすると，$g=\gcd{(X, M)}=4$ なので，$M'=M/g=6$ です．この時，$\gcd{(X,M')}=2$ です．

---

この解法全体の計算量は $\gcd{(X,M)}=1$ となるまでの再帰の計算量と再帰の終端における Baby-step Giant-step の計算量を合わせたものになります．

再帰するのは $\gcd{(X,M)}\neq 1$ の場合で，このとき $(X, M)$ は $(X, M/\gcd{(X, M)})$ へと移ります．従って，再帰の度に $M$ は半分以下になるので再帰の深さは $O(\log M)$ で抑えられます．再帰の度にユークリッドの互除法を行うので，再帰パートの計算量は $O((\log M)^2)$ となります．

Baby-step Giant-step の計算量は先に述べた通り $O(\sqrt{M})$ なので，全体の計算量は $O(\sqrt M)$ です．

### 再帰の除去について (追記: 2021/02/27)

上の議論からもう少し考えてみると，再帰を行わずに任意 mod で離散対数問題が解けることが分かります．

再帰する必要があるのは $\gcd{(X, M/\gcd{(X,M)})}\neq 1$ となるケース，即ち，$X$ と $M$ の共通の素因数 $p$ であって，$p$ が $M$ を割り切る回数が $X$ を割り切る回数よりも多いようなものが存在する場合です．従って，合同式を割る値 $g$ として，もう少し大きな $M$ の約数を持ってくることができれば良さそうだということになります．

そこで，突然ですが，ある正整数 $d$ が存在して，$K=0,\ldots,d-1$ は式 $(1)$ を満たさないと仮定します．このとき，$(1)$ を満たす $K$ が存在するならば
$$
X^d\cdot X^t\equiv Y\pmod M
$$
を満たす非負整数 $t$ が存在します．そして，$g=\gcd{(X^d,M)}$ とおいて同様の式変形を行うことで
$$
X^t\equiv (Y/g)\cdot\left(X^d/g\right)^{-1}\pmod{M/g}
$$
を得ます ($Y$ が $g$ で割り切れる必要があることに注意してください)．$d$ が十分大きければ $\gcd{(X,M/g)}=1$ となり，Baby-step Giant-step を適用して $t$ を求めることが出来るため，再帰することなく問題が解けることになります．

では，「$d$ が十分大きければ」とは具体的にどれくらいあれば十分でしょうか？ $X^d$ と $M$ に共通する全ての素因数 $p$ に対して $X^d$ を割り切る回数が $M$ を割り切る回数以上であればよいので，$M$ の素因数分解における指数の最大値以上あれば十分ということになります．そして，この最大値は $\lfloor\log_2 M\rfloor$ で抑えられます．

以上より，次のアルゴリズムで任意 mod における離散対数問題を解くことができ，全体の計算量は同じく $O(\sqrt{M})$ です．

1. $d=\lfloor\log_2 M\rfloor$ を求める．$K=0,\ldots,d-1$ に対して式 $(1)$ を満たすかを確認し，満たすものが存在すればそれを出力して終了する．存在しない場合は，2 へ進む．
2. $g=\gcd{(X^d,M)}$ を計算し，式 $(6)$ を満たす $t$ を Baby-step Giant-step で求めて $d+t$ を出力する．ただし，$t$ が存在しない場合や $Y$ が $g$ で割り切れない場合は式 $(1)$ を満たす $K$ は存在しない．

```math
\begin{align}
X^t&\equiv (Y/g)\cdot\left(X^d/g\right)^{-1}\pmod{M/g}\\
&\equiv Y\cdot\left(X^d\right)^{-1}\pmod{M/g}\tag{6}
\end{align}
```

## 実装例

<details><summary>再帰解法 (Java)</summary><div>

```java:DiscreteLogarithm.java
public class DiscreteLogarithm {
    /**
     * 任意 mod で離散対数問題を解く．つまり，x^k=y mod m を満たす最小の非負整数 k を求める．
     * @return x^k=y mod m を満たす最小の非負整数 k．存在しない場合は -1 を返す．
     */
    public static long log(long x, long y, long m) {
        if ((x %= m) < 0) x += m;
        if ((y %= m) < 0) y += m;
        // m = 1 のケースを処理
        if (m == 1) return 0;
        // y = 1 (k = 0) のケースを処理
        if (y == 1) return 0;
        long[] gi = gcdInv(x, m);
        // g=gcd(x,m), ix*x=g mod m
        long g = gi[0], ix = gi[1];
        if (g != 1) {
            // y は gcd(x, g) の倍数である必要がある
            if (y % g != 0) return -1;
            m /= g;
            y = (y / g) * ix;
            // 小さいサイズの問題を解く
            long log = log(x, y, m);
            // 小さいサイズの問題に解が存在しなければ，元の問題にも解は存在しない
            if (log < 0) return -1;
            // 小さいサイズの問題の答えを l とすると，元の問題の答えは l + 1
            return log + 1;
        }
        // gcd(x, g) = 1 であれば Baby-step Giant-step で解ける
        return coprimeLog(x, y, m, ix);
    }
    /**
     * x と m が互いに素な場合に，Baby-step Giant-step により離散対数問題を解く．
     * @param ix mod {@code m} での {@code x} の乗法逆元 x^{-1}
     * @return x^k=y mod m を満たす最小の非負整数 k．存在しない場合は -1 を返す．
     */
    private static long coprimeLog(long x, long y, long m, long ix) {
        long p = ceilSqrt(m);
        java.util.HashMap<Long, Long> memo = new java.util.HashMap<>();
        for (long i = 0, powx = 1; i < p; i++) {
            memo.putIfAbsent(powx, i);
            powx = powx * x % m;
        }
        long a = pow(ix, p, m);
        for (long i = 0, ypowa = y; i < p; i++) {
            long q = memo.getOrDefault(ypowa, -1L);
            if (q >= 0) return p * i + q;
            ypowa = ypowa * a % m;
        }
        return -1L;
    }
    /**
     * (x-1)^2 < n <= x^2 を満たす x を二分探索により求める．
     * @return (x-1)^2 < n <= x^2 を満たす x
     */
    private static long ceilSqrt(long n) {
        long l = -1, r = n;
        while (r - l > 1) {
            long m = (r + l) >> 1;
            if (m * m >= n) r = m;
            else l = m;
        }
        return r;
    }
    /**
     * a^b mod m を二分累乗法により計算する．
     * @return a^b mod m
     */
    private static long pow(long a, long b, long m) {
        if (m == 1) return 0;
        if ((a %= m) < 0) a += m;
        long pow = 1;
        for (long p = a, c = 1; b > 0;) {
            long lsb = b & -b;
            while (lsb != c) {
                c <<= 1;
                p = (p * p) % m;
            }
            pow = (pow * p) % m;
            b ^= lsb;
        }
        return pow;
    }
    /**
     * g = gcd(x, m), a * x = g mod m を満たす組 (g, x) をユークリッドの互除法により計算する．
     * @param a
     * @param m mod
     * @return g = gcd(x, m), a * x = g mod m を満たす組 (g, x)
     */
    private static long[] gcdInv(long a, long m) {
        if ((a %= m) < 0) a += m;
        if (a == 0) return new long[]{m, 0};
        long s = m, t = a;
        long m0 = 0, m1 = 1;
        while (t > 0) {
            long u = s / t;
            s -= t * u;
            m0 -= m1 * u;
            long tmp;
            tmp = s; s = t; t = tmp;
            tmp = m0; m0 = m1; m1 = tmp;
        }
        if (m0 < 0) m0 += m / s;
        return new long[]{s, m0};
    }
}
```

</div></details>

[問題リンク](https://judge.yosupo.jp/problem/discrete_logarithm_mod)
[提出 URL](https://judge.yosupo.jp/submission/23711)

<details><summary>非再帰解法 (Python)</summary><div>

```Python
import math

def ceil_sqrt(n):
    l, r = -1, n
    while r - l > 1:
        m = (l + r) >> 1
        if m * m >= n:
            r = m
        else:
            l = m
    return r

def discrete_logarithm_coprime_mod(x, y, m):
    x %= m
    y %= m
    if y == 1 or m == 1:
        return 0
    if x == 0:
        if y == 0:
            return 1
        else:
            return None
    p = ceil_sqrt(m)
    a = pow(x, -p, m)
    mp = {}
    pow_x = 1
    for j in range(p):
        mp.setdefault(pow_x, j)
        pow_x = pow_x * x % m
    ya = y
    for i in range(p):
        if ya in mp:
            j = mp[ya]
            return p * i + j
        ya = ya * a % m
    return None


def discrete_logarithm_arbitrary_mod(x, y, m):
    if m == 1:
        return 0
    x %= m
    y %= m
    d = len(format(m, 'b'))
    pow_x = 1
    for i in range(d):
        if pow_x == y:
            return i
        pow_x = pow_x * x % m
    g = math.gcd(pow_x, m)
    if y % g != 0:
        return None
    m //= g
    z = y * pow(pow_x, -1, m)
    t = discrete_logarithm_coprime_mod(x, z, m)
    return None if t is None else d + t


if __name__ == '__main__':
    t = int(input())
    for _ in range(t):
        x, y, m = map(int, input().split())
        k = discrete_logarithm_arbitrary_mod(x, y, m)
        print(-1 if k is None else k)
```

</div></details>

[問題リンク](https://judge.yosupo.jp/problem/discrete_logarithm_mod)
[提出 URL](https://judge.yosupo.jp/submission/40389)