Trait nekolib::math::DLog

pub trait DLog: Sized {
    // Required method
    fn dlog(self, a: Self, n: Self) -> Option<Self>;
}

離散対数。

$n$ を法とする $\mathrm{dlog}_b(a)$。 $b^z \equiv a \pmod{n}$ なる $z\ge 0$ が存在すれば、そのうち最小のものを返す。 $0\le b, a\lt n$ とする。 また、$0^0 = 1$ とする。

コーナーケースに注意。 $n=1$ なら $0^0=1\equiv 0$ より $z=0$ を返す。 以下、$n\gt 1$ とする。

• $\mathrm{dlog}_\bullet(1) = 0$
• $\mathrm{dlog}_0(0) = 1$

上記に該当しないとき、$\mathrm{dlog}_0(\bullet)$ と $\mathrm{dlog}_1(\bullet)$ は存在しない。残りのケースについては、以下の方針に基づいて求める。

Idea

$1\to b\to b^2\to\dots\pmod{n}$ でできる functional graph を考えると、$\rho$ の字状になっていることに気づく。しっぽの長さ $\mu$ は高々 $\log(n)$、頭の長さ $\lambda$ は $\varphi(n)$ の約数になっていることが示せる。

そこで、まずしっぽの部分に解があるかを愚直に $O(\mu)$ 時間で調べる。 見つからなければ、頭の部分に解があるかを $o(\lambda)$ 時間で調べる。 $O(\sqrt{\lambda}\log(\lambda))$ 時間のアルゴリズムを紹介する。

$\rho$ の頭の部分は $b^\mu\to b^{\mu+1}\to\dots\to b^{\mu+\lambda}\equiv b^\mu$ である。このとき、${}^\forall i\ge\mu$ について $b^{i+1}\cdot b^{\lambda-1} \equiv b^i \pmod{n}$ が成り立つ。 すなわち、$b^{\lambda-1}$ を掛けることで、$\to$ をひとつ戻ることができる¹。

また、頭の部分に $a$ が含まれるためには $\gcd(b^\mu, n)=\gcd(a, n)$ が必要であり、これが成り立たなければ解なしを報告すればよい。 よって、以下ではこれが成り立つとする。

$z = \mu+i\cdot\lfloor\sqrt{\lambda}\rfloor+j$ ($0\le j\lt\sqrt{\lambda}$) とし、$(i, j)$ を探すことを考える。 $$\begin{aligned} b^{\mu+i\cdot\lfloor\sqrt{\lambda}\rfloor+j} &\equiv a \\ b^{\mu+i\cdot\lfloor\sqrt{\lambda}\rfloor+j}\cdot (b^{\lambda-1})^j &\equiv a\cdot (b^{\lambda-1})^j \\ b^{\mu+i\cdot\lfloor\sqrt{\lambda}\rfloor} &\equiv a\cdot (b^{\lambda-1})^j \\ \end{aligned}$$ これより、$i$ と $j$ を分離でき、$a\cdot(b^{\lambda-1})^j\mapsto j$ の連想配列を作っておくことで、$b^{\mu+i\cdot\lfloor\sqrt{\lambda}\rfloor}$ に対応する要素があれば返せる。ただし、複数の $j$ について $a\cdot(b^{\lambda-1})^j$ が同じ値を取りうるので注意する必要がある。 $j$ による小さい幅で連想配列を作り、$i$ による大きい幅でそれにアクセスする様子から baby-step giant-step algorithm と呼ばれている。

About $\mu$ and $\lambda$

$\mu$ と $\lambda$ を求める方法について述べる。 自然数 $x$ が素数 $p$ で割り切れる回数、すなわち $x$ の素因数分解における $p$ の指数を $\pi_p(x)$ で表すとする²。 $b = \prod_{p_i:\text{ prime}} p_i^{e_i}$ とすると、次が成り立つ。 $$\mu = \max_i \left\lceil\frac{\pi_{p_i}(n)}{e_i}\right\rceil.$$ また、$p = \prod_i p_i^{\pi_{p_i}(n)}$, $n' = n / p$ とする。すなわち、$n$ を $b$ の各素因数で割り切れるだけ割ったものを $n'$ とする。このとき、$\lambda$ は $\varphi(n')$ の約数³であり、次の式に基づき愚直に求める⁴。 $$\lambda = \min_{i\sqsubseteq \varphi(n')} i\quad\text{s.t.}\quad b^\mu\cdot b^i \equiv b^\mu \pmod{n}.$$ $\varphi(n')$ の約数の個数は $o(\sqrt{n}/\log(n))$ なので、繰り返し二乗法を用いて $O(\sqrt{n})$ 時間でできる。 実際には、解が見つかった時点で打ち切れるので、$O(\lambda + \sum_{i\sqsubseteq \varphi(n')}^{i\le\lambda} \log(i))$ 時間になる。前半の $\lambda$ は試し割り法に由来する項なので、約数列挙を $O(1)$ delay でできるなら消去できる。

Proof

