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 Fac[1100];
int num[1100];
int len;//0-len
*/
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;
}