hdu 3930 Broot x^k=b%m 已知k,b,m,求所有小于m的解x m是素數

Problem Description There is a way of encryption in the Digital planet. If the number x, who lives in M area, K layer, then it will change its name to x ^ K mod M, that is newx = x ^ k mod m. Now the Digital Kingdom wants to make a program, which can find all the original number of a name living in some area, some layer. If there is no solution, output -1.  

Input There are multiply test cases. Each test case contains there integers k, m , newx ,(0 <= newx , m ,k <= 1.5*10^15) , m is a prime number.  

Output For each test case, you should output some lines, the format of the first line is: “caseN:” (N is the label of test case). Then next lines each contain a number x, output x as ascending order. (0 <= x < m)  

Sample Input

1 5 4
2 13 8
3 13 8
        

Sample Output

case1:
4
case2:
-1
case3:
2
5
6
        

//

題目要求我們求出x^a mod b = c 其中所有 x 的值。因為其中規定了 b 總是素數，是以我們可以很容易的找到一個b的原根broot，然後利用babystep求出broot^t1 mod b = c,然後我們能利用同餘原理輕易得出若t2=b-1，broot^(t1+n*t2) mod b = c ，那麼我們下一步需要求的就是 n*t2 + m*a = t1 其中的 n 和 m ，然後枚舉所有符合條件的m，求出broot^m1，broot^m2......broot^mi 一系列數字就是我們要求的答案了。

#include<iostream>

#include<cstdio>

#include<cstring>

#include<cmath>

#include<algorithm>

#include<ctime>

#define bignum __int64

using namespace std;

//x^k=b%m 已知k,b,m,求所有小于m的解x  m是素數

//大整數分解

//求a，b的最大公約數

bignum gcd(bignum a,bignum b)

{

    return b==0?a:gcd(b,a%b);

}

//求a*b%c,因為a，b很大，是以要先将b寫成二進制數，再加：例如3*7=3*(1+2+4);

bignum mulmod(bignum a,bignum b,bignum c)

{

    bignum cnt=0,temp=a;

    while(b)

    {

        if(b&1) cnt=(cnt+temp)%c;

        temp=(temp+temp)%c;

        b>>=1;

    }

    return cnt;

}

//求a^b%c,再次将b寫成二進制形式，例如：3^7=3^1*3^2*3^4;

bignum powmod(bignum a,bignum b,bignum c)

{

    bignum cnt=1,temp=a;

    while(b)

    {

        if(b&1) cnt=mulmod(cnt,temp,c);//cnt=(cnt*temp)%c;

        temp=mulmod(temp,temp,c);//temp=(temp*temp)%c;

        b>>=1;

    }

    return cnt;

}

//Miller-Rabin測試n是否為素數,1表示為素數，0表示非素數

int pri[10]={2,3,5,7,11,13,17,19,23,29};

bool Miller_Rabin(bignum n)

{

    if(n<2) return 0;

    if(n==2) return 1;

    if(!(n&1)) return 0;

    bignum k=0,m;

    m=n-1;

    while(m%2==0) m>>=1,k++;//n-1=m*2^k

    for(int i=0;i<10;i++)

    {

        if(pri[i]>=n) return 1;

        bignum a=powmod(pri[i],m,n);

        if(a==1) continue;

        int j;

        for(j=0;j<k;j++)

        {

            if(a==n-1) break;

            a=mulmod(a,a,n);

        }

        if(j<k) continue;

        return 0;

    }

    return 1;

}

//pollard_rho 大整數分解，給出n的一個非1因子，傳回n是為一次沒有找到

bignum pollard_rho(bignum C,bignum N)

{

    bignum I, X, Y, K, D;

    I = 1;

    X = rand() % N;

    Y = X;

    K = 2;

    do

    {

        I++;

        D = gcd(N + Y - X, N);

        if (D > 1 && D < N) return D;

        if (I == K) Y = X, K *= 2;

        X = (mulmod(X, X, N) + N - C) % N;

    }while (Y != X);

    return N;

}

//找出N的最小質因數

bignum rho(bignum N)

{

    if (Miller_Rabin(N)) return N;

    do

    {

        bignum T = pollard_rho(rand() % (N - 1) + 1, N);

        if (T < N)

        {

              bignum A, B;

              A = rho(T);

              B = rho(N / T);

              return A < B ? A : B;

        }

    }

    while(1);

}

//N分解質因數，這裡是所有質因數，有重複的

bignum AllFac[1100];

int Facnum;

void findrepeatfac(bignum n)

{

    if(Miller_Rabin(n))

    {

        AllFac[++Facnum]=n;

        return ;

    }

    bignum factor;

    do

    {

        factor=pollard_rho(rand() % (n - 1) + 1, n);

    }while(factor>=n);

    findrepeatfac(factor);

    findrepeatfac(n/factor);

}

//求N的所有質因數，是除去重複的

bignum Fac[1100];

int num[1100];

int len;//0-len

void FindFac(bignum n)

{

    len=0;

    //初始化

    memset(AllFac,0,sizeof(AllFac));

    memset(num,0,sizeof(num));

    Facnum=0;

    findrepeatfac(n);

    sort(AllFac+1,AllFac+1+Facnum);

    Fac[0]=AllFac[1];

    num[0]=1;

    for(int i=2;i<=Facnum;i++)

    {

        if(Fac[len]!=AllFac[i])

        {

            Fac[++len]=AllFac[i];//important

        }

        num[len]++;

    }

}

