Hdu 2815 Mod Tree + Poj 3243 Clever Y 扩展Baby Step Giant Step 解决离散对数问题

两道题除了输入输出略有不同，其余基本一样

这类问题用哈希表效率比较高，本题加入哈希表是不判重回WA，渣渣我还没想明白为什么……

以下内容转自 POJ 3243 Clever Y（babyStepGiantStep） - 将狼踩尽 19891101 - 博客园

题意：给定A、B、P，求x满足A^x%P=B。

思路：首先，这里的(A,P)不一定等于1，即A与P不一定互质，不能直接利用babyStepGiantStep。下面我们进行变换。首先我们暴力x到100，若在100以内有解则直接返回。下面我们讨论的是x大于等于100的情况。

以下内容节选自：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 <cstdio>
#include <cstring>
#include <cmath>

#define i64 __int64

const int N=65535;
int head[N],e;

struct Node
{
    int id,val,next;
}edge[N];

void Insert (int val,int id)
{
	int k=val%N;
	for (int i=head[k];i!=-1;i=edge[i].next)
		if (edge[i].val==val)
			return ;   //不重复加边
	edge[e].id=id;
	edge[e].val=val;
	edge[e].next=head[k];
	head[k]=e++;
}

int Find (int val)
{
	int k=val%N;
	for (int i=head[k];i!=-1;i=edge[i].next)
		if (edge[i].val==val)
			return edge[i].id;
	return -1;
}

__int64 modPow (__int64 s,__int64 index,__int64 mod)  
//蒙哥马利幂模算法
//快速幂
//返回值(s^index)%mod
{  
    __int64 ans=1;  
    s%=mod;  
    while (index>=1)  
    {  
        if ((index&1)==1)   //奇数  
            ans=(ans*s)%mod;  
        index>>=1;  
        s=s*s%mod;  
    }  
    return ans;  
} 

__int64 gcd (__int64 a,__int64 b)
{
    return b==0?a:gcd(b,a%b);
}

i64 Extended_Euclid (i64 a,i64 b,i64 &x,i64 &y)
{//扩展欧几里得算法,求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
    int d;
    if (b==0)
    {
        x=1;y=0;
        return a;
    }
    d=Extended_Euclid(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

int Inval (int a,int b,int n)
{// 求解方程a*x+b*y=n;
    i64 e,x,y;  
    Extended_Euclid (a,b,x,y);
    e=x*n%b;
    return e<0?e+b:e;  
}

int exBSGS (int a,int b,int n)
{//A^x = B (mod n)  n无限制
	memset(head,-1,sizeof(head));
	memset(edge,0,sizeof(edge));
	e=0;
	i64 D=1,p=1;
	int d=0,i,w;
    for (i=0;i<=100;i++)
    {
        if (p==b) return i;
		p=(i64)p*a%n;
	}

	int tmp=gcd(a,n);
	while (tmp!=1)
	{
		if (b%tmp) return -1;
		d++;
		n/=tmp;
		b/=tmp;
		D=D*a/tmp%n;
		tmp=gcd(a,n);
	}
	int m=ceil(sqrt(1.0*n));
	for (i=0 , p=1;i<=m;i++, p = p*(i64)a%n )  //p对应a^i
		Insert(p,i);

	int mm=modPow(a,m,n); //a^m%C
	for (i=0;i<=m;i++, D=D*mm%n)   // mm^i
	{
		tmp=Inval(D,n,b);
		if (tmp>=0 && (w=Find(tmp))!=-1)
			return i*m+w+d;
	}
    return -1;
}
/*
int main ()  //Hdu2815
{
	int k,p,n;
	while (~scanf("%d%d%d",&k,&p,&n))
	{
		if (n>=p)
		{
			printf("Orz,I can’t find D!\n");
			continue;
		}
        if (p==1)
        {
            if (n==0) printf("0\n");
            else printf("Orz,I can’t find D!\n");
            continue;
        }
        int ans=exBSGS(k,n,p);
        if (ans==-1)
			printf("Orz,I can’t find D!\n");
        else
			printf("%d\n",ans);
    }
    return 0;
}*/

int main () //Poj 3243
{
#ifdef ONLINE_JUDGE
#else
	freopen("read.txt","r",stdin);
#endif
	int x,z,k;
	while (scanf("%d%d%d",&x,&z,&k) , k||z||k)
	{
		k%=z;
        int ans=exBSGS(x,k,z);
        if (ans==-1)
			printf("No Solution\n");
        else
			printf("%d\n",ans);
    }
    return 0;
}