POJ3243 clever Y

求\(a^x\equiv b(\mod p)\)的最小整数解，其中p不一定为质数。

普通的BSGS只能解决a,p互质的问题，当a,p不互质怎么办呢？ 扩展BSGS的思想：

每次取一个a，如果\(\gcd(a,p)\not = 1\)，则进行判定:如果\(b\not | \gcd(a,p)\)则无解，否则b和p全部\(/\gcd(a,p)\)。 取到\(\gcd(a,p)=1\)时结束,记录取的次数s。 取时记录d为每一次\(\frac{a}{\gcd(a,p)}\)的积。 取完以后b要\(/ d\)，也就是说新的b是原来的\(b/a^s\)

接着跑BSGS就行，计算出的答案要+取的次数。 这样就能算出一个\(\ge s\)的答案，但是答案如果\(\le s\)怎么办呢？

因为每次p最少\(/2\)，所以\(s \le log_2{p}\)，显然\(s \le 50\)，所以一开始先计算x在\([0,50]\)之间有没有可能的解。

注意\(b/d\)时，计算逆元不能用费马小定理！要用扩欧！！！因为p不一定是个质数。 代码：(一开始做的是数据水p为质数的BZOJ1467)

/*
Author: CNYALI_LK
LANG: C++
PROG: 1467.cpp
Mail: cnyalilk@vip.qq.com
*/
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<map>
#define debug(...) fprintf(stderr,__VA_ARGS__)
#define DEBUG printf("Passing [%s] in LINE %lld\n",__FUNCTION__,__LINE__)
#define Debug debug("Passing [%s] in LINE %lld\n",__FUNCTION__,__LINE__)
#define all(x) x.begin(),x.end()
using namespace std;
const double eps=1e-8;
const double PI=acos(-1.0);
typedef long long ll;
template<class T>ll chkmin(T &a,T b){return a>b?a=b,1:0;}
template<class T>ll chkmax(T &a,T b){return a<b?a=b,1:0;}
template<class T>T sqr(T a){return a*a;}
template<class T>T mmin(T a,T b){return a<b?a:b;}
template<class T>T mmax(T a,T b){return a>b?a:b;}
template<class T>T aabs(T a){return a<0?-a:a;}
#define min mmin
#define max mmax
#define abs aabs
typedef pair<ll,ll> pii;
#define x first
#define y second
pii aa[2333333];
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll find(ll l,ll r,ll s){
    ll mid;
    while(l<=r){
        mid=(l+r)>>1;
        if(aa[mid].x<=s)l=mid+1;
        else r=mid-1;
    }
    --l;
    if(aa[l].x==s)return aa[l].y;
    else return -1;
}
void exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){x=1;y=0;return;}
    exgcd(b,a%b,y,x);
    y=y-(a/b)*x;
}
ll inv(ll a,ll p){
    ll x,y;
    exgcd(a,p,x,y);
    x=(x%p+p)%p;
    return x;
}
ll fpm(ll a,ll b,ll p){
    ll c=1;
    while(b){
        if(b&1)c=c*a%p;
        a=a*a%p;
        b>>=1;
    }
    return c;
}
ll exbsgs(ll a,ll p,ll b){
    b%=p;

    for(ll i=0,j=1%p;i<=50;++i,j=j*a%p){
        if(j==b)return i;
    }
    ll s=0,g,f=1%p;
    while((g=gcd(a,p))>1){
//      printf("!%lld\n",g);
        if(b%g){return -1;}
        b/=g;
        p/=g;
        f=f*a/g%p;
        ++s;
    }
    b=b*inv(f,p)%p;
    a%=p;

    ll m=ceil(sqrt((long double)p));
    for(ll i=1;i<=m;++i){
        b=b*a%p;
        aa[i]=make_pair(b,i);
    }
    sort(aa+1,aa+m+1);
//  for(ll i=1;i<=m;++i)printf("%lld%c",aa[i].x,i%10?' ':'\n');
    ll c=1,fs;
    for(ll i=1;i<=m;++i)c=c*a%p;
    b=1;
    for(ll i=1;i<=m;++i){
        b=b*c%p;
        ll ss=find(1,m,b);
        if(ss+1)return i*m-ss+s;
    }
    return -1;
}
int main(){
#ifdef cnyali_lk
    freopen("1467.in","r",stdin);
    freopen("1467.out","w",stdout);
#endif
    ll x,z,k,ans;
    while(cin>>x>>z>>k&&(x||z||k)){
        ans=exbsgs(x,z,k);
        if(ans+1)
        printf("%lld\n",ans);
        else printf("No Solution\n");
    }
    return 0;
}