联系：http://acm.hdu.edu.cn/showproblem.php?pid=2815
意甲冠军：
思路：与上题不同。这道题不要求m是素数。是利用扩展Baby Step Giant Step算法求离散对数。
下面转载自：AekdyCoin
【扩展Baby Step Giant Step】
【问题模型】
求解
A^x = B (mod C) 中 0 <= x < C 的解,C 无限制(当然大小有限制……)
【写在前面】
这个问题比較麻烦。眼下网络上流传很多版本号的做法，只是大部分已近被证明是全然错误的！
这里就不再累述这些做法。以下是我的做法(有问题欢迎提出)
以下先给出算法框架。稍后给出具体证明:
(0) for i = 0 -> 50 if(A^i mod C == B) return i O(50)
(1) d<- 0 D<- 1 mod C
while((tmp=gcd(A,C))!=1)
{
if(B%tmp)return -1; // 无解!
++d;
C/=tmp;
B/=tmp;
D=D*A/tmp%C;
}
(2) m = Ceil ( sqrt(C) ) //Ceil是必要的 O(1)
(3) for i = 0 -> m 插入Hash表(i, A^i mod C) O( m)
(4) K=pow_mod(A,m,C)
for i = 0 -> m
解 D * X = B (mod C) 的唯一解 (假设存在解。必定唯一!)
之后Hash表中查询,若查到(如果是 j)，则 return i * m + j + d
否则
D=D*K%C,继续循环
(5) 无条件返回 -1 ;//无解!
以下是证明:
推论1:
A^x = B(mod C)
等价为
A^a * A^b = B ( mod C) (a+b) == x a,b >= 0
证明:
A^x = K * C + B (模的定义)
A^a * A^b = K*C + B( a,b >=0, a + b == x)
所以有
A^a * A^b = B(mod C)
推论 2:
令 AA * A^b = B(mod C)
那么解存在的必要条件为: 能够得到至少一个可行解 A^b = X (mod C)
使上式成立
推论3
AA * A^b = B(mod C)
中解的个数为 (AA,C)
由推论3 不难想到对原始Baby Step Giant Step的改进
For I = 0 -> m
For any solution that AA * X = B (mod C)
If find X
Return I * m + j
而依据推论3,以上算法的复杂度实际在 (AA,C)非常大的时候会退化到差点儿O(C)
归结原因，是由于(AA,C)过大，而就是(A,C)过大
于是我们须要找到一中做法，能够将(A,C)降低。并不影响解
以下介绍一种“消因子”的做法
一開始D = 1 mod C
进行若干论的消因子,对于每次消因子
令 G = (A,C[i]) // C[i]表示经过i轮消因子以后的C的值
假设不存在 G | B[i] //B[i]表示经过i轮消因子以后的B的值
直接返回无解
否则
B[i+1] = B[i] / G
C[i+1] = C[i] / G
D = D * A / G
详细实现仅仅须要用若干变量,细节參考代码
如果我们消了a'轮(如果最后得到的B,C分别为B',C')
那么有
D * A^b = B' (mod C')
于是能够得到算法
for i = 0 -> m
解 ( D* (A^m) ^i ) * X = B'(mod C')
因为 ( D* (A^m) ^i , C') = 1 (想想为什么？)
于是我们能够得到唯一解
之后的做法就是对于这个唯一解在Hash中查找
这样我们能够得到b的值,那么最小解就是a' + b !!
如今问题大约已近攻克了，但是细心看来，事实上还是有BUG的。那就是
对于
A^x = B(mod C)
假设x的最小解< a',那么会出错
而考虑到每次消因子最小消 2
故a'最大值为log(C)
于是我们能够暴力枚举0->log(C)的解，若得到了一个解，直接返回
否则必定有 解x > log(C)
PS.以上算法基于Hash 表,假设使用map等平衡树维护，那么复杂度会更大
代码：
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <cstdlib>
#include <queue>
#include <stack>
#include <vector>
#include <ctype.h>
#include <algorithm>
#include <string>
#include <set>
#define PI acos(-1.0)
#define maxn 1000005
#define INF 0x7fffffff
#define eps 1e-8
typedef long long LL;
typedef unsigned long long ULL;
using namespace std;
int extend_gcd(int a, int b, int &x, int &y)
{
if(b==0)
{
x=1;
y=0;
return a;
}
int r=extend_gcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
return r;
}
LL pow_mod(LL aa,LL ii,LL nn)
{
if(ii==0)
return 1%nn;
LL temp=pow_mod(aa,ii>>1,nn);
temp=temp*temp%nn;
if(ii&1)
temp=temp*aa%nn;
return temp;
}
struct b_step
{
int i,m;
} bb[maxn];
int inval(int a,int b,int n)
{
int e,x,y;
extend_gcd(a,n,x,y);
e=((LL)x*b)%n;
return e<0?e+n:e;
}
bool cmp(b_step a,b_step b)
{
return a.m==b.m?a.i<b.i:a.m<b.m;
}
int BiSearch(int m,LL num)
{
int low=0,high=m-1,mid;
while(low<=high)
{
mid=(low+high)>>1;
if(bb[mid].m==num)
return bb[mid].i;
if(bb[mid].m<num)
low=mid+1;
else
high=mid-1;
}
return -1;
}
int giant_step_baby_step(LL b,LL n,LL p)
{
LL tt=1%p;
for(int i=0; i<100; i++)
{
if(tt%p==n)
return i;
tt=((LL)tt*b%p);
}
LL D=1%p;
int d=0,temp;
while((temp=__gcd(b,p))!=1)
{
if(n%temp)
return -1;
d++;
n/=temp;
p/=temp;
D=((b/temp)*D)%p;
}
int m=(int)ceil(sqrt((double)(p)));
bb[0].i=0,bb[0].m=1%p;
for(int i=1; i<=m; i++)
{
bb[i].i=i;
bb[i].m=bb[i-1].m*b%p;
}
sort(bb,bb+m+1,cmp);
int top=1;
for(int i=1; i<=m; i++)
if(bb[i].m!=bb[top-1].m)
bb[top++]=bb[i];
int bm=pow_mod(b,m,p);
for(int i=0; i<=m; i++)
{
int tmp=inval(D,n,p);
if(tmp>=0)
{
int pos=BiSearch(top,tmp);
if(pos!=-1)
return i*m+pos+d;
}
D=((LL)(D*bm))%p;
}
return -1;
}
int main()
{
int b,p,n;
while(~scanf("%d%d%d",&b,&p,&n))
{
if(n>=p)
{
puts("Orz,I can’t find D!");
continue;
}
int ans=giant_step_baby_step(b,n,p);
if(ans==-1)
puts("Orz,I can’t find D!");
else printf("%d
",ans);
}
return 0;
}
版权声明：本文博主原创文章。博客，未经同意不得转载。