hdu 2815 baby-step

转载请注明出处，谢谢   http://blog.csdn.net/ACM_cxlove?viewmode=contents            by---cxlove

对于A^X==B（MOD C）的情况，如果A与C不互质，那么普通的baby_step giant_step 便不能解决，AC给出了消因子的做法，并且

有证明 

【扩展Baby Step Giant Step解决离散对数问题】

原创帖！转载请注明作者 AekdyCoin ! http://hi.baidu.com/aekdycoin/item/236937318413c680c2cf29d4 【普通Baby Step Giant Step】 【问题模型】 求解 A^x = B (mod C) 中 0 <= x < C 的解,C 为素数 【思路】 我们可以做一个等价 x = i * m + j  ( 0 <= i < m, 0 <=j < m) m = Ceil ( sqrt( C) ) 而这么分解的目的无非是为了转化为: (A^i)^m * A^j = B ( mod C) 之后做少许暴力的工作就可以解决问题： (1) for i = 0 -> m, 插入Hash (i, A^i mod C) (2) 枚举 i ,对于每一个枚举到的i,令  AA = (A^m)^i mod C 我们有 AA * A^j = B (mod C) 显然AA,B,C均已知，而由于C为素数,那么 (AA,C)无条件为1 于是对于这个模方程解的个数唯一(可以利用扩展欧几里得或 欧拉定理来求解) 那么对于得到的唯一解X,在Hash表中寻找,如果找到，则返回  i * m + j  注意: 由于i从小到大的枚举,而Hash表中存在的j必然是对于某个剩余系内的元素X 是最小的(就是指标) 所以显然此时就可以得到最小解 如果需要得到 x > 0的解，那么只需要在上面的步骤中判断 当 i * m + j > 0 的时候才返回 到目前为止，以上的算法都不存在争议,大家实现的代码均相差不大。可见当C为素数的时候，此类离散对数的问题可以变得十分容易实现。 【扩展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<algorithm> #include<cmath> #define LL __int64 #define N 1000000 using namespace std; struct Node { LL idx; LL val; } baby[N]; bool cmp(Node n1,Node n2) { return n1.val!=n2.val?n1.val<n2.val:n1.idx<n2.idx; } LL gcd(LL a,LL b) { return b==0?a:gcd(b,a%b); } LL extend_gcd(LL a,LL b,LL &x,LL &y) { if(b==0) { x=1; y=0; return a; } LL gcd=extend_gcd(b,a%b,x,y); LL t=x; x=y; y=t-a/b*y; return gcd; } LL inval(LL a,LL b,LL n) { LL e,x,y; extend_gcd(a,n,x,y); e=((LL)x*b)%n; return e<0?e+n:e; } LL PowMod(LL a,LL b,LL MOD) { LL ret=1%MOD,t=a%MOD; while(b) { if(b&1) ret=((LL)ret*t)%MOD; t=((LL)t*t)%MOD; b>>=1; } return ret; } LL BinSearch(int num,int m) { int low=0,high=m-1,mid; while(low<=high) { mid=(low+high)>>1; if(baby[mid].val==num) return baby[mid].idx; if(baby[mid].val<num) low=mid+1; else high=mid-1; } return -1; } LL BabyStep(LL A,LL B,LL C) { LL tmp,D=1%C; LL temp; for(LL i=0,tmp=1%C; i<100; i++,tmp=((LL)tmp*A)%C) if(tmp==B) return i; LL d=0; while((temp=gcd(A,C))!=1) { if(B%temp) return -1; d++; C/=temp; B/=temp; D=((A/temp)*D)%C; } LL m=(LL)ceil(sqrt((long double)C)); for(LL i=0,tmp=1%C; i<=m; i++,tmp=((LL)tmp*A)%C) { baby[i].idx=i; baby[i].val=tmp; } sort(baby,baby+m+1,cmp); int cnt=1; for(int i=1; i<=m; i++) if(baby[i].val!=baby[cnt-1].val) baby[cnt++]=baby[i]; LL am=PowMod(A,m,C); for(LL i=0; i<=m; i++,D=((LL)(D*am))%C) { LL tmp=inval(D,B,C); if(tmp>=0) { LL pos=BinSearch(tmp,cnt); if(pos!=-1) return i*m+pos+d; } } return -1; } int main() { LL A,B,C; while(scanf("%lld%lld%lld",&A,&C,&B)!=EOF) { if(B>=C) { puts("Orz,I can’t find D!"); continue; } LL ans=BabyStep(A,B,C); if(ans==-1) puts("Orz,I can’t find D!"); else printf("%lld\n",ans); } return 0; }