x^a=b(mod c)求解x在[0,c-1]上解的个数模板+原根求法

/*************************************
 求解x^a=b(mod c) x在[0,c-1]上解的个数模板
 输入：1e9>=a,b>=1,1e9>=c>=3.
 返回：调用xaeqbmodc(a,b,c),返回解的个数
 复杂度: 找原根的复杂度很低，所以总的复杂度为O(c^0.5)
 ************************************/

typedef long long ll;
#define HASH_N 100007

struct hashnode
{
    int next;
    ll key;
    int j;
}HashLink[ HASH_N ];

int hashpre[ HASH_N ],hashcnt;

void hash_insert(ll x,ll key,int j)
{
    for(int p=hashpre[x];p!=-1;p=HashLink[p].next)
    {
        if(HashLink[p].key==key) return ;
    }
    HashLink[ hashcnt ].key = key;
    HashLink[ hashcnt ].j = j;
    HashLink[ hashcnt ].next = hashpre[x];
    hashpre[x] = hashcnt++;
}

int hash_get(ll key)
{
    ll x=key%HASH_N;
    for(int p=hashpre[x];p!=-1;p=HashLink[p].next)
    {
        if( HashLink[p].key == key ) return HashLink[p].j;
    }
    return -1;
}

ll gcd(ll a,ll b)
{
    if(b==0) return a;
    return gcd(b,a%b);
}

//ax + by = gcd(a,b)
//传入固定值a,b.放回 d=gcd(a,b), x , y
void extendgcd(long long a,long long b,long long &d,long long &x,long long &y)
{
    if(b==0){d=a;x=1;y=0;return;}
    extendgcd(b,a%b,d,y,x);
    y-=x*(a/b);
}

//Ax=1(mod M)
//返回x的范围是[0,M-1]
ll GetNi(ll A,ll M)
{
    ll rex=0,rey=0;
    ll td=0;
    extendgcd(A,M,td,rex,rey);
    return (rex%M+M)%M;
}

//a^b%mod 快速幂
long long Quk_Mul(long long a,long long b,long long mod)
{
    long long qsum=1;
    while(b)
    {
        if(b&1) qsum=(qsum*a)%mod;
        b>>=1;
        a=(a*a)%mod;
    }
    return qsum;
}

//测试x较小的情况，必须！
ll firsttest(ll A,ll B,ll C)
{
    ll tmp=1;
    if(B==1) return 0;
    for(int i=1;i<100;i++)
    {
        tmp = (tmp*A)%C;
        if(tmp==B) return i;
    }
    return -1;
}

ll BabyStep(ll A,ll B,ll C,ll OC)
{
    if(0==A || 0==C) return -1;
    if(C==1) return 0;
    B = B%C;
    ll ans = firsttest(A,B,C);//为了防止x比较小的时候
    if(ans != -1) return ans;
    ll D=1;
    int c=0;
    ll d;
    while( (d=gcd(A,C)) != 1 )
    {
        if( B%d !=0 ) return -1;//无解的情况
        B /= d;
        C /= d;
        D = D*A/d%C;
        c++;
    }
    
    //得到了 D*A^(x-c)=B (mod C) ，gcd(A,C)=1 , gcd(D,C)=1
    ll D_1=GetNi(D,C);//求D的逆元
    B = B*D_1%C;
    //求A^x=B (mod C),然后返回x+c
    ll m = ceil( sqrt(C+0.0) );
    
    memset(hashpre,-1,sizeof(hashpre));
    hashcnt=0;
    ll hashnum=1;
    hash_insert(1, 1, 0);
    for(int i=1;i<m;i++)
    {
        hashnum = (hashnum*A)%C;
        hash_insert(hashnum%HASH_N, hashnum ,i);
    }
    
    ll ol=OC;//这一步可以省略
    ll mA=Quk_Mul(A, m, C);
    ll ta=1;
    
    ll tmpni = Quk_Mul(mA, ol-1, C);
    
    for(int i=0;i<m;i++,ta=ta*tmpni%C)
    {
        //解线性同余方程 tx=B*ta (%C) ,ta直接用费马小定理求逆元
        ll tx = ta;
        tx = (tx*B)%C;
        int j=hash_get(tx);
        if(j!=-1)//找到解了
        {
            return i*m+j+c;
        }
    }
    
    return -1;
}

ll mypow(ll a,ll b)
{
    ll sum=1;
    for(int i=1;i<=b;i++)
        sum*=a;
    return sum;
}

ll slove(ll a,ll b,ll c,ll p,int a1)
{
    //第0种情况a==1
    if(a==1) return 1;
    
    //第一种情况b==c
    b %= c;
    if(b==0)
    {
        ll tmp = mypow(p,ceil( (double)a1/a ) );
        return c/tmp;
    }
    ll d=gcd(b,c);
    //第二种情况 gcd(b,c)!=1
    if( d != 1 )
    {
        return 0;
    }
    
    //第三种情况 gcd(b,c)==1
    
    //第一步找出原根x0,
    ll save[55];
    int cnt=0;
    
    ll tc = mypow(p,a1-1)*(p-1);
    for(ll i=2;i*i<=tc;i++)
    {
        if(tc%i == 0)
        {
            save[ cnt++ ] = i;
            while(tc%i==0) tc/=i;
        }
    }
    if(tc != 1)
    {
        save[ cnt++ ] = tc;
    }
    tc= mypow(p,a1-1)*(p-1);
    for(int i=0;i<cnt;i++)
    {
        save[i] = tc/save[i];
    }
    
    ll x0=0;
    for(int i=2;i<c;i++)
    {
        int flag=0;
        for(int j=0;j<cnt;j++)
        {
            if( Quk_Mul(i, save[j], c) ==1)
            {
                flag=1;
                break;
            }
        }
        if(flag==0)
        {
            x0=i;
            break;
        }
    }
    //找到原根x0后,然后找x0^a0 = b (mod c)
    ll a0 = BabyStep(x0 , b, c,mypow(p,a1-1)*(p-1));
    ll ord = mypow(p,a1-1)*(p-1);
    d = gcd(a,ord);
    if( a0%d != 0 ) return 0;
    return d;
}

ll xaeqbmodc(ll a,ll b,ll c)
{
    ll ans=1;
    b%=c;
    ll tc=c;
    //然后对c进行因式分解
    for(ll i=2;i*i<=tc;i++)
    {
        if(tc%i==0)
        {
            ll tmpyz=1;
            int a1=0;
            while(tc%i==0)
            {
                a1++;
                tmpyz *= i;
                tc/=i;
            }
            //然后对这个因子进行处理
            ans *= slove(a,b,tmpyz,i,a1);
            if(ans==0) break;
        }
    }
    if(tc!=1)
    {
        ans *= slove(a,b,tc,tc,1);
    }
    return ans;
}