まず、$\lambda$ が $\varphi(n')$ の約数になることを示す。 $b^\mu$ における $p_i$ の指数は $e_i\cdot\max_j \lceil\pi_{p_j}(n)/e_j\rceil$ であり、特に $e_i\cdot \lceil\pi_{p_i}(n)/e_i\rceil \ge \pi_{p_i}(n)$ 以上である。 よって、$b^\mu$ は $p$ で割り切れ、$b^\mu \equiv 0\pmod{p}$ である。これより、$b^{\mu+i} \equiv 0\pmod{p}$ ($i\ge 0$) である。

また、$b$ と $n'$ は互いに素なので、Euler の定理より $b^{\varphi(n')} \equiv 1\pmod{n'}$ が成り立つ。よって、$b^{\mu+\varphi(n')}\equiv b^\mu \pmod{n'}$ である。上記の $p$ に関する結果と合わせて $b^{\mu+\varphi(n')}\equiv b^\mu \pmod{n}$ を得る。よって、$\lambda$ が $\varphi(n')$ の約数であることが示された。

次に、$\mu$ の値が正しいことを示す。$b^\mu$ が頭の部分に含まれることは上で示した $b^{\mu+\varphi(n')}\equiv b^\mu\pmod{n}$ より従うので、$\mu\gt 0$ のとき $b^{\mu-1}$ が頭の部分に含まれないことを示せば十分である。 $\mu$ の定義より、ある $i$ に対して $e_i\cdot(\mu-1)\lt \pi_{p_i}(n)$ が成り立つ。 よって $b^{\mu-1}$ は $p$ は割り切れず、$b^{\mu-1} \not\equiv 0\pmod{n}$ である。 これにより、頭の部分に含まれるどの要素とも等しくないことがわかる。

最後に、$\gcd(b^\mu, n)=\gcd(a, n)$ であることが、頭の部分に $a$ が存在する、すなわち ${}^\exists i\ge 0: b^{\mu+i} \equiv a\pmod{n}$ であることの必要条件であることを示す。

• $\gcd(b^{\mu+i})\equiv 0\pmod{p}$ for $i\ge 0$
• $\gcd(b, n') = 1$ より $\gcd(b^{\mu+i}, n') = 1$ for $i\ge 0$
• $\gcd(p, n') = 1$
• $n=p\cdot n'$

より、$\gcd(b^{\mu+i}, n) = p$ ($i\ge 0$) で一定となる。 すなわち、$\gcd(a, n) = p$ が必要となる。

Complexity

$O(\sqrt{n} + \sqrt{\lambda}\cdot H(\sqrt{\lambda}))$ time.

ただし、$H(n)$ は要素数 $n$ の HashMap の insert と get にかかる時間とする。 BS/GS パートでは BTreeMap を用いるよりも HashMap を用いた方が高速だったので、とりあえずそうした。

ここで、$\lambda$ は最悪 $\Theta(n)$ である。 また、篩の前計算などで $n$, $b$ の素因数分解と $\varphi(n)$ の約数列挙を $O(1)$ delay でできるなら $\sqrt{n}$ の項をそれらの個数に落とせる。

Implementation notes

$\lambda$ を求める際に $b^{\mu+i}\equiv b^\mu$ なる最小の $i\sqsubseteq \varphi(n')$ を探したが、実際には Euler の定理ではなく Carmichael の定理に基づいて $\lambda(n')$ の約数を探せばよい。関数値を求めるの自体は $\varphi(\bullet)$ よりも $\lambda(\bullet)$ の方が若干複雑だが、全体では高速に動作したので後者を選択した。

Examples

use nekolib::math::DLog;

assert_eq!(6_u64.dlog(5, 13), Some(9));
assert_eq!(27_u64.dlog(3, 30), Some(3));
assert_eq!(2_u64.dlog(0, 4), Some(2));
assert_eq!(0_u64.dlog(1, 2), Some(0));
assert_eq!(2_u64.dlog(3, 10), None);

References

• https://divinejk.hatenablog.com/entry/2021/02/07/133155

Suggestions

内部で Euler の $\varphi$ 関数 Carmichael の $\lambda$ 関数と約数列挙と素因数分解が必要となるので、 篩などを用いて高速化が図れる場合に対応しやすいようにしてみるとよい？

pub fn dlog<I, J>(
    b: u64, a: u64, n: u64,
    // euler_phi: impl Fn(u64) -> u64,
    carmichael_lambda: impl Fn(u64) -> u64,
    divisors: impl Fn(u64) -> I,
    factors: impl Fn(u64) -> J,
) -> Option<u64>
where
    I: Iterator<Item = u64>,
    J: Iterator<Item = (u64, u32)>,
{ ... }

一方で、そうでない場合にわざわざ関数を渡す必要があるので単に面倒かも。 お手軽パターンと二つ用意するといいかも？（デフォルト引数欲しいね）

とはいえ、結局 BS/GS パートが律速になってしまいそう。

Naming

関数名は discrete logarithm、引数 $b$ は base、$a$ は antilogarithm から。

---

1. $\mu\gt 0$ なら、 $b^\mu$ は $b^{\mu-1}$ と $b^{\mu+\lambda-1}$ ($\neq b^{\mu-1}$) のふたつから $\to$ が入ってくるが、$\rho$ の頭に含まれる後者に戻る。 ↩

2. 標準的な記法があったら知りたいです。 ↩

3. $\varphi(n)$ と $n'$ の定義から、$\varphi(n')$ は $\varphi(n)$ の約数であることに注意。 ↩

4. $x\sqsubseteq y$ で $x$ が $y$ を割り切ることを表す。$x\mid y$ は方向がわかりにくくて嫌い。素因数の多重集合の包含関係のつもり。$\subseteq$ でもいいかも？ ↩

Required Methods

fn dlog(self, a: Self, n: Self) -> Option<Self>

Implementations on Foreign Types

impl DLog for u32

fn dlog(self, a: Self, n: Self) -> Option<Self>

impl DLog for u8

fn dlog(self, a: Self, n: Self) -> Option<Self>

impl DLog for usize

fn dlog(self, a: Self, n: Self) -> Option<Self>

impl DLog for u128

fn dlog(self, a: Self, n: Self) -> Option<Self>

impl DLog for u64

fn dlog(self, a: Self, n: Self) -> Option<Self>

impl DLog for u16

fn dlog(self, a: Self, n: Self) -> Option<Self>

Implementors