学习笔记——同余数论
前置芝士
同余概念，逆元，费马小定理，Exgcd，欧拉函数及定理。
（扩展）中国剩余定理
个人喜欢使用非公式形式的合并方式来求解同余方程组。
我们需要求解的问题是求合法的解 \(x\)，使之满足：
此时我们考虑一个人类智慧来合并两个方程组，如果我们能做到这一点，那么我们最终能把方程组变成一个方程，从而得到通解。所以现在我们对：
进行合并。对 \(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 的方便，这里采用一种不需要求逆元的方法。上面的等式其实可以写成：
那我们把 \(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)\)，则：
如果 \(b\) 不整除 \(g\)，那就无解。
然后我们看 \(\gcd(\dfrac{p}{g},a)\) 是否为 \(1\)，否则继续做上面这个，直到 \(\gcd(\dfrac{p}{\prod g},a)=1\)。我们记 \(G=\prod_{i=1}^k g_i\)，则这个式子变成了：
然后利用 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 按钮
得到方程：
然后直接 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;
}
