Luogu4195 【模板】exBSGS（exBSGS）

　　如果a和p互质，用扩欧求逆元就可以直接套用普通BSGS。考虑怎么将其化至这种情况。

　　注意到当x>=logp时gcd(a^x,p)是一个定值，因为这样的话每个存在于a中的质因子，其在a^x中的出现次数一定比在p中的多。

　　于是对x<logp的情况暴力验证。对x>=logp的情况，设d=gcd(a^x,p)，剩下的问题变为求a^x/d≡b/d(mod p/d)，这里a^x和p/d显然就是互质的了。

　　要求解这个方程，显然不能把d直接乘过去（好像也说不清为啥）。首先b%d>0时无解。然后考虑从a^x中分离一部分，使该部分能整除d，再将该部分除以d后移到式子右边。直接分离的话会爆long long，每次分离一个a即可。剩下的就是普通BSGS了。

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<map>
using namespace std;
#define ll long long
int gcd(int n,int m){return m==0?n:gcd(m,n%m);}
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
int a,p,b;
map<int,int> f;
void exgcd(int a,int b,int &x,int &y)
{
    if (b==0)
    {
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    int t=x;x=y;y=t-a/b*x;
}
int inv(int a)
{
    int x,y;exgcd(a,p,x,y);
    x=(x%p+p)%p;
    return x;
}
int BSGS(int a,int b,int p)
{
    int block=sqrt(p),t=1;//cout<<a<<' '<<b<<' '<<p<<endl;
    f.clear();
    for (int i=0;i<block;i++)
    {
        if (f.find(t)==f.end()) f[t]=i;
        if (t==b) return i;
        t=1ll*t*a%p;
    }
    int v=t;
    for (int i=1;i<=(p-1)/block;i++)
    {
        if (f.find(1ll*b*inv(t)%p)!=f.end()) return i*block+f[1ll*b*inv(t)%p];
        t=1ll*t*v%p;
    }
    return -32;
}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("exbsgs.in","r",stdin);
    freopen("exbsgs.out","w",stdout);
#endif
    a=read(),p=read(),b=read();
    while (a)
    {
        if (p==1) printf(b==0?"0
":"No Solution
");
        else if (b==1) printf("0
");
        else
        {
            int t=1,ans=0;
            for (int i=1;i<=31;i++)
            {
                t=1ll*t*a%p;
                if (t==b) {ans=i;break;}
            }
            if (!ans)
            {
                int u=gcd(p,t);
                if (b%u==0)
                {
                    //a^x/u=b/u (%p/u)
                    int P=p;b/=u;p/=u;
                    for (int i=1;i<=31;i++)
                    if (P==p) {ans=i-1;break;}
                    else
                    {
                        int u=gcd(a,P);
                        b=1ll*b*inv(a/u)%p;
                        P/=u;
                    }
                    ans+=BSGS(a,b,p);ans=max(ans,0);
                }
            }
            if (ans) printf("%d
",ans);
            else printf("No Solution
");
        }
        a=read(),p=read(),b=read();
    }
    return 0;
}