🔖 math数论扩展欧几里得算法中国剩余定理Baby Step Gaint Step

扩展欧几里得算法

问题描述

已知 $a,b$ 为常整数，求方程 $ax+by=\gcd(a,b)$ 的一对整数解?

问题简析

令 $\displaystyle \left\lbrace\begin{aligned}&A=a \mod b\\&B=b\end{aligned}\right.$，构造方程

$$\begin{align} AX+BY=\gcd(A,B). \tag{1}\end{align}$$

由于 $\gcd(a,b)=\gcd\big(a \mod b,b\big)=\gcd(A,B)$，所以

$$\begin{align} ax+by &=\left(a-\left\lfloor \frac{a}{b} \right\rfloor \times b+\left\lfloor \frac{a}{b} \right\rfloor \times b\right)\times x+by \notag \\ &=\left(a-\left\lfloor \frac{a}{b} \right\rfloor \times b\right)\times x+\left(by +\left\lfloor \frac{a}{b} \right\rfloor \times bx\right) \notag \\ &=\big(a (\hskip -1em \mod b)\big)\times x+b\times\left(y+\left\lfloor \frac{a}{b} \right\rfloor \times x\right) \notag \\ &=Ax+B\times\left(y+\left\lfloor \frac{a}{b} \right\rfloor \times x\right) \notag \end{align}$$

因此有 $\displaystyle \left\lbrace\begin{aligned} &X=x \\ &Y=y+\left\lfloor \frac{a}{b} \right\rfloor \times x \end{aligned} \right. \quad\Rightarrow\quad \left\lbrace\begin{aligned} &x=X \\ &y=Y- \left\lfloor \frac{a}{b} \right\rfloor \times X \end{aligned}\right.$

方程 (1) 的构造与辗转相除有异曲同工之妙，需要说明的是， $\gcd(x,0)=x$。所以：

  1. 当 $b=0$ 时，我们得到一组解 $\displaystyle \left\lbrace\begin{aligned} &x=1 \\ &y=0 \end{aligned}\right.$

  2. 当 $b\neq 0$ 时，套用辗转相除的框架求出 $(X,Y)$，再得到 $(x,y)$。

程序实现

1
2
3
4
5
6
7
8
9
10
11
typedef long long LL;
void extendEuclid(LL a, LL& x, LL b, LL &y, LL& d) { // d=gcd(a,b)
  if (!b) { 
    d = a;
    x = 1;
    y = 0; 
  } else { 
    extendEuclid(b, y, a%b, x, d);
    y -= (a / b) * x; 
  }
}

上述代码中，将 d = a 的执行顺序放在最前面是一个小技巧，因为这样就可以像下面这样调用函数，此时不需要预先将 $gcd(a,b)$ 求出来，而是在整个算法中顺带求出的：

1
2
3
extendEuclid(a, x, b, d, d);
// or
extendEuclid(a, d, b, d, y);

值得一提的是，扩展欧几里得算法解出的一组解是所有解中$|x|+|y|$ 最小的；

问题扩展

求方程 $ax+by+c=0$ 有多少对整数解，满足 $x \in [x_1,x_2],~y \in [y_1,y_2]$？

由前文可知，当 $\gcd(a,b) \mid c$ 时，方程有解。因为方程 $ax+by=\gcd(a,b)$ 的解乘上 $\frac{c}{\gcd(a,b)}$ 就是此方程的一组解。由此，我们可以推出 $ax+by=k{\bf\cdot}\gcd(a,b)$，所以，当 $\displaystyle \gcd(a,b) \nmid c$ 时，方程无解。

那有多少个解呢？

如果 $(x_0,y_0),\;(x_1,y_1)$ 是方程的两个解，那么 $a(x_1-x_0)+b(y_1-y_0) \equiv 0$.
令 $\displaystyle a'=\frac{a}{\gcd(a,b)},~b'=\frac{b}{\gcd(a,b)}$，显然，$\gcd(a',b')=1$
所以，$(x_1-x_0) \mid b',~(y_0-y_1) \mid a'$，并且 $\displaystyle \frac{x_1-x_0}{b'} \equiv \frac{y_0-y_1}{a'}$。

一般地，若 $(x_0,y_0)$ 是方程 $ax+by+c=0$ 一组解，那么 $(x_0+kb',y_0-ka')$ 也是一组解。其中，$\displaystyle b'=\frac{b}{\gcd(a,b)},~a'=\frac{a}{\gcd(a,b)}$， $k$ 为任意整数。

逆元

对于任意整数 $A,~N$，如果存在一个整数 $B$ 使得 $A{\bf\cdot}B \equiv 1 \mod C$，那么，$A$ 和 $B$ 称作模 $C$ 意义下的互为逆元，记做 $A^{-1}=B$ 或 $B^{-1}=A$。

求解逆元

  1. 如果 $\gcd(A,N)=1$，根据欧拉定理，有 $A^{-1} = N^{\varphi(N)-1}$，特别地，当 $N$ 是素数时，$\varphi(N)=N-1$。另外，由于 $\gcd(A,N)=1$，构造方程 $Ax+Ny=1$，用扩展欧几里得算法也可求出 $A$ 在模 $N$ 意义下的逆元 $A^{-1}=x$.

  2. 如果 $\gcd(A,N)\neq 1$，由前文可知，方程 $Ax+Ny=1 \neq k\cdot\gcd(A,N)$ 不存在整数解，故不存在逆元。

程序实现

  1. 欧拉定理求逆元。当 $N$ 是素数时，$\varphi(N)=N-1$，$A$ 的逆元为 $A^{N-2}$

    inv.euler.cpp 
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    typedef long long LL;
    LL powerMod(LL A, LL x, LL mod) { // 返回 A^x%mod
      LL ans = 1;
      for(; x > 0; x >>= 1, A = A * A % mod) {
        if (x & 1) ans = ans * A % mod;
      }
      return ans;
    }
    

    LL inv(LL A, LL mod) {
      return powerMod(A, mod-2, mod);   
    }

  2. 扩展欧几里得求逆元

    inv.euclid.cpp 
    1
    2
    3
    4
    5
    6
    7
    typedef long long LL;
    LL inv(LL A, LL mod) {
      LL x, y, d;
      extendEuclid(A, x, N, y, d);
      // 若 A 和 mod 不互质，则不存在逆元
      return d == 1 ? x : -1;
    }

中国剩余定理

若正整数 $m_1,m_2,\cdots,m_k$ 两两互素，且 $0 \leqslant a_1 < m_1, 0 \leqslant a_2 < m_2,\cdots, 0 \leqslant a_k < m_k$. 求下列同余方程组的整数解：

$$\begin{align} \begin{aligned} &x \equiv a_1(\hskip -1em \mod m_1) \\ &x \equiv a_2(\hskip -1em \mod m_2) \\ &\hskip 2.5em \cdots \\ &x \equiv a_k(\hskip -1em \mod m_k) \end{aligned} \tag{2} \end{align}$$

问题简析

令 $\displaystyle w_i=\prod_{j \neq i} m_j$，显然有：

$$\begin{align} \gcd(w_i,m_j) = \left\lbrace \begin{aligned} &1, &i = j \\ &m_j, &i \neq j \end{aligned} \right. \tag{3} \end{align}$$

构造方程 $w_i x_i+m_i y_i=\gcd(w_i,m_i)=1$。那么，由方程组 (3) 有，

$$\begin{align} w_i x_i ~(\hskip -1em \mod m_j) \equiv \left\lbrace \begin{aligned} &1, &i=j \\ &0, &i\neq j \end{aligned} \right. \tag{4} \end{align}$$

结合方程组 (2) 和方程组 (4)，进而有，

$$\begin{align} a_i w_i x_i ~(\hskip -1em \mod m_j) \equiv \left\lbrace \begin{aligned} &a_i, &i=j \\ &0, &i\neq j \end{aligned} \right. \tag{5} \end{align}$$

所以，$\displaystyle x=\sum_{i=1}^k (a_i\cdot w_i\cdot x_i)$ 为同余方程组 (2) 的解。更一般地，方程的解为

$$x=\sum_{i=1}^k (a_i\cdot w_i\cdot x_i) + k\prod_{i=1}^k m_i$$

剩下的问题仅需利用扩展欧几里得算法求出 $x_i$。

程序实现

crt.cpp 
1
2
3
4
5
6
7
8
9
10
LL crt(int N, int* a, int* m) {
  LL M = 1, ans = 0;
  for (int i = 0; i < N; ++i) M *= m[i];
  for (int i = 0; i < N; ++i) {
    LL w = M / m[i], x, y, d;
    extendEuclid(w, x, m[i], y, d);
    ans = (ans + x * w * a[i]) % M;
  }
  return ans;
}

Baby Step Gaint Step

求一个整数 $x \in [0,C-1]$ 使得 $A^x \equiv B(\hskip -.6em \mod C)$。其中，$C$ 为素数。

问题简析

注意到 $C$ 是一个素数，那么对于任意 $X$ 都存在一个模 $C$ 意义下的逆元 $X^{-1}$，使得 $X \cdot X^{-1} \equiv 1 \mod C$。取任意正整数 $m$，如果事先知道所有的 $A^{x_0} \mod C,\;(0 \leqslant x_0 < m)$ 的值，如何判断是否存在一个 $x'$ 满足 $km \leqslant x'=x_0+km < (k+1)m$，使得 $A^{x'} = A^{x_0+km} \equiv B \mod C$ 成立呢？

注意到，令方程两边同乘以 $A^{km}$ 的逆元 $A^{-km}$，得到 $A^{x_0} \equiv B\cdot A^{-km} \mod C$。所以，我们仅需将映射关系 $x_0 \rightarrow f(x_0)=A^{x_0} \mod C,\;(0 \leqslant x_0 < m)$ 存进哈希表或平衡树中，再从 $0 \sim m$ 枚举 $k$ 的值，对于每个 $k$ 检查哈希表或平衡树中是否存在 $B\cdot A^{-km} \mod C$，若存在即可取得对应的 $x_0$ 的值。

$m$ 应该取多少合适呢？注意到这样做的时间复杂度是 $\displaystyle O\left( m+\frac{C}{m} \right)$ 的 [1]，不难证明，当 $m$ 取 $\sqrt{C}$ 时复杂度最优。

程序实现

bsgs.cpp 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
typedef long long LL;
unordered_map<LL, int> emp;


int bsgs(LL A, LL B, LL C) {
  emp.clear();
  int M = ceil(sqrt(1.0 * C));


  LL AA = 1;
  for (int x0 = 0; x0 < M; ++x0) {
    if (!emp.count(AA)) emp[AA] = x0;
    AA = AA * A % C;
  }


  LL Am = inv(powerMod(A, M, C), C);
  for (int k = 0; k < M; ++k) {
    if (emp.count(B)) return emp[B] + k * M;
    B = B * Am % C;
  }
  return -1;
}

Extend Baby Step Gaint Step

求一个整数 $x \in [0,C-1]$ 使得 $A^x \equiv B(\hskip -.6em \mod C)$。其中，$C$ 为合数。

问题简析

在分析之前，先看一个定理：若 $k\cdot A \equiv k\cdot B\mod k\cdot C$，则有，$A \equiv B \mod C$.
这个定理不难证明，事实上，仅需构造方程 $k\cdot A - k\cdot B=k'\cdot k\cdot C$，在等式两边同除以 $k$ 即可得证。

再回过头来看普通版的 Baby Step Gaint Step 是如何工作的。注意到， Gaint Step 能够成功走出去的前提是 $A^m$ 存在模 $C$ 意义下的逆元 $A^{-m}$。由前文可知，$A^{-m}$ 存在的充要条件为 $\gcd(A^m,C)=1$，此时必有 $\gcd(A,C)=1$ [2]。设 $\displaystyle G_t =\gcd(A^t,C), \; C'=\frac{C}{G_t}, \; D=\frac{A^t}{G_t}$，根据上面的定理，不难想到，我们只要找出某个整数 $0 \leqslant t < C$ 满足 $G_t=G_{t+1}$，则上述方程可转化为

$$A^x = G_t \times \frac{A^t}{G_t} \times A^{x-t} = G_t\times D\times A^{x-t}\\C=G_t\times \frac{C}{G_t}=G_t\times C'$$

不难发现，$D$ 和 $C'$ 是互质的；因为 $G_t=G_{t+1}$，所以 $A$ 和 $C'$ 也是互质的 [3] 。

接下来讨论 $B$ 和 $G_t$ 的关系：

  • 若 $G_t \mid B$，结合前面给出的定理，有

    $$A^x \equiv B \mod C \quad\Rightarrow\quad D\cdot A^{x-t} \equiv \frac{B}{G_t} \mod C'$$

    又因为 $D$ 和 $C'$ 是互质的，所以存在模 $C'$ 意义下的逆元 $D^{-1}$；进一步将同余方程化为

    $$\displaystyle A^{x-t} \equiv \left(D^{-1}\cdot \frac{B}{G_t}\right) \mod C'$$

    由于 $A$ 和 $C'$ 是互质的，根据前面的分析，直接套用普通版的 Baby Step Gaint Step 即可解决，这部分算出的答案加上 $t$ 就是原同余方程的解。

  • 若 $G_t \nmid B$，由于 $A^x \mod C = (G_t\times D\times A^{x-t}) \mod (G_t\times C') = K\cdot G_t$，所以无解。

程序实现

ebsgs.cpp 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
typedef long long LL;
unordered_map<LL, int> emp;


int ebsgs(LL A, LL B, LL C) {
  LL D = 1, cnt = 0;
  for (LL G; (G = gcd(A, C)) != 1; cnt += 1) {
    if (B % G) return -1;
    B /= G;
    C /= G;
    D = D * A / G % C;
    if (D == B) return cnt;
  }


  emp.clear();
  int M = ceil(sqrt(1.0 * C));


  LL AA = 1;
  for (int x0 = 0; x0 < M; ++x0) {
    if (!emp.count(AA)) emp[AA] = x0;
    AA = AA * A % C;
  }


  B = B * inv(D, C) % C;
  LL Am = inv(powerMod(A, M, C), C);
  for (int k = 0; k < M; ++k) {
    if (emp.count(B)) return emp[B] + cnt + k * M;
    B = B * Am % C;
  }
  return -1;
}

练习

#ProblemCategorySolution / Code
151Nod/X^A Mod PBSGSsolution.cpp
251Nod/X^A Mod B扩展BSGSsolution.cpp
3SPJ/MOD扩展 BSGSsolution.cpp
4UVa/Code Feat中国剩余定理solution.cpp

小结

拖沓了很久，终于耐下性子整理了一下。
鉴于本人能力有限，有误之处欢迎指正。

Related

footnote-definitions

  • ↑ [1]: 

    使用哈希表，查询均摊 $O(1)$

  • ↑ [2]: 

    若 $\gcd(A,C) \neq 1$，则 $gcd(A^m,C) \geqslant \gcd(A,C) > 1$

  • ↑ [3]: 

    $G_t=G_{t+1}$ 说明 $\displaystyle \gcd \left(\frac{C}{G_t}, A \right)=1$

© 2017-2025 光和尘有花满渚、有酒盈瓯

Comments