洛谷 P4195 【模板】exBSGS/Spoj3105 Mod
已知数a,p,b，求满足a^x≡b(mod p)的最小自然数x。
扩展BSGS的板子题
回忆\(BSGS\)算法，给定整数\(a,b,p\)，其中\(a,p\)互质，求方程\(a^x\equiv b\ (mod\ p)\)的最小整数解\(x\)
做法：设\(x=i\times m-j,m=\left \lceil \sqrt p \right \rceil,1\le j\le m,1\le i\le m\)
方程变为\(a^{im-j}\equiv b\ (mod\ p)\)
化一下\((a^m)^i\equiv b\times a^j\ (mod\ p)\)
这时只要把\((a^m)^i\)和\(b\times a^j\)预处理出来丢到哈希表里就做完了
而现在的\(a,p\)不互质了，所以我们要考虑其他的方法，也就是\(\text{Ex\_BSGS}\)
那么我们设\(g=gcd(a,p)\)
根据模的分配率，方程变为\(\frac{a^x}{g}\equiv \frac{b}{g}\ (mod\ \frac{p}{g})\)
无解情况就是\(g\nmid b\)并且\(b\ne 1\)
我们来证明一下
设\(a'=\frac{a}{g},p'=\frac{p}{g}\)，那么\(a=a'g,p=p'g\)
代入到原方程中变为\((a'g)^x\equiv b\ (mod\ p'g)\)
这样子\(g\)就是\(b\)的因子，而只有在\(b=1\)时，\(a^0=1\)，其余情况若\(g\nmid b\)，方程无解
证毕
那么我们再把上面的式子化一下变为\(a^{x-1}\times\frac{a}{g}\equiv \frac{b}{g}\ (mod\ \frac{p}{g})\)
而\(\frac{p}{g}\)一定是比\(p\)小的，所以可以一直约到\(a,\frac{p}{g}\)互质
设\(na=\prod_{i=1}^k\frac{a}{g_i}\)
原式就可以写成\(a^{x-k}\equiv \frac{b}{\prod_{i=1}^k g\times na}(mod\ \frac{p}{\prod_{i=1}^kg})\) ，\(k\)是做了几次化简
这样子就可以用\(BSGS\)求解啦
有一点需要注意，因为是在模意义下运算，所以除以\(na\)时要乘其逆元
因为最后一直化到了\(\frac{a}{\prod_{i=1}^kg}\)和\(\frac{p}{\prod_{i=1}^kg}\)互质，所以\(na\)和\(\frac{p}{\prod_{i=1}^kg}\)看起来也很互质
就可以用扩欧求逆元了
Code
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <map>
#include <cmath>
#define int long long
using namespace std;
int a,p,b;
map <int,int> f;
int gcd(int a,int b) //最大公约数
{
if (!b)return a;
return gcd(b,a%b);
}
void exgcd(int a,int b,int &x,int &y) //扩欧
{
if (!b)x=1,y=0;
else
{
exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
}
}
int inv(int a,int b) //逆元
{
int x,y;
exgcd(a,b,x,y);
return (x%b+b)%b;
}
int mypow(int a,int x,int p)
{
int s=1;
while (x)
{
if (x&1)s=s*a%p;
a=a*a%p;
x>>=1;
}
return s;
}
int bsgs(int a,int b,int p) //BSGS算法
{
f.clear();
int m=ceil(sqrt(p));
b%=p;
for (int i=1;i<=m;i++)
{
b=b*a%p;
f[b]=i;
}
int tmp=mypow(a,m,p);
b=1;
for (int i=1;i<=m;i++)
{
b=b*tmp%p;
if (f[b])return (i*m-f[b]+p)%p;
}
return -1;
}
int exbsgs(int a,int b,int p)
{
if (b==1||p==1)return 0; //特殊情况，x=0时最小解
int g=gcd(a,p),k=0,na=1;
while (g>1)
{
if (b%g!=0)return -1; //无法整除则无解
k++;b/=g;p/=g;na=na*(a/g)%p;
if (na==b)return k; //na=b说明前面的a的次数为0，只需要返回k
g=gcd(a,p);
}
int f=bsgs(a,b*inv(na,p)%p,p);
if (f==-1)return -1;
return f+k;
}
signed main()
{
cin>>a>>p>>b;
while(a||b||p)
{
a%=p;b%=p;
int t=exbsgs(a,b,p);
if (t==-1)cout<<"No Solution"<<endl;
else cout<<t<<endl;
cin>>a>>p>>b;
}
return 0;
}
\(BSGS\)写挂了调了半天QAQ
