HDU 2815 Mod Tree 离散对数 扩张Baby Step Giant Step算法

联系：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\n",ans); } return 0; }