HDU 2815 Mod Tree【扩展Baby Step Giant Step解决离散对数问题】

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

又是AC神，数论帝，太强了。

对于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{
	int idx;
	int val;
}baby[N];
bool cmp(Node n1,Node n2){
	return n1.val!=n2.val?n1.val<n2.val:n1.idx<n2.idx;
}
int gcd(int a,int b){
	return b==0?a:gcd(b,a%b);
}
int extend_gcd(int a,int b,int &x,int &y){
	if(b==0){
		x=1;
		y=0;
		return a;
	}
	int gcd=extend_gcd(b,a%b,x,y);
	int t=x;
	x=y;
	y=t-a/b*y;
	return gcd;
}
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;
}
int PowMod(int a,int b,int 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 (int)ret;
}
int 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;
}
int BabyStep(int A,int B,int C){
	LL tmp,D=1%C;
	int temp;
	for(int i=0,tmp=1%C;i<100;i++,tmp=((LL)tmp*A)%C)
		if(tmp==B)
			return i;
	int d=0;
	while((temp=gcd(A,C))!=1){
		if(B%temp) return -1;
		d++;
		C/=temp;
		B/=temp;
		D=((A/temp)*D)%C;
	}
	int m=(int)ceil(sqrt((double)C));
	for(int 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];
	int am=PowMod(A,m,C);
	for(int i=0;i<=m;i++,D=((LL)(D*am))%C){
		int tmp=inval(D,B,C);
		if(tmp>=0){
			int pos=BinSearch(tmp,cnt);
			if(pos!=-1)
				return i*m+pos+d;
		}
	}
	return -1;
}
int main(){
	int A,B,C;
	while(scanf("%d%d%d",&A,&C,&B)!=EOF){
		if(B>=C){
			puts("Orz,I can’t find D!");
			continue;
		}
		int ans=BabyStep(A,B,C);
		if(ans==-1)
			puts("Orz,I can’t find D!");
		else
			printf("%d\n",ans);
	}
	return 0;
}