题解 P4195 【【模板】exBSGS/Spoj3105 Mod】

eee_hoho · 2019-06-27 09:37:13 · 题解

扩展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;
}