upd: 应评论区,我自己查了一些资料之后发现原来写的东西确实有些不对的地方,造成的误解请各位海涵,本人蒟蒻,如果对于某个描述还是有疑问的话欢迎继续指出,谢谢各位 ^_^
特别鸣谢
NaVi_Awson 学长 以及
落汐 大佬/qq
目录
- 概念讲解
- 基础算法-解单个线性同余方程
- 进阶-解模数互质的线性同余方程组
- 再进阶-解一般的线性同余方程组
- 入坑-解未知数在指数位置且模数为质数的同余方程
- 入土-解未知数在指数位置的一般模数同余方程
(为什么没了呢?因为我不会了TAT)
概念
1.带余除法
实际上就是小学学的那种“
3 除以
2 商
1 余
1”的除法,书面一点的描述就是
a=q×b+r 其中
a,b,q,r 都是整数,
a 叫做被除数,
b 叫做除数
(b=0),
q 是商,
r 是余数(
0≤r<b),C++中,整型之间的除法是整数除法,设
a=q×b+r ,那么
a/b=q
本文中的除法都是整数除法
2.同余方程
同余方程就是长成这样的柿子:
ax≡c(modb)
这个式子什么意思? 就是字面意思,同余即指余数相同(
axmodb=cmodb)
算法
Day 1.Exgcd (扩展欧几里得算法)(求解二元一次不定方程/单个线性同余方程)
不定方程
ax+by=c
(或者
ax≡c(modb) 设
ax=c+(−y)×b 移项可得不定方程
ax+by=c )
我们先看一下简单一点的一个方程
ax+by=gcd(a,b)
证明 gcd(a,b)=gcd(b,amodb)
设
g=gcd(a,b)
则
a=g×k1 ,
b=g×k2,
则
amodb=(g×k1)mod(g×k2)=(k1modk2)×g,
b=g×k2
为了使
gcd(a,b)=gcd(b,amodb)
k1modk2 与
k2 应该没有公因子
假设
k1modk2 与
k2 有公因子为
r
则
k1modk2=n1×r,k2=n2×r
所以
k1=x×k2+n1×r=(x×n2+n1)×r
此时
gcd(k1,k2)=r 不符合
gcd(a,b)=g (否则
gcd(a,b)=r×g)
所以
k1modk2 与
k2 没有公因子
即
gcd(a,b)=gcd(b,amodb)
证毕
证明 当a,b不全为0时ax+by=gcd(a,b)一定有解
已证:
gcd(a,b)=gcd(b,amodb)
所以我们可以得到推论:
ax+by=gcd(a,b) 可以转化为
bx+(amodb)y=gcd(b,amodb)
所以当
b=0 时(欧几里得算法最后一步)
方程变成
ax=gcd(a,b)=a
所以一定有一组解
x=1,y=0
证毕
利用这个方法我们可以递归到边界(
b=0),然后将
x 设为
1,
y 设为
0
然后在递归返回的时候更新
x 和
y 的值
由于方程
ax+by=gcd(a,b) 转化为了方程
bx+(amodb)y=gcd(b,amodb)
而
amodb=a−b×(a/b)除号当然是整除啦
方程左边
bx+(a−b×(a/b))y=ay+(x−(a/b)×y)b
所以转化后新的
x=y,新的
y=x−(a/b)×y
但是代码实现可以更简单一些,利用传引用在访问下一层递归的时候就交换
x 和
y
Code
CPPint g,x,y;
void Exgcd(int a,int b,int &x,int &y)
{
if(!b)
{
x=1;y=0;g=a;
return ;
}
Exgcd(b,a%b,y,x);
y-=(a/b)*x;
return ;
}
那么解决了方程
ax+by=gcd(a,b) 之后
ax+by=c 也很好解决了
我们解出
ax+by=gcd(a,b) 后,如果
c 不是
gcd(a,b) 的整数倍,那么这个同余方程是没有解的。
如果有解,那么解为
x×(c/gcd(a,b))
通解为
xi=x×(c/gcd(a,b))+(b/gcd(a,b))×k (
k∈Z)
利用 Exgcd 可以解决较为基本的不定方程和同余方程问题
时间复杂度
O(logn) (证明可以自己去了解)
Day 2.CRT(中国剩余定理/孙子定理)(求解模数两两互质的线性同余方程组)
(线性的意思就是一次的方程)
同余方程组长这样:
⎩⎨⎧x≡r1(modm1)x≡r2(modm2)⋯x≡ri(modmi)⋯x≡rn(modmn)
其中
m1,m2,⋯,mn 两两互质
中国剩余定理:对于如上的线性同余方程组,若其中
m1,m2,⋯,mn 两两互质
设
M=i=1∏nmi,
Mi=M/mi
设
Ki 为
MiKi≡1(modmi)的解
则对于任意的
r1,r2,⋯,rn 方程组有整数解
x=i=1∑nriMiKi
证明中国剩余定理
因为
Mi=M/mi 是除了
mi 以外的所有的
m 的倍数(模数两两互质)
所以
∀k=i,
riMiKi≡0(modmk)(都说了是倍数那么右边肯定是
0 了)
又因为
Ki 为
MiKi≡1(modmi)的解,左右同时乘以
ri 得到
riMiKi≡ri(modmi)
所以当
x=i=1∑nriMiKi 时,对于任意的
x≡ri(modmi),因为
∀k=i,
riMiKi≡0(modmk) 而且
riMiKi≡ri(modmi),所以方程成立。
即方程组的解为
x=i=1∑nriMiKi
证毕
当题目保证给出的模数两两互质的时候,
CRT it is!
步骤:
- 求出M=i=1∏nmi
- 利用 Exgcd 求出 MiKi≡1(modmi)的解 Ki
- 逐步累加 riMiKi (即利用 x=i=1∑nriMiKi)得到解(模 M 意义下)
Code:
CPPtypedef long long LL;
LL CRT()
{
LL pro=1;
for(LL i=1;i<=N;i++)
pro*=m[i];
LL x,y,Mi;LL ans=0;
for(LL i=1;i<=N;i++)
{
Mi=pro/m[i];
Exgcd(Mi,m[i],x,y);
ans=(ans+Mi*r[i]*x)%pro;
}
return (ans+pro)%pro;
}
时间复杂度
O(nlogn)
Day 3.ExCRT(扩展中国剩余定理/孙子定理)(求解一般的线性同余方程组)
ExCRT 实际上
superset 了
CRT. By connect
不是所有的方程组模数都两两互质(
正如不是所有牛奶...算了)
那么现在所有的
m 之间都没有什么关系了
这样的话方程
MiKi≡1(modmi) 不一定有解,无法通过
x=i=1∑nriMiKi 得到解
好的那么现在你面对的是?
⎩⎨⎧x≡r1(modm1)x≡r2(modm2)⋯x≡ri(modmi)⋯x≡rn(modmn)
Q:现在你会什么? A:解一个方程
假设我们解出来了第一个方程,将其解设为
ans ,我们知道第一个方程的通解
xv=ans+k×m1,k∈Z
我们现在的目的就是要求一个同时满足第一第二两个方程的
x ,那么我们把第一个方程的通解代入第二个方程
ans+k×m1≡r2(modm2) 移项得
k×m1≡r2−ans(modm2)
即求解不定方程
m1×k+m2×y=r2−ans
解出来
k 之后得到第一第二两个方程的公共解
ans′=ans+k×lcm(m1,m2)
那么前两个同余方程等价为:
x≡ans′(modlcm(m1,m2))
为什么等价?
因为向方程中代入
xv=ans′+k×lcm(m1,m2),k∈Z
你会发现这两个方程和上面的等价方程同时成立,所以这里可以等价
那么你就会发现一件事情,你成功地把两个方程整成了一个方程
那么下一步,继续从方程组里面拿一个方程出来,这时候你手上又有了两个方程
{x≡ans(modlcm)x≡ri(modmi)
由于第一个方程的通解是
xv=ans+k×lcm,k∈Z
所以我们就是要求一个满足条件的
xv 使
xv≡ri(modmi) 成立
所以就是求解
ans+k×lcm≡ri(modmi)
变形得
k×lcm≡ri−ans(modmi)
即求解不定方程
lcm×k+mi×y=ri−ans
如果方程无解,则同余方程组无解
方程有解,那么解
ans′=ans+k×lcm
循环求解即可。
形式化描述过程就是:
设
ans 为前
i−1 个同余方程的公共解,
lcm 是前
i−1 个方程模数的最小公倍数
那么对于第
i 个方程我们求解不定方程
lcm×k+mi×y=ri−ans
不定方程无解则同余方程组无解
得到
k 后我们知道了前
i 个同余方程的公共解
ans′=ans+lcm×k
前
i 个方程模数的最小公倍数
lcm′=lcm×(mi/gcd(lcm,mi))
此时可求解第
i+1 个方程,循环即可。
CPPinline LL ExCRT()
{
LL ans=0,lcm=1;
LL k1,k2,c;
for(int i=1;i<=N;i++)
{
Exgcd(lcm,m[i],k1,k2);
c=((r[i]-ans)%m[i]+m[i])%m[i];
if(c%g) return -1;
k1=((k1*c/g)%(m[i]/g)+(m[i]/g))%(m[i]/g);
ans=ans+k1*lcm;
lcm=lcm*m[i]/g;
ans=(ans%lcm+lcm)%lcm;
}
return (ans%lcm+lcm)%lcm;
}
时间复杂度
O(nlogn) (常数比
CRT 大一点)
Day 4.BSGS(BabyStepGiantStep)(解未知数在指数位置且模数为质数的同余方程)
有时候事情并不会都是那么简单,经过之前几个小节内容的学习你应该已经可以熟练地解决线性同余方程的相关问题
但是啊
Ax≡B(modC)
其中
A,
B,
C 已知 且
A,C 互质,需要求解最小的
x 这个该怎么办呢?
BSGS it is!
首先我们想一下对于这个问题暴力怎么搞
当然直接枚举指数啊qwq
但是这样枚举的尽头是什么呢?
首先我们知道
A0≡1(modC) ,然后又由费马小定理(不知道这个的同学可以先去了解一下逆元再了解费马小定理)
AC−1≡1(modC)
所以我们发现实际上枚举指数的时候枚举
C 个就会开始重复,那么暴力枚举
[0,C−1) 即可。
CPPint a,m,r;
scanf("%d%d%d",&a,&m,&r);
int base=1;
for(int i=0;i<m-1;i++)
{
if(base==r) {printf("%d",i);return 0;}
base=(base*a)%m;
}
printf("no solution");
return 0;
但是假如现在
C=19260817 但是只给你
200ms 这个暴力算法是会超时的(而且题目要是想卡暴力质数可以出的很大)
那么我们要对暴力算法进行改进,核心肯定是缩小枚举的范围
这就是
BSGS 的雏形了
设
x=i×C−j
对于同余方程
Ai×C−j≡B(modC)
左右两边同时除以
B 得
Aj×BAi×C≡1(modC)
预处理分母,然后当模
C 意义下有
Ai×C=Aj×B 时
Aj×BAi×C=1, 方程有解
i×C−j
分母怎么预处理呢?
因为我们的解形式是
i×C−j ,所以
j∈[0,C)
预处理出
(Aj×B)modC,j∈[0,C) ,将其幂指数存在哈希表里面
暴力枚举
i,i∈[1,C] ,对于每一个
Ai×C 在哈希表里面查询是否有对应的
Aj×B
若存在那么最小解就是
i×C−哈希表中存储的幂指数
那么现在是不是能理解
BabyStepGiantStep 这个名字了呢?
我们在预处理哈希表的时候就是在走
BabyStep (一次加
1)
枚举
i×C 中的
i 时就是在走
GiantStep (一次加
C )
CPPLL BSGS()
{
Hash.clear();
LL D=B%C,sqrtm=ceil(sqrt(C));
for(int i=0;i<sqrtm;i++)
{
Hash.insert(D,i);
D=D*A%C;
}
LL t=ti=fast_pow(A,sqrtm),ans;
for(int i=1;i<=sqrtm;i++)
{
if((ans=Hash.find(ti))!=-1)
return i*sqrtm-ans;
ti=ti*t%C;
}
return -1;
}
时间复杂度
O(C) (快速幂的
logn 在
C 下忽略不计)
Day5.ExBSGS(ExtendedBabyStepGiantStep)(解未知数在指数位置的一般模数同余方程)
就像我们在
ExCRT 的时候说过的一样,不是所有的模数都那么美好让
A,C 互质
不互质的时候分式
Aj×BAi×C 会因为
Aj 没得逆元导致无意义
我们只会互质的解法,怎么办?
没有条件创造条件!
对于同余方程
Ax≡B(modC),A,C 不互质
那么我们让他互质!
设
g=gcd(A,C) ,我们需要对方程执行消去因子
设
A=a×g,C=c×g 我们必须要把方程中的幂底数和模数化为互质的,如果途中发现
Bmodg=0 那么这个方程无解
证明
若A=a×g,C=c×g,g=1.而Bmodg=0,方程无解
设
D=Ax−1 即有同余方程
A×D≡B(modC)
化为不定方程
ax+by=c 形式:
A×D+C×y=B
在学习
Exgcd 的时候我们就已经知道了 上述方程有解的充分必要条件就是
B 可以被
gcd(A,C) 也就是
g 整除,而现在
Bmodg=0 所以方程无解
证毕
为了使模数互质我们执行消去因子达到
gcd(A,C)=1 为止
依据
BSGS 的经验我们再来
对于方程
Ax≡B(modC)
设
A=a×g,C=c×g,g=1 ,全除以
g 得到
a×(a×g)x−1≡B/g(modC) 假如
gcd(A,c)=1 继续消去因子
最后可以得到这样一个柿子 :
D×Ax−cnt≡B′(modC′)
其中的
D 是每消去一次因子左边的值除掉一个
g 后得到的商的乘积 (即
D=i=1∏cnt(A/gi),gi为第
i 次消去因子时的
gcd(A,C) )
设
x=i×C′−j+cnt ,两边同时除以
B′ 得
Aj×B′D×Ai×C′≡1(modC′)
现在
A,C′ 互质了,套上
BSGS 求解出对应的
i×C′−j 然后加上
cnt 即可
但是这样得到的
x>cnt 我们需要对小于
cnt 的解特判,消去因子的时候做即可
CPPinline LL BSGS()
{
LL sqrtm=ceil(sqrt(C));
H.clear();
LL t=B%C;
for(int i=0;i<sqrtm;i++)
{
H.insert(t,i);
t=t*A%C;
}
t=fast_pow(A,sqrtm,C);
LL ti=D*t%C,ans;
for(int i=1;i<=sqrtm;i++)
{
if((ans=H.find(ti))!=-1&&i*sqrtm-ans!=0)
return i*sqrtm-ans;
ti=ti*t%C;
}
return -1;
}
inline LL ExBSGS()
{
LL cnt=0;D=1;
if(B==1) return 0;
while((g=gcd(A,C))!=1)
{
if(B%g) return -1;
B/=g;C/=g;
D=(D*A/g)%C;
cnt++;
if(D==B) return cnt;
}
LL ans=BSGS();
if(ans==-1) return -1;
else return ans+cnt;
}
时间复杂度
O(C) (消去因子部分复杂度是
log 级别,忽略)
Congratulations!Learning Completed
如果还有任何对算法不理解的问题珂以私信窝,力所能及的范围内窝一定尽力解答 ^_^