//求n的歐拉函數值

bignum oula(bignum n)

{

    FindFac(n);

    bignum cnt=n;

    for(int i=0;i<=len;i++)

    {

        cnt-=cnt/Fac[i];

    }

    return cnt;

}

//枚舉n的所有因子  cnt

bignum yinzi[100000];

bignum yinzinum;//初始化在main中(0-yinzinum-1)

void dfs(int id,bignum cnt)

{

    yinzi[yinzinum++]=cnt;

    if(id==len+1)

    {

        return ;

    }

    bignum temp=1;

    for(int i=0;i<=num[id];i++)

    {

        dfs(id+1,cnt*temp);

        temp*=Fac[id];

    }

}

//求原根 p是素數

int primitive(bignum p)

{

    //枚舉n的所有因子

    if(p==2) return 1;

    FindFac(p-1);

    yinzinum=0;//因子個數

    dfs(0,1);

    sort(yinzi,yinzi+yinzinum);

    yinzinum=unique(yinzi,yinzi+yinzinum)-yinzi;

    for(int i=2;i<p;i++)

    {

        int flag=1;

        for(int j=0;j<yinzinum-1;j++)//i^(p-1)=1(mod p)

        {

            if(powmod(i,yinzi[j],p)==1)

            {

                flag=0;break;

            }

        }

        if(flag) return i;

    }

    return -1;

}

//baby-step

struct Node

{

    int index;//序号

    __int64 val;//值

};

Node rnum[6600000];

__int64 _pow(__int64 a,__int64 b,__int64 m)

{

    if(b==0) return 1;

    if(b==1) return a%m;

    __int64 t=_pow(a,b/2,m);

    t=(t*t)%m;

    if(b&1) t=(t*a)%m;

    return t;

}

bool cmp(Node h,Node k)

{

    if(h.val!=k.val) return h.val<k.val;

    return h.index<k.index;

}

int _binary_search(__int64 key,__int64 len)//二分查找

{

    int left = 0,right = len;

    int result = - 1,mid;

    while(left <= right)

    {

        mid = (left + right)>>1;

        if(rnum[mid].val == key)

        {

            result = rnum[mid].index;

            right = mid -1;

        }

        else if(rnum[mid].val < key)  left = mid + 1;

        else   right = mid -1;

    }

return result;

}

int Baby_step_Giant_step(__int64 a,__int64 b,__int64 n)//計算最小的x使得a^x=b(mod n) n是素數

{

    __int64 m=(__int64)(sqrt(n+0.5)+1);

    //計算 (j,a^j%n) (0<=j<m)

    rnum[0].index=0,rnum[0].val=1;

    for(int i=1;i<m;i++) rnum[i].index=i,rnum[i].val=rnum[i-1].val*a%n;

    sort(rnum,rnum+m,cmp);

    //計算 a^-m

    __int64 a_m=_pow(a,n-1-m,n);//因為n是素數,是以a^(n-1)=1(mod n)

    __int64 y=b;

    for(int i=0;i<m;i++)

    {

        int res=_binary_search(y,m);

        if(res!=-1) return i*m+res;

        y=y*a_m%n;

    }

    return -1;

}

//exgcd

bignum exgcd(bignum a,bignum b,bignum&x,bignum &y)

{

    if(b==0) return x=1,y=0,a;

    bignum r=exgcd(b,a%b,x,y);

    bignum t=x;

    x=y;

    y=t-a/b*y;

    return r;

}

//求(a/b)%p,gcd(B,p)=1

__int64 calc(__int64 a,__int64 b,__int64 mod)

{

    __int64 x,y;

    __int64 r=exgcd(b,mod,x,y);

    x*=a;

    x=(x%mod+mod)%mod;

    return x;

}

//x^k=b(mod m)  m是素數

bignum ans[1000];

int ansnum;

//x^k=b(mod m)  m是素數  k>0,b>0,m>1

void solve(bignum k,bignum b,bignum m)

{

    bignum root=primitive(m);//存在

    bignum t1=Baby_step_Giant_step(root,b,m);

    bignum t2=m-1;

    //n*t2 + m*k = t1

    bignum x,y;

    bignum r=exgcd(t2,k,x,y);

    if(t1%r)

    {

        printf("-1\n");

        return ;

    }

    x*=t1/r,y*=t1/r;

    ansnum=0;

    while(y<0) x-=k/r,y+=t2/r;

    for(;;)

    {

        bignum t=powmod(root,y,m);

        int flag=1;

        for(int j=0;j<ansnum;j++)

        {

            if(ans[j]==t)

            {

                flag=0;

                break;

            }

        }

        if(!flag) break;

        ans[ansnum++]=t;

        x-=k/r,y+=t2/r;

    }

    sort(ans,ans+ansnum);

    for(int i=0;i<ansnum;i++) printf("%I64d\n",ans[i]);

}

int main()

{

    srand(time(NULL));

    int pl=1;

    __int64 x,k,b,m;

    while(scanf("%I64d%I64d%I64d",&k,&m,&b)==3)

    {

        printf("case%d:\n",pl++);

        solve(k,b,m);

    }

    return 0;

}