zoukankan      html  css  js  c++  java
  • 数论之BSGS算法

    UPDATE ON 2020.12.23 更新了原代码中的错误

    BSGS算法

    (BSGS(baby-step giant-step)),即大步小步算法。常用于求解离散对数问题。

    形式化地说,该算法可以在 (O sqrt{p}) 的时间内求解

    [a^x equiv b (mod p) ]

    其中 (gcd(a,p)=1)

    根据欧拉定理,当 (gcd(a,p)=1) 时,(a^{varphi(p)}mod p=1)

    (a^0mod p=1)

    所以从 (1 sim varphi(p)) 出现了一个循环节

    因此,我们只需要求出 (a^1 sim a^{varphi(p)})

    然后分别带入左右两边判断是否相等

    一般的题目中 (p) 为质数,所以 (varphi(p)=p-1)

    一般就用 (p) 替代

    但是这样做时间复杂度太高,可以用分块的思想优化

    (a^x=a^{km-t}),那么原式就变成了

    [a^{km} equiv ba^t (mod p) ]

    我们只需要分别求出 (t=1 sim m)(ba^t) 的值,将其存到哈希表里

    然后从 (1 sim p/m) 枚举 (k) ,判断哈希表中有没有和 (a^{km}) 相等的数

    (m=sqrt{p}) 时最优

    代码

    模板题

    #include<cstdio>
    #include<algorithm>
    #include<cmath>
    #include<cstring>
    #define rg register
    inline int read(){
       rg int x=0,fh=1;
       rg char ch=getchar();
       while(ch<'0' || ch>'9'){
       	if(ch=='-') fh=-1;
       	ch=getchar();
       }
       while(ch>='0' && ch<='9'){
       	x=(x<<1)+(x<<3)+(ch^48);
       	ch=getchar();
       }
       return x*fh;
    }
    const int maxn=1e5+5,mod=1e5+3;
    int p,a,n,blo,h[maxn],tot=1;
    struct asd{
       int nxt,val,num;
    }b[maxn];
    void ad(int val,int id){
       rg int now=val%mod;
       for(rg int i=h[now];i!=-1;i=b[i].nxt){
       	rg int cs=b[i].val;
       	if(cs==val){
       		if(b[i].num<id) b[i].num=id;
       		return;
       	}
       }
       b[tot].val=val;
       b[tot].nxt=h[now];
       b[tot].num=id;
       h[now]=tot++;
    }
    int cx(int val){
       rg int now=val%mod;
       for(rg int i=h[now];i!=-1;i=b[i].nxt){
       	rg int cs=b[i].val;
       	if(cs==val) return b[i].num;
       }
       return -1;
    }
    int main(){
       memset(h,-1,sizeof(h));
       p=read(),a=read(),n=read();
       blo=sqrt(p)+1;//不要忘了加上一个1保证正确性
       rg int now=n,cs=1,haha;
       for(rg int i=1;i<=blo;i++){
       	now=1LL*now*a%p;
       	ad(now,i);
       	cs=1LL*cs*a%p;
       }
       now=1;
       for(rg int i=1;i<=blo;i++){
       	now=1LL*now*cs%p;
       	haha=cx(now);
       	if(haha!=-1){
       		printf("%d
    ",i*blo-haha);
       		return 0;
       	}
       }
       printf("no solution
    ");
       return 0;
    }
    

    扩展BSGS算法

    求解

    [a^x equiv b (mod p) ]

    其中 (a,p) 不一定互质

    这个时候就不能再找 (1 sim varphi(p)) 的循环节了,因为不满足欧拉定理的条件

    我们把式子变一下,可以写成

    [a imes a^{x-1}+kp=b ]

    的形式

    根据裴蜀定理,当 (gcd(a,p)|b) 时方程有解,否则无解

    所以在有解的情况下,我们可以把方程左右两边同时除以 (gcd(a,p)),等号不变

    此时原方程变为

    [frac{a}{gcd(a,p)} imes a^{x-1}+k frac{p}{gcd(a,p)}=frac{b}{gcd(a,p)} ]

    (frac{p}{gcd(a,p)}) 看成新的 (p)

    (frac{b}{gcd(a,p)}) 看成新的 (b)

    同时记录一下左边的系数之积 (frac{a}{gcd(a,p)}),并不断地提出 (a)

    这样在最后一定会出现 (gcd(a,p)=1) 的情况,用普通的 (bsgs) 即可

    不要忘了把提出去的 (a) 的系数加上

    代码

    模板题

    #include<cstdio>
    #include<algorithm>
    #include<cmath>
    #include<cstring>
    #define rg register
    inline int read(){
    	rg int x=0,fh=1;
    	rg char ch=getchar();
    	while(ch<'0' || ch>'9'){
    		if(ch=='-') fh=-1;
    		ch=getchar();
    	}
    	while(ch>='0' && ch<='9'){
    		x=(x<<1)+(x<<3)+(ch^48);
    		ch=getchar();
    	}
    	return x*fh;
    }
    const int maxn=1e5+5,mod=1e5+3;
    int h[maxn],tot=1;
    struct asd{
    	int nxt,val,num;
    }b[maxn];
    void ad(int val,int id){
    	rg int now=val%mod;
    	for(rg int i=h[now];i!=-1;i=b[i].nxt){
    		rg int cs=b[i].val;
    		if(cs==val){
    			if(b[i].num<id) b[i].num=id;
    			return;
    		}
    	}
    	b[tot].val=val;
    	b[tot].nxt=h[now];
    	b[tot].num=id;
    	h[now]=tot++;
    }
    int cx(int val){
    	rg int now=val%mod;
    	for(rg int i=h[now];i!=-1;i=b[i].nxt){
    		rg int cs=b[i].val;
    		if(cs==val) return b[i].num;
    	}
    	return -1;
    }
    int bsgs(rg int a,rg int p,rg int n,rg int xs){
    	memset(h,-1,sizeof(h));
    	tot=1;
    	rg int blo=sqrt(p)+1;
    	rg int ac1=n,ac2=1,ac3=xs;
    	for(rg int i=1;i<=blo;i++){
    		ac1=1LL*ac1*a%p;
    		ad(ac1,i);
    		ac2=1LL*ac2*a%p;
    	}
    	for(rg int i=1;i<=blo;i++){
    		ac3=1LL*ac3*ac2%p;
    		ac1=cx(ac3);
    		if(ac1!=-1) return i*blo-ac1;
    	}
    	return -1;
    }
    int gcd(rg int aa,rg int bb){
    	return bb==0?aa:gcd(bb,aa%bb);
    }
    int main(){
    	rg int a,p,n,haha=1,now,ans,cnt;
    	rg bool jud=0;
    	while(1){
    		a=read(),p=read(),n=read();
    		if(a==0 && p==0 && n==0) break;
    		haha=1,jud=0,cnt=0;
    		a%=p,n%=p;
    		if(n==1){
    			printf("0
    ");
    		} else if(a==0 && n==0){
    			printf("0
    ");
    		} else if(a==0){
    			printf("No Solution
    ");
    		} else if(n==0){
    			while(1){
    				if(p==1){
    					printf("%d
    ",cnt);
    					break;
    				}
    				now=gcd(a,p);
    				if(now==1 && p!=-1){
    					printf("No Solution
    ");
    					break;
    				}
    				cnt++,p/=now;
    			}
    		} else {
    			while(1){
    				now=gcd(a,p);
    				if(n%now!=0){
    					jud=1;
    					break;
    				}
    				if(now==1){
    					ans=bsgs(a,p,n,haha);
    					if(ans==-1) jud=1;
    					ans+=cnt;
    					break;
    				}
    				cnt++,p/=now,n/=now;
    				haha=1LL*haha*(a/now)%p;
    				if(haha==n){
    					ans=cnt;
    					break;
    				}
    			}
    			if(jud){
    				printf("No Solution
    ");
    			} else {
    				printf("%d
    ",ans);
    			}
    		}
    	}
    	return 0;
    }
    
  • 相关阅读:
    JS清除IE浏览器缓存的方法
    大数据基础2
    CI/CD
    手机连接fiddler
    npm run build 报错
    django.core.exceptions.ImproperlyConfigured: mysqlclient 1.3.13 or newer is required; you have 0.9.2
    读取ini文件的方法
    ES小知识
    svn连接pycharm
    创建python文件时添加相关信息
  • 原文地址:https://www.cnblogs.com/liuchanglc/p/14140291.html
Copyright © 2011-2022 走看看