BZOJ 1420: Discrete Root【原根+BSGS+exgcd】

本文介绍了一种解决模意义下方程x^k=a(modp)的方法,通过使用原根和离散对数技巧,结合BSGS算法和扩展欧几里得算法,找到所有可能的根。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

题面:

Description
已知k,a,p,求x ^ k=a (mod p)的所有根(根的范围[0,p-1],P为质数
Input
三个整数p,k,a。0 < = a <p < = 10^9, 2 < = k < = 100000
Output
第一行一个整数,表示符合条件的x的个数。
第二行开始每行一个数,表示符合条件的x,按从小到大的顺序输出。

题目分析:


设p的原根为g,那么g的1到φ(p)\varphi(p)φ(p)次幂模p互不相等。
因为p是质数,所以任意一个数在模p意义下都可以用g的几次幂来表示。
(更一般的,一个与m互质的数都可以用m的原根的几次幂来表示)
I(x)I(x)I(x)表示x是g的几次幂,即:
gI(x)≡x (mod p)g^{I(x)}\equiv x~(mod~p)gI(x)x (mod p)gI(a)≡a (mod p)g^{I(a)}\equiv a~(mod~p)gI(a)a (mod p)
那么原方程就是 (gI(x))k≡gI(a) (mod p)(g^{I(x)})^k\equiv g^{I(a)}~(mod~p)(gI(x))kgI(a) (mod p)
就是要求k∗I(x)=I(a) (mod φ(p)=p−1)k*I(x)=I(a)~(mod~\varphi(p)=p-1)kI(x)=I(a) (mod φ(p)=p1)
用BSGS求出I(a)I(a)I(a),然后用exgcd求出所有满足条件的I(x)I(x)I(x)即可。

Upd:
重做了cf新年赛,然后发现我上面傻逼了
原方程就是(gI(x))k≡a (mod p)(g^{I(x)})^k\equiv a~(mod~p)(gI(x))ka (mod p)
换一换就是(gk)I(x)≡a (mod p)(g^{k})^{I(x)}\equiv a~(mod~p)(gk)I(x)a (mod p)
这不就妥妥的BSGS直接求出I(x)I(x)I(x)吗?!。。。

Upd2:
又做了一场模拟赛,发现一件有趣的事情:
k∗I(x)=I(a)(modp−1)k*I(x)=I(a)\pmod{p-1}kI(x)=I(a)(modp1),如果gcd(k,p−1)=1gcd(k,p-1)=1gcd(k,p1)=1,那么I(x)=I(a)∗k−1(modp−1)I(x)=I(a)*k^{-1}\pmod {p-1}I(x)=I(a)k1(modp1)
那么 x=ak−1(modp)x=a^{k^{-1}} \pmod px=ak1(modp)

数学题做起来真是有种异样的感觉。。

Code:

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<map>
#define maxn 100005
using namespace std;
int pr[maxn],cnt;
inline int exgcd(int a,int b,int &x,int &y){
	if(!b) {x=1,y=0;return a;}
	else {int d=exgcd(b,a%b,y,x);y-=a/b*x;return d;}
}
inline int ksm(int a,int b,int p){
	int s=1;
	for(;b;b>>=1,a=1ll*a*a%p) if(b&1) s=1ll*s*a%p;
	return s;
}
void getprime(int p){
	for(int i=2;i*i<=p;i++) if(p%i==0){
		pr[++cnt]=i;
		while(p%i==0) p/=i;
	}
	if(p>1) pr[++cnt]=p;
}
int getrg(int p){
	if(p==2) return 1;
	getprime(p-1);
	for(int i=2;i<p;i++){
		for(int j=1;j<=cnt;j++) if(ksm(i,(p-1)/pr[j],p)==1) goto he;
		return i; he:;
	}
}
int BSGS(int a,int b,int p){
	if(b==1) return 0;
	int m=int(sqrt(p)+1),base=b;
	map<int,int>has;
	for(int i=0;i<m;i++) has[base]=i,base=1ll*base*a%p;
	int now=1;base=ksm(a,m,p);
	for(int i=1;i<=m;i++)
		if(has.count(now=1ll*now*base%p)) return i*m-has[now];
	return -1;
}
int ans[maxn],tot;
int main()
{
	int k,a,p;
	scanf("%d%d%d",&p,&k,&a);
	if(a==0) return printf("1\n0"),0;
	int g=getrg(p);
	int A=BSGS(g,a,p),x,y;
	int d=exgcd(k,p-1,x,y);
	if(A==-1||A%d) return puts("0"),0;
	int sp=(p-1)/d; x=(1ll*x*A/d%sp+sp)%sp;
	if(!x) x=sp;//x->phi[p]=p-1
	for(int i=x;i<p;i+=sp) ans[++tot]=ksm(g,i,p);
	sort(ans+1,ans+1+tot);
	printf("%d\n",tot);
	for(int i=1;i<=tot;i++) printf("%d%c",ans[i],i==tot?10:32);
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值