洛谷 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)\)

\[a'^xg^x+yp'g=b \]

\[g(a'^xg^{x-1}+yp')=b \]

这样子\(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