离散对数之BSGS(扩展BSGS)与Pohlig-Hellman(知识点整理+板子总结)

心得

离散对数这里还是不怎么熟,

很早就学过,然而不怎么会套板子,

现在总算是把板子整理下来了

思路来源

https://blog.csdn.net/qq_34921856/article/details/79794335

https://blog.csdn.net/litble/article/details/73252928

知识点整理

1.BSGS

a^{x}\equiv b(mod\ p)的最小x值,p仅为素数

x=i*m+j,j<m,则需解a^{i*m}\equiv b*a^{-j}(mod\ p)

j从0到m-1预处理等式右边的值v,令map[v]=j存下标,

再左端枚举i从0到m-1,去查map表即可,找到最小的i满足j在表中,i*m+j即为答案

复杂度最低时,令m=\sqrt pO(\sqrt p)解决一个x

/*==================================================*\
| Baby-Step-Giant-Step 大步小步算法
| 求 a^x === b (mod p) 中的 最小 x值 -- 此处p仅为素数
| 实际运用中用自己的hash表代替map可防TLE
\*==================================================*/
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll a,b,p;
ll extgcd(ll a,ll b,ll &x,ll &y)
{
	ll d=a;
	if(b)d=extgcd(b,a%b,y,x),y-=(a/b)*x;
	else x=1,y=0;
	return d;
}
ll BSGS(ll a, ll b, ll p) 
{ 
    a %= p; b %= p;
    map<ll, ll> h;
    ll m = ceil(sqrt(p)), x, y, d, t = 1, v = 1;
    for(ll i = 0; i < m; ++i)
	{
        if(h.count(t)) h[t] = min(h[t], i);
        else h[t] = i;
        t = (t*a) % p;
    }
    for(ll i = 0; i < m; ++i) 
	{
        d = extgcd(v, p, x, y);//vx==1(mod p)即v在模p意义下逆元 
        x = (x* b/d % p + p) % (p);
        if(h.count(x)) return i*m + h[x];
        v = (v*t) % p;
    }
    return -1;
}
int main()
{
	scanf("%lld%lld%lld",&a,&b,&p);//a^x==b(mod p)
	printf("%lld\n",BSGS(a,b,p));
	return 0;
}

2.扩展BSGS

a^{x}\equiv b(mod\ p)的最小x值,p不为素数

a^{x-1}*a\equiv b(mod\ p)

考虑扩展欧几里得,在求ax\equiv b(mod\ p)时,如果a、p不互质则a、b、p同除gcd(a,p)

以上方法也同理,a^{x-1}*\frac{a}{d}{}\equiv \frac{b}{d}(mod\ \frac{p}{d})

\frac{p}{d}与a互质,直接解x-1,所得加1即答案

否则再提一个a出来,a^{x-2}*\frac{a}{d_{1}}*\frac{a}{d_{2}}\equiv \frac{b}{d_{1}*d_{2}}(mod\ \frac{p}{d_{1}*d_{2}})

重复该过程,直至a与\frac{p}{\prod_{i=1}^{k}d_{i}}互质,

去解a^{x-k}*\frac{a^{k}}{\prod_{i=1}^{k}d_{i}} \equiv \frac{b}{\prod_{i=1}^{k}d_{i}}(mod\ \frac{p}{\prod_{i=1}^{k}d_{i}})

注意解出来的是x-k,再加k即答案,

这里整理的板子,是手写哈希的,开大于预处理部分的空间的最小素数值

#include<iostream>
#include<cstdio>
#include<climits>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
#define ll long long
#define mod 99991
int head[mod],tot;
ll Next[mod<<1],id[mod<<1],has[mod<<1];
ll gcd(ll x,ll y)
{
	ll r=x%y;
	while(r){x=y;y=r;r=x%y;}
	return y;
}
void add(ll x,ll num)
{
	ll k=num%mod;
	id[++tot]=x;has[tot]=num;
	Next[tot]=head[k];head[k]=tot;
}
ll find(ll num)
{
	ll k=num%mod;
	for(ll i=head[k];i!=-1;i=Next[i])
	if(has[i]==num)return id[i];
	return -1;
}
ll exbsgs(ll a,ll b,ll p)
{
	a%=p;b%=p;
	if(b==1)return 0;
	ll tmp=1,d=1,cnt=0;
	while((tmp=gcd(a,p))!=1)
	{
		if(b%tmp)return -1;
		cnt++;b/=tmp;p/=tmp;d=d*(a/tmp)%p;
		if(b==d)return cnt;//注意这里
	}
	ll m=ceil(sqrt(p)),q=1,j;
	//一般m开为根号p 视具体询问而改变 
	//看预处理和询问的次数 如只询问一次则应根号p 
	for(ll i=0;i<m;i++)
	{add(i,(q*b)%p);q=(q*a)%p;}
	for(ll i=1;i<=m+1;i++)
	{
		d=(d*q)%p;j=find(d);
		if(j!=-1)return i*m-j+cnt;
	}
	return -1;
}
void init()
{
	memset(head,-1,sizeof head);
	tot=0;
}
int main()
{
	ll a,b,p,ans;
	while(~scanf("%lld%lld%lld",&a,&p,&b))
	{
		if(!a&&!b&&!p)break;
		init();
		ans=exbsgs(a,b,p);
		if(ans==-1)puts("No Solution");
		else printf("%lld\n",ans);
	}
	return 0;
}

3.Pohlig-Hellman

http://www-math.ucdenver.edu/~wcherowi/courses/m5410/phexam.html

啃了一天,贴一下网页局部,描述一下算法流程,讲点心得

a^{x}\equiv b(mod\ p)的最小x值,p仅为素数但p<=1e18,此时O(\sqrt p)解决不了x,

但p的特点在于,p-1可以表示为小质数的乘积,如例p-1=2^{2}*3^{4}*5^{2}

如果第一步能先求出x模每个素因子的幂p_{i}^{k_{i}}的值,届时第二步CRT合并即可

 

考虑第一步,

要求最小循环节,初始循环节cycle=p-1

应当检查每个素因子是否满足,对于当前循环节cycle而言,

a^{\frac{cycle}{p_{i}}}\equiv 1(mod\ p)是否成立,如果成立,

说明cycle可以除掉p_{i},使循环节更小,

那就除掉p_{i},直至该素因子p_{i}不满足a^{\frac{cycle}{p_{i}}}\equiv 1(mod\ p)

再去检查下一素因子,直至检查完所有素因子

 

经验证,这里的p的最小循环节为p-1,

先预处理a^{1*\frac {p-1}{p_{i}}},...a^{(p_{i}-1)*\frac {p-1}{p_{i}}},对于x2而言,只有a^{1*\frac{p-1}{2}}\equiv 8100(mod\ 8101)

然后去求(a^{x})^{\frac{p-1}{2}}\equiv7531^{\frac{p-1}{2}}\equiv8100(mod\ p)

x_{2}是x除以4的余数,可以表示为x_{2}=c_{0}+c_{1}*2

x=c_{0}+c_{1}*2+c_{2}*4+...,显然后面包含*2的项是循环节的倍数,被循环节消掉了,

(a^{x})^{\frac{p-1}{2}}\equiv(a^{c_{0}+c_{1}*2+c_{2}*4+...})^{\frac{p-1}{2}}\equiv a^{c_0*\frac{p-1}{2}},查表知c_{0}=1

对原余数7531乘以a^{c_{0}}的逆元,新的x=c_{1}*2+c_{2}*4+...,重复该过程,

(a^{x})^{\frac{p-1}{4}}\equiv(a^{c_{1}*2+c_{2}*4+...})^{\frac{p-1}{4}}\equiv a^{c_1*2*\frac{p-1}{4}}\equiv a^{c_1*\frac{p-1}{2}}\equiv 1,查表知c_{1}=0

由于x_{2}=c_{0}+c_{1}*2,代入解得x_{2}=1,该过程终止

 

对于每个素因子p_{i}^{k_{i}},都类似地解出其余数x_{i},得到同余方程组\begin{bmatrix} x \equiv x_{2}(mod\ 2^{2})\\ x \equiv x_{3}(mod\ 3^{4}) \\ x \equiv x_{5}(mod\ 5^{2}) \end{bmatrix}

由于素因子幂次之间显然互质,裸的CRT,合并构造一下即得x,

 

附:hdu6632 discrete logarithm problem

 求最小的x满足a^{x}\equiv b(mod\ p),p-1的素因子只由2、3组成

2<=a,b<=p-1,65537<=p<=10^{18}

据说是CTF界人尽皆知sb题,然而我一个打ACM的怎么可能会知道

(抓周树人和我鲁迅有什么关系)

仿着网页的过程敲,写了一波又丑常数又大的代码

#include<bits/stdc++.h>
using namespace std;
const int maxn=1e4+10;
const int N=60;
const int MAXP=20;//最大素因子的大小 
typedef long long ll;
bool ok[maxn],yes;
int prime[maxn],cnt;
int t;
int v[N],num[N],tot;//素数因子v[i] 出现了num[i]次 
ll X[N],P[N],sz;//返回(x',pi^ki) x'=x(mod pi^ki) 用于CRT 
ll a,b,p;
void sieve()
{
	for(ll i=2;i<maxn;++i)
	{
		if(!ok[i])prime[cnt++]=i;
		for(int j=0;j<cnt;++j)
		{
			if(i*prime[j]>=maxn)break;
			ok[i*prime[j]]=1;
			if(i%prime[j]==0)break; 
		}
	}
}
int fac(ll x)
{
	int tot=0;
	for(int i=0;i<cnt;++i)
	{
		if(1ll*prime[i]*prime[i]>x)break;
		if(x%prime[i]==0)
		{
			v[++tot]=prime[i];
			num[tot]=0;
			while(x%prime[i]==0)
			{
				num[tot]++;
				x/=prime[i];
			}
		}
	}
	if(x>1)
	{
		v[++tot]=x;
		num[tot]=1;
	}
	return tot;
}
ll extgcd(ll a,ll b,ll &x,ll &y)
{
	ll d=a;
	if(b)d=extgcd(b,a%b,y,x),y-=(a/b)*x;
	else x=1,y=0;
	return d;
}
ll mul(ll u,ll v,ll p)
{
	return (u*v-(ll)((long double)u*v/p)*p+p)%p;
}
ll modpow(ll x,ll n,ll mod)
{
	ll res=1;
	for(;n;n/=2,x=mul(x,x,mod)%mod)
	if(n&1)res=mul(res,x,mod)%mod;
	return res;
}
ll Pholig_Hellman(ll a,ll B,ll p)//返回(x',pi^ki) x'=x(mod pi^ki) 
{
	ll r=p-1;
	//找a的阶(最小)r p-1不一定是最小 类似判原根的过程 
	//最后r应该化为最小的阶
	tot=fac(r);
	for(int i=1;i<=tot;++i)
	while(r%v[i]==0&&modpow(a,r/v[i],p)==1)r/=v[i];
	tot=fac(r);//重新求素因子 
	ll x,y,inv,tmp[MAXP],sz=0;
	extgcd(a,p,x,y);
	inv=(x%p+p)%p;
	for(int i=1;i<=tot;++i)
	{
		int c=v[i],d=num[i];
		ll b=B,f=0,g=1;
		ll l=1,now=r,h=modpow(a,now/c,p);
		//printf("h:%lld\n",h);
		for(int k=0;k<c;++k,l=mul(l,h,p))
		{
			tmp[k]=l; 
			//printf("tmp[%d]:%lld\n",k,tmp[k]);
		}
		for(int j=1;j<=d;++j)
		{
			now/=c;//(p-1)/(c^j)
			ll e=modpow(b,now,p);
			ll co=-1,dif;//系数co 
			for(int k=0;k<c;++k)
			{
				if(e==tmp[k])
				{
					co=k; 
					break;
				}
			}
			if(co==-1)return -1;//无解
			//printf("co:%lld\n",co);
			co*=g;//实际系数 
			f+=co;
			dif=modpow(inv,co,p);
			b=mul(b,dif,p); 
			g*=c;
		}
		//printf("f:%lld g:%lld\n",f,g);
		X[++sz]=f;
		P[sz]=g;
	}
	return sz;
}
ll CRT(ll r[],ll m[],int n)//x'==x(mod p)方程组 共n项 要求pi互素 
{
	ll ans=0,M=1,x,y;
	for(int i=1;i<=n;++i)
	M*=m[i];
	for(int i=1;i<=n;++i)
	{
		ll mi=M/m[i];
		extgcd(mi,m[i],x,y);
		x=(x%m[i]+m[i])%m[i];
		ans=(ans+mul(mul(r[i],mi,M),x,M))%M;
	} 
	return (ans+M)%M;
}
int main()
{
	sieve();
	scanf("%d",&t);
	while(t--)
	{
		scanf("%lld%lld%lld",&p,&a,&b); 
		sz=Pholig_Hellman(a,b,p);
		if(sz!=-1)printf("%lld\n",CRT(X,P,sz));
		else puts("-1");
	}
	return 0;
}
//p a b
/*
105
65537 2 4
8101 6 7531
*/

 

相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页