模拟赛 T1 费马小定理+质因数分解+exgcd

求：$a^{bx \%p}equiv 1(mod p)$ 的一个可行的 $x$.

根据欧拉定理，我们知道 $a^{phi(p)}equiv 1(mod p)$

而在 $a^xequiv 1(mod p)$ 这个式子中 $x$ 是存在很多个解的.

这些解之间存在着循环节，使得任意解 $x$ 可以被表示成循环节的倍数.

我们设这个循环节为 $cir$.

由于已知 $phi(p)$ 一定是一个可行解，所以最小循环节一定是 $phi(p)$ 的约数.

然后我们就可以对 $phi(p)$ 进行质因数分解来求这个最小循环节 $cir$.

求出来 $cir$ 后，可得 $bx\%p=cir imes y$

你发现因为有这个 $\%p$ 操作，所以会导致这个方程解不出来.

但是好在我们发现，左面那个式子可以变为 $bx+pk=l$，而这个 $l$ 可以被表示为 $i imes gcd(b,p)$

故我们可以将式子变为 $i imes gcd(b,p)=y imes cir$.

然后我们可以对这个求最小通解（因为如果最小解小于 $p$，则一定可以被 $bx\%p$ 表示出来，而且满足 $xleqslant p-1$ )

这个最小通解是 $c=frac{cir imes gcd(b,p)}{gcd(gcd(b,p),cir)}$

然后用 exgcd 求一下 $bx+pk=c$ 的 $x$ 的最小正整数解就可以了.

这里可以证明一下为什么只要存在 $bx+pk=c$ 就能保证 $x<p$：

我们可以将 $x$ 表示成 $p+d$ 的形式，那么原式为 $b(p+d)+pk=c$

$Rightarrow bp+bd+pk=c$

$Rightarrow b imes d+p imes (k+d)=c$

所以，一旦 $x>p$，我们就可以将一些部分导到 $p imes k$ 那里，以此来实现 $x<p$   

#include <string>  
#include <ctime>
#include <cmath> 
#include <cstdio> 
#include <vector>
#include <algorithm>  
#define N 10000060  
#define ll long long  
using namespace std;     
void setIO(string s) {
    string in=s+".in"; 
    string out=s+".out"; 
    freopen(in.c_str(),"r",stdin); 
    freopen(out.c_str(),"w",stdout); 
}
int H=0;  
int fr[N];  
int prime[N],vis[N],phi[N]; 
int answer[N];
int MN[N];      
struct data {
    int a,b,p,id;      
    data(int a=0,int b=0,int p=0,int id=0):a(a),b(b),p(p),id(id){} 
};   
vector<data>G[N];  
inline int qpow(int x,int y,int mod) { 
    int tmp=1;  
    while(y) {       
        if(y&1) {   
            tmp=(ll)tmp*x%mod; 
        }
        x=(ll)x*x%mod; 
        y>>=1; 
    }
    return tmp; 
}  
int exgcd(int a,int b,int &x,int &y) {
    if(!b) {
        x=1,y=0;   
        return a;   
    }     
    int gcd=exgcd(b,a%b,x,y);       
    int tmp=x;    
    x=y,y=tmp-(a/b)*y;           
    return gcd;  
} 
void Linear_shaker() {        
    int i,j,cnt=0;    
    fr[1]=1;   
    for(i=2;i<N;++i) {
        if(!vis[i]) {
            prime[++cnt]=i;   
            phi[i]=i-1;  
            fr[i]=1;   
        }
        for(j=1;j<=cnt&&prime[j]*i<N;++j) {
            vis[i*prime[j]]=1;   
            if(i%prime[j]) {
                phi[i*prime[j]]=phi[i]*(prime[j]-1);   
                fr[i*prime[j]]=i;  
            }
            else {
                fr[i*prime[j]]=i;  
                phi[i*prime[j]]=phi[i]*prime[j]; 
                break; 
            }
        }
    }   
}
int main() {   
    // setIO("mod");        
    Linear_shaker();  
    int i=0,j=0,t1,t2,tp=0;                    
    int a,b,p,Mx=0,k;        
    while(scanf("%d",&a)!=EOF) { 
        scanf("%d%d",&b,&p);  
        ++i;                  
        ++tp;             
        if(__gcd(a,p)!=1) {   
            answer[i]=-1;    
        }
        else {
            Mx=max(Mx,phi[p]);       
            G[phi[p]].push_back(data(a,b,p,i));  
            MN[i]=phi[p];    
            int w=phi[p],tmp,c=1;   
            while(w!=1) {
                while(w!=1&&qpow(a,fr[w]*c,p)==1) {   
                    w=fr[w];   
                }                            
                tmp=w/fr[w];          
                while(w%tmp==0&&w!=1) {    
                    c*=tmp;  
                    w/=tmp;   
                }
            }
            MN[i]=c;  
        }
    }                     
    for(i=1;i<=Mx;++i) {
        for(j=0;j<G[i].size();++j) { 
            a=G[i][j].a; 
            b=G[i][j].b;  
            p=G[i][j].p;   
            int x=0,y=0; 
            int id=G[i][j].id;  
            int delta=MN[id]; 
            int gcd=__gcd(p,b);      
            int tmp=(1ll*gcd*delta)/(__gcd(gcd,delta));    
            if(tmp>=p) {
                answer[id]=-1; 
            }
            else {  
                gcd=exgcd(b,p,x,y);  
                x=(1ll*x*(tmp/gcd)%(p/gcd)+(p/gcd))%(p/gcd);  
                answer[id]=x;                       
            }                 
        }
    }       
    for(i=1;i<=tp;++i) printf("%d
",answer[i]);         
    return 0; 
}