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

深藏阁楼爱情的钟 2022-09-17 10:25 11阅读 0赞

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

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

以下内容转自 [POJ 3243 Clever Y（babyStepGiantStep） - 将狼踩尽 19891101 - 博客园][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的情况。

![27123222-78bd7fbae111496ca9d4129c1e44cad8.jpg][]

以下内容节选自：[AekdyCoin的空间][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;
    }

[POJ 3243 Clever Y_babyStepGiantStep_ - _ 19891101 -]: http://www.cnblogs.com/jianglangcaijin/archive/2013/04/26/3046061.html
[27123222-78bd7fbae111496ca9d4129c1e44cad8.jpg]: /images/20220828/b9a9001799f64ed98cbacaabe07877f0.png
[AekdyCoin]: http://hi.baidu.com/aekdycoin/item/236937318413c680c2cf29d4