学习笔记——同余数论

前置芝士#

同余概念，逆元，费马小定理，Exgcd，欧拉函数及定理。

（扩展）中国剩余定理#

个人喜欢使用非公式形式的合并方式来求解同余方程组。

我们需要求解的问题是求合法的解 $x$，使之满足：

$$\begin{cases} x\equiv a_1\pmod {p_1}\\ x\equiv a_2\pmod {p_2}\\ \vdots\\ x\equiv a_n\pmod {p_n} \end{cases}$$

此时我们考虑一个人类智慧来合并两个方程组，如果我们能做到这一点，那么我们最终能把方程组变成一个方程，从而得到通解。所以现在我们对：

$$\begin{cases} x\equiv a_1\pmod {p_1}\\ x\equiv a_2\pmod {p_2} \end{cases}$$

进行合并。对 $x\equiv a_1\pmod{p_1}$，我们容易改写成：$x=a_1+k_1p_1$。那换句话说，我们希望找到一个 $k_1$，使得 $a_1+k_1p_1\equiv a_2\pmod{p_2}$，那我们再次故技重施，得到：$a_1+k_1p_1=a_2-k_2p_2$，移项就得到了：$k_1p_1+k_2p_2=a_2-a_1$，这玩意儿直接扩欧解一下就行了。

在此之前，我们实际上是构造了一个 $x\equiv 0\pmod 1$ 的方程作为边界来合并的。

可喜可贺的是这做法不需要保证 $\gcd(p_1,p_2,\cdots,p_n)=1$。

• P4777 【模板】扩展中国剩余定理（EXCRT）

#define int long long
using namespace std;
const int MAXN=1e5;
int exgcd(int a,int b,int &x,int &y){
	if(b==0){x=1,y=0;return a;}
	int g=exgcd(b,a%b,y,x);
	y-=a/b*x;return g;
}
int a[MAXN],p[MAXN];
signed main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	int n;cin>>n;
	rep(i,1,n) cin>>p[i]>>a[i];
	int ans=0,lcm=1,x,y,g;
	rep(i,1,n){
		int A=lcm,B=p[i],C=((a[i]-ans)%p[i]+p[i])%p[i];
		g=exgcd(A,B,x,y);x=(x%B+B)%B;
		ans+=(__int128)(C/g)*x%(B/g)*lcm%(lcm*=p[i]/g);
		ans%=lcm;
	}if(ans<0) ans+=((-ans-1)/lcm+1)*lcm;
	cout<<ans<<'\n';
	return 0;
}

• P4774 [NOI2018] 屠龙勇士

这题在求出了每次选用的剑之后，和例题的差别就是 $x$ 多了系数，其实本质是相同的，重新推一遍就行了：

我们肯定是想着把系数消掉的，所以我们假装已经消掉了此前的系数，即 $x\equiv a\pmod p$（可以看成是系数为 $1$，你可以发现上面我们构造的边界条件正是这样的）。假如说，它要与方程 $b_1x\equiv a_1\pmod{p_1}$ 合并，那么我们故技重施，得到 $b_1(a+k_1p)=a_1-k_2p_1$，移项得到 $b_1pk_1+p_1k_2=a_1-ab_1$，那么 $k_1,k_2$ 又可以用扩欧来解了。

My Code

#define int long long
using namespace std;
const int MAXN=1e5+10;
int exgcd(int a,int b,int &x,int &y){
	if(b==0){x=1,y=0;return a;}
	int g=exgcd(b,a%b,y,x);
	y-=a/b*x;return g;
}
int a[MAXN],p[MAXN],b[MAXN],r[MAXN];
multiset<int> s;
void solve(){
	int n,m;cin>>n>>m;
	s.clear();
	rep(i,1,n) cin>>a[i];
	rep(i,1,n) cin>>p[i];
	rep(i,1,n) cin>>r[i];
	int atk;
	rep(i,1,m){
		cin>>atk;
		s.insert(atk);
	}
	int mx=0;
	rep(i,1,n){
		auto it=s.upper_bound(a[i]);
		if(it!=s.begin()) b[i]=*(--it);
		else b[i]=*it;
		s.erase(it);s.insert(r[i]);
		mx=max(mx,(a[i]-1)/b[i]+1);
	}
	int ans=0,lcm=1;
	rep(i,1,n){
		int A=(__int128)b[i]*lcm%p[i],B=p[i],C=(a[i]-b[i]*ans%p[i]+p[i])%p[i];
		int x,y,g=exgcd(A,B,x,y);x=(x%B+B)%B;
		if(C%g){cout<<"-1\n";return;}
		ans+=(__int128)(C/g)*x%(B/g)*lcm%(lcm*=B/g);
		ans%=lcm;
	}
	if(ans<mx) ans=ans+((mx-ans-1)/lcm+1)*lcm;
	cout<<ans<<'\n';
}
signed main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	int T;for(cin>>T;T--;)solve();
	return 0;
}

离散对数问题#

这是这样一类问题：求解 $a^x\equiv b\pmod p$。由 $a^{\varphi(p)}\equiv 1\pmod p$，可知必然有 $\varphi(p)$ 以内的解。

BSGS#

复杂度 $O(\sqrt{p})$，求解离散对数的最小自然数解。考虑分块思想，令 $B=\lceil\sqrt{p}\rceil$，那么假如说我们有一个解是 $x$，将其分解成 $kB+j$，那么肯定有 $k,j\le \lceil\sqrt{p}\rceil$。于是我们带入方程：$a^{kB+j}\equiv b\pmod p$，移项得到：$a^j\equiv b\times a^{-kB}\pmod p$。基于此，我们先预处理出 $a^j,j\in[0,B-1]$，存入哈希表中，然后枚举计算 $a^{-0},a^{-B},a^{-2B}\cdots$。每次带入看哈希表中是否有对应的值，然后就能求出解了。

为了 exBSGS 的方便，这里采用一种不需要求逆元的方法。上面的等式其实可以写成：

$$a^{kB}\equiv b\cdot a^j\pmod p$$

那我们把 $b\cdot a^j$ 存入哈希表，然后查询就可以了，最后答案是 $kB-j$。

int BSGS(int a,int b,int p){
	unordered_map<int,int> J;
	//这个哈希表可以用 pbds 代替，可能安全一点
	int A=1,B=(int)ceil(sqrt(p));
	rep(i,0,B) J[A*b%p]=i,A=A*((i==B)?1:a)%p;
	int BA=A;A=1;
	rep(k,0,B){
		if(J.count(A)&&k*B>=J[A]) return k*B-J[A];
		A=A*BA%p;
	}return -1;
}

exBSGS#

如果说 $p$ 不是质数，那么假如 $a$ 不互质了就没法逆元，我们需要转化成 $a,p$ 互质的情况。我们考虑到 $a$ 可以有很多很多次，我们不妨每次把整个式子除以 $\gcd(a,p)$，然后把 $a$ 拆下来，即如果有 $g=\gcd(a,p)$，则：

$$\dfrac{a}{g}\cdot a^{x-1}\equiv \dfrac{b}{g}\pmod{\dfrac{p}{g}}$$

如果 $b$ 不整除 $g$，那就无解。

然后我们看 $\gcd(\dfrac{p}{g},a)$ 是否为 $1$，否则继续做上面这个，直到 $\gcd(\dfrac{p}{\prod g},a)=1$。我们记 $G=\prod_{i=1}^k g_i$，则这个式子变成了：

$$a^{x-k}\equiv \dfrac{b}{G}\cdot(\dfrac{a^k}{G})^{-1}\pmod {\dfrac{p}{G}}$$

然后利用 BSGS，可以求出 $x-k$，再回带就可以了。注意这里的逆元可能需要扩欧。

int exgcd(int a,int b,int &x,int &y){
	if(b==0){x=1;y=0;return a;}
	int g=exgcd(b,a%b,y,x);
	y-=a/b*x;return g;
}
int inv(int a,int p){
	int x,y;exgcd(a,p,x,y);
	x=(x%p+p)%p;return x;
}
int BSGS(int a,int b,int p){
	unordered_map<int,int> J;
	int A=1,B=(int)ceil(sqrt(p));
	rep(i,0,B) J[A*b%p]=i,A=A*((i==B)?1:a)%p;
	int BA=A;A=1;
	rep(k,0,B){
		if(J.count(A)&&k*B>=J[A]) return k*B-J[A];
		A=A*BA%p;
	}return -1;
}
int exBSGS(int a,int b,int p){
	b%=p;
	if(b==1||p==1) return 0;
	int k=0,g,A=1;
	while((g=__gcd(p,a))!=1){
		if(b%g) return -1;
		b/=g;p/=g;(A*=a/g)%=p;
		k++;if(A==b) return k;
	}
	int res=BSGS(a,b*inv(A,p)%p,p);
	if(res==-1) return -1;
	return res+k;
}

• P4861 按钮

得到方程：

$$K^x\equiv 1\pmod M$$

然后直接 exBSGS 就能求出最小解了。

My Code

#define int long long
using namespace std;
int exgcd(int a,int b,int &x,int &y){
	if(b==0){x=1;y=0;return a;}
	int g=exgcd(b,a%b,y,x);
	y-=a/b*x;return g;
}
int inv(int a,int p){
	int x,y;exgcd(a,p,x,y);
	x=(x%p+p)%p;return x;
}
int BSGS(int a,int b,int p){
	unordered_map<int,int> J;
	int A=1,B=(int)ceil(sqrt(p));
	rep(i,0,B) J[A*b%p]=i,A=A*((i==B)?1:a)%p;
	int BA=A;A=1;
	rep(k,0,B){
		if(J.count(A)&&k*B>J[A]) return k*B-J[A];
		A=A*BA%p;
	}return -1;
}
int exBSGS(int a,int b,int p){
	b%=p;
	int k=0,g,A=1;
	while((g=__gcd(p,a))!=1){
		if(b%g) return -1;
		b/=g;p/=g;(A*=a/g)%=p;
		k++;if(A==b) return k;
	}
	int res=BSGS(a,b*inv(A,p)%p,p);
	if(res==-1) return -1;
	return res+k;
}
signed main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	int m,k;cin>>m>>k;
	int res=exBSGS(k,1,m);
	if(res==-1) cout<<"Let's go Blue Jays!\n";
	else cout<<res<<'\n';
	return 0;
}