【瞎口胡】基础数学 2 裴蜀定理 exgcd 同余相关

写在前面

在第一篇当中我们介绍了一点基本知识。

第二篇的内容较第一篇稍有难度，但是也非常基础，相信只要学过小学数学就能看明白。

裴蜀定理与扩展欧几里得算法

• 裴蜀定理

  对于任意非负整数 (a,b)：

  1. 任意整数 (x,y) 都满足 (ax+by) 是 (gcd(a,b)) 的倍数。

     证明：(gcd(a,b) mid a,gcd(a,b) mid b)，则一定有 (gcd(a,b) mid ax +by)。

  2. 存在整数 (x,y) 使得 (ax+by = gcd(a,b))。

     证明：数学归纳法。回忆欧几里得算法的最后一步，(a_0 = gcd (a,b),b_0 = 0)，取 (x'=1,y'=0)，即有 (a_0 x' + b_0y' = gcd(a_0,b_0))。

     根据欧几里得算法的流程，(a_0 = b,b_0 = amod b)。

     [a_0 x' + b_0y' = gcd(a_0,b_0) ]
     [bx' + (a mod b)y' = gcd(a,b) ]
     [bx' + (a- leftlfloor dfrac{a}{b} ight floor b)y'=gcd(a,b) ]
     [ay'+b(x'- leftlfloordfrac{a}{b} ight floor y')=gcd(a,b) ]

     因此有 (x=y',y=x'-leftlfloordfrac{a}{b} ight floor y')。对于每一层的 (a,b) 都是如此，在回溯完成之后就会得到一组正确的解。

• 扩展欧几里得算法

  简称扩欧，exgcd。

  按照证明裴蜀定理的方式，将欧几里得算法稍作修改：

  inline int exgcd(int _a,int _b,int &_x,int &_y){
  	if(!_b){
  		_x=1,_y=0; // 最后一步
  		return _a;
  	}
  	int g=exgcd(_b,_a%_b,_x,_y),Temp; // 先算出下一层的答案
  	Temp=_x,_x=_y,_y=Temp-(_a/_b)*_y; // 计算这一层的答案 注意顺序是先算 x 再算 y，还需要一个辅助变量
  	return g;
  }
  

• 应用

  • 求方程 (ax+by=gcd(a,b)) 的解集

    考虑我们 exgcd 求出了该方程的一组特解 (x_0,y_0)。我们希望在此基础上将 (x_0,y_0) 加上一个值，设其为 (Delta x,Delta y)。

    那么它们必须满足 (aDelta x + b Delta y = 0)。即

    [Delta y = -dfrac{a}{b}Delta{x} ]

    将 (dfrac ab) 约分，设其最简形式为 (dfrac{a'}{b'})，显然 (a' = dfrac{a}{gcd(a,b)},b' = dfrac{b}{gcd(a,b)})。此时有

    [Delta y = -dfrac{a'}{b'}Delta{x} ]

    因为 (Delta y) 是整数，因此 (Delta{x}) 必须是 (b') 的倍数。设 (Delta x = l cdot b')，则 (Delta y = - l cdot a')，其中 (l) 是任意整数。

    综上，(ax + by = gcd(a,b)) 的解集为 (x = x_0 + k dfrac{b}{gcd(a,b)}, y = y_0 - k dfrac{a}{gcd(a,b)})。

  • 求方程 (ax + by = c) 的解集

    由裴蜀定理，(ax + by) 一定是 (gcd(a,b)) 的倍数。因此如果 (gcd(a,b) mid c)，该方程无解。

    令 (g = gcd(a,b)，c = t imes g)。

    [ax + by = c ]
    [ax + by = g cdot t ]
    [adfrac xt + b dfrac yt = g ]

    此时用扩欧求出该方程的解 (x'_0,y'_0)，则原方程的一组特解为 (x_0 = t x'_0 ,y_0= t y'_0)。

    利用求 (ax + by = gcd(a,b)) 时的证明方法，可以知道该方程的解集为 (x = x_0 + k dfrac{b}{gcd(a,b)}, y = y_0 - k dfrac{a}{gcd(a,b)})。

同余

对于整数 (a,b) 和正整数 (p)，若 (a mod p = b mod p)，则成 (a equiv b pmod p)，读作「(a) 在模 (p) 意义下与 (b) 同余」。

基本性质

• 性质 (1)：自反性

  [a equiv a pmod p ]

• 性质 (2)：对称性

  [a equiv b Rightarrow b equiv a pmod p ]

• 性质 (3)：传递性

  [a equiv b,b equiv c Rightarrow a equiv c pmod p ]

• 性质 (4)：可加、可乘性

  [a equiv b Rightarrow an equiv bn ,a pm n equiv b pm n pmod p ]

• 性质 (5)

  当 (c) 与 (p) 互质的时候：

  [ac equiv bc Rightarrow a equiv b pmod p ]

  证明：考虑 (a equiv b) 的充要条件是 (p mid (a-b))。

  (ac equiv bc) 可以推出 (p mid (ac-bc))，即 (p mid (a-b)c)。当 (c) 与 (p) 互质的时候，(a-b) 才一定被 (p) 整除。换句话说，(gcd(c,p) = 1) 是上式成立的充分条件。

  性质 (5) 非常重要，一定要牢记只有 (gcd(c,p) = 1) 的时候该性质才成立。

剩余系与剩余类

所有对 (p) 同余的整数构成模 (p) 的一个剩余类。考虑一个整数 (mod p) 的取值范围是 ([0,p-1])，于是剩余类显然有 (p) 个。若某个剩余类中的整数模 (p) 的结果与 (p) 互质，则称该剩余类为互质剩余类。

注意，模 (p) 为 (0) 的整数构成的剩余类不是互素剩余类。

完全剩余系：在所有同余类中各任取一个整数，它们构成了模 (p) 的一个完全剩余系。

简化剩余系（缩系）：在所有互质剩余类中任取一个整数，它们构成了模 (p) 的一个缩系。

完全剩余系和缩系的个数是无限的。

费马小定理

对于质数 (p) 和任意满足 (gcd(a,p) = 1) 的 (a)，有

[a^{p-1} equiv 1 pmod p ]

• 引理

  对于满足上述条件的 (a)，(a,2a,3a,...,(p-1)a) 构成模 (p) 的一个缩系。

  证明：显然上面的 (p-1) 个数与 (p) 互质。如果存在 (0 leq i eq j < p) 使得 (ai equiv aj pmod p)，由性质 (5) 将两边同时除以 (a)，得 (i equiv j pmod p)，因为 (0 < i,j < p)，即 (i = j)，与假设矛盾。

显然

[1 imes 2 imes ... imes (p-1) equiv a imes 2a... imes (p-1) a pmod p ]

[1 imes 2 imes.. imes (p-1) equiv 1 imes 2 imes ... imes (p-1) imes a^{p-1} pmod p ]

由性质 (5)，将两边化简并调换位置得

[a^{p-1} equiv 1 pmod p ]

欧拉定理

对于任意模数 (p) 和任意满足 (gcd(a,p) = 1) 的 (a)，有

[a^{varphi(p)} equiv 1 pmod p ]

其中 (varphi(x)) 表示小于 (x) 的正整数中与 (x) 互质的个数。

证明和证明费马小定理类似，不再赘述。

注意到费马小定理是欧拉定理的特例。

乘法逆元

对于整数 (a) 和模数 (p)，称 (ax equiv 1 pmod p) 的整数 (x) 为 (a) 在模 (p) 意义下的乘法逆元。(a) 的乘法逆元记为 (a^{-1})。

求法 (1)：对于 (p) 是质数的情况下，由费马小定理有 (a imes a^{p-2} equiv 1 pmod p)，则 (a) 的乘法逆元为 (a^{p-2})。

求法 (2)：对于 (p) 不是质数的情况下，有 (ax - py = 1)，此时 (a,p) 为定值，使用 exgcd 求出 (x) 的最小非负整数解即可。

上述求法的时间复杂度均为 (O(log n))，由求法 (2) 可知，(a) 在模 (p) 意义下存在逆元的充要条件是 (gcd(a,p)=1)。

求法 (3)：在已知 (forall 1 leq i leq n,gcd(i,p)=1) 的情况下，怎样快速求出 (1 sim n) 的每一个数在模 (p) 意义下的逆元？

我们有一个线性的做法。对于 (i)，如果我们知道 (1sim i-1) 中每一个数的逆元，那么将 (p) 写成 (ki + r(0 leq r<i)) 的形式，则有：

[ki + r equiv 0 pmod p ]

左右两侧同乘 (i^{-1} imes r^{-1})：

[k imes r^{-1} + i^{-1} equiv 0 pmod p ]

整理得

[i^{-1} = (-left lfloor dfrac{p}{i} ight floor imes (p mod i)^{-1}) mod p ]

。

	long long inv[MAXN];
	......
	n=read(),m=read(); // m 为模数
	inv[1]=1;
	for(rr int i=2;i<=n;++i){
		inv[i]=(m-m/i)*inv[m%i]%m; // m-m/i 为模意义下的 -p/i
	}
	for(rr int i=1;i<=n;++i){
		printf("%lld
",inv[i]);
	}

记忆技巧：用 (i imes dfrac pi + p mod i) 表示 (p)，再乘上逆元。

阶乘逆元

如果题目的模数 (p) 较大且为质数，而 (n<p)，则 (n!) 一定和 (p) 互质，于是根据费马小定理，(n!) 一定存在模 (p) 意义下的逆元。

方法 (1)：我们可以使用 (n) 次费马小定理，利用快速幂进行计算，从而得到每一个 (i!(1 leq i leq n))，时间复杂度 (O(n log n))。

方法 (2)：观察到 ((i+1)!=i! imes (i+1))。我们可以用 (O(n)) 的时间计算出 ([1,n]) 中每一个数的逆元和 (n!)，然后倒着递推。

线性同余方程组

中国剩余定理 / CRT

给定同余方程组

[left{egin{aligned}xequiv a_1pmod{{m_1}}\ xequiv a_2pmod{{m_2}} \ xequiv a_3 pmod{{m_3}} \ cdots \ xequiv a_npmod{{m_n}}end{aligned} ight. ]

其中 (a_i) 是任意整数，(m_i) 两两互质。

不失一般性，令 (0 leq a_i < m_i)。设 (p = prod limits_{i=1}^{n} m_i)，则该方程存在一解 (x_0 = sum limits_{i=1}^{n} a_it_iM_i)，其中 (M_i = dfrac{p}{m_i})（即：(M_i = prod limits _{1 leq r leq n and r eq i} m_r)），(t_i) 是 (M_i) 在模 (m_i) 意义下的乘法逆元。

证明：对于每一个 (i(1 leq i leq n)，)考虑和式中的任意 (a_jt_jM_j (j eq i))。由 (M_j) 的定义知它是 (m_i) 的倍数，于是 (a_jt_jM_j equiv 0 pmod{m_i})。再来考虑和式中 (a_it_iM_i) 一项。由乘法逆元和 (t_i) 的定义只 (t_iM_i = 1)，于是 (a_it_iM_i equiv a_i pmod {m_i})。综上，(x_0 equiv a_i pmod{m_i})，满足该约束条件。

中国剩余定理断言该方程的通解是 (x= x_0 + kp(k in mathbb Z))。易证这一定是原方程组的解。接下来证明不存在其它任意解。

考虑有同余方程组

[left{egin{aligned}xequiv a_1pmod{{m_1}}\xequiv a_2 pmod{m_2}end{aligned} ight. ]

那么对于任意满足 (m_1j_1+a_1=m_2j_2+a_2) 的整数 (j_1,j_2)，方程的解 (x=m_1j_1+a_1)。

整理得 (m_1j_1-m_2j_2=a_2-a_1)。令 (j_2'=-j_2)，即 (m_1j_1+m_2j_2'=a_2-a_1)。容易观察到这是 exgcd 的形式。观察到如果 (gcd(m_1,m_2)mid a_2-a_1)，那么一定存在一组 (j_1,j_2') 的特解 (j_{1,0},j_{2,0}')。同时，我们也知道了，该方程有解当且仅当 (gcd(m_1,m_2) mid a_2 - a_1)（不失一般性，我们假设 (a_2 geq a_1)）。由裴蜀定理，(j_1 = j_{1,0} + kdfrac{m_2}{gcd(m_1,m_2)},j_2'=j_{2,0}'-kdfrac{m_1}{gcd(m_1,m_2)}(k in mathbb Z))。

则

[x = m_1j_1+a_1 ]

[= m_1(j_{1,0} + kdfrac{m_2}{gcd(m_1,m_2)} +a_1 ]

[=m_1j_{1,0}+km_1m_2gcd(m_1,m_2) + a_1 ]

[=koperatorname{lcm}(m_1,m_2) + m_1j_{1,0} + a_1 ]

于是我们将之前的同余方程组变成了一个同余方程

[x equiv m_1j_{1,0} + a_1 pmod{operatorname{lcm}(m_1,m_2)} ]

这证明了我们可以将这两个同余方程合并成一个模数为 (operatorname{lcm}(m_1,m_2)) 的新同余方程。同时我们得到了一个结论：不存在该方程组的两个解 (x_1 = k_1operatorname{lcm}(m_1,m_2) + r_1,x_2 =k_2operatorname{lcm}(m_1,m_2) +r_2) 使得 (r_1 eq r_2)，因为 (x) 的所有解在模 (operatorname{lcm}(m_1,m_2)) 意义下同余。

将该结论推广到原同余方程组

[left{egin{aligned}xequiv a_1pmod{{m_1}}\ xequiv a_2pmod{{m_2}} \ xequiv a_3 pmod{{m_3}} \ cdots \ xequiv a_npmod{{m_n}}end{aligned} ight. ]

我们可以说，不存在该方程组的两个解 (x_1 = k_1operatorname{lcm}(m_1,m_2,cdots,m_n) + r_1,x_2 =k_2operatorname{lcm}(m_1,m_2,cdots,m_n) r_2) 使得 (r_1 eq r_2)。

由于 (m_i) 两两互质，即不存在该方程组的两个解 (x_1 = k_1p + r_1,x_2 =k_2p+r_2) 使得 (r_1 eq r_2)。不失一般性，设 (0 leq x_0 < p)。假设存在解 (k'p+r'(0 leq r'< p)) 但 (r' eq x_0)，则我们可以找到另外一个解 (k'p+x_0)。(r' eq x_0)，于是这和我们得到的结论矛盾。

扩展中国剩余定理 / exCRT

问题与 CRT 的应用场景基本一致，但 (m_i) 不一定两两互质。中国剩余定理不再适用。

考虑求出前 (k-1) 个方程的一个解 (x)。记 ( ext {lcm} (m_1,m_2,...,m_{k-1})=M)。我们要做的是找到一个正整数 (t)，使得 (x + tM equiv a_i (mod m_i))。这很容易理解，因为前 (k-1) 个方程合并之后的模数就是 (operatorname{lcm}(m_1,m_2,cdots,m_{k-1}))。

移项，得 (tM equiv a_i - x (mod m_i))，可以 exgcd 求出，合并前 (k) 个方程即可。

该方程的通解是 (x + k imes ext{lcm}(m_1,m_2,...,m_n) (k in mathbb Z))。这很容易证明，像证明 CRT 那样就可以了。

注意到当 (m_i) 不一定两两互质时，可能会出现无解的情况。

高次同余方程 / BSGS 算法

BSGS 算法，全程 Baby Step, Giant Step 算法，一译「小步大步算法」。

其作用是求解形如 (a^xequiv b pmod m) 的方程，其中 (gcd(a,m) =1) 且 (m) 是质数。

由费马小定理，(a^{m-1}=1 pmod m)，于是我们只需要考虑在 ([0,m-1)) 范围的解。

BSGS 算法的名字「小步大步」解释了该算法的流程：

考虑将某个指数拆成 (ik-r) 的形式，其中 (0 leq r < k,k>1)。如果有 (a^{ik-r} equiv b pmod m)，那么也有 (a^{ik} equiv a^r imes b pmod m)。

我们可以用 (O(k)) 的时间计算出 (a^r mod m (0 leq r < k))（并存放在 hash 表或 map 中。然后，我们用 (O(dfrac{m}{k})) 的时间检查：对于每一个满足 $0 leq i $，hash 表或 map 中是否存在 (a^r imes b mod m) 和 (a^{ik} mod m) 相等。

总时间复杂度为 (O(k+dfrac{m}{k}))，当 (k= sqrt m) 时最优，为 (O(sqrt m))。如果使用了 map，则还需要带一个 (log) 的复杂度。

模板题 POJ2417

# include <cstdio>
# include <cmath>
# include <cstring>
# define int long long

const int MOD=9997;
const int N=100010,INF=0x3f3f3f3f;
int p,b,n;

struct HashTable{
	struct Hash_Edge{
		int key,value,next;
	}edge[100010];
	int head[10010],sum;
	inline void add(int key,int value){
		edge[++sum].key=key;
		edge[sum].value=value;
		edge[sum].next=head[key%MOD];
		head[key%MOD]=sum;
		return;
	}
	inline int query(int key){
		for(int j=head[key%MOD];j;j=edge[j].next){
			if(edge[j].key==key){
				return edge[j].value;
			}
		}
		return -1;
	}
}HashT;

inline int read(void){
	int res,f=1;
	char c;
	while((c=getchar())<'0'||c>'9')
		if(c=='-')f=-1;
	res=c-48;
	while((c=getchar())>='0'&&c<='9')
		res=res*10+c-48;
	return res*f;
}
inline int qpow(int d,int p,int mod){
	int ans=1;
	while(p){
		if(p&1){
			ans=ans*d%mod;
		}
		p>>=1,d=d*d%mod;
	}
	return ans;
}
# undef int
int main(void){
# define int long long
	while(~scanf("%lld%lld%lld",&p,&b,&n)){
		memset(HashT.head,0,sizeof(HashT.head));
		HashT.sum=0; 
		if(n==1){ //特判 x^0 = 1 (x != 0)
			printf("0
");
			continue;
		}
		int m=ceil(sqrt(double(p)));
		int now=1;
		for(int i=0;i<m;++i){ // baby step
			if(!i){
				HashT.add(now*n%p,i);
				continue;
			}
			now=now*b%p;
			HashT.add(now*n%p,i);
		}
		int base=qpow(b,m,p),ans=1;
		for(int i=1;i<=m;++i){ // giant step
			ans=ans*base%p;
			if(HashT.query(ans)!=-1){
				printf("%lld
",i*m-HashT.query(ans));
				goto END;
			}
		}
		puts("no solution");
		END:;		
	}

	return 0;
}