BZOJ2219数论之神——BSGS+中国剩余定理+原根与指标+欧拉定理+exgcd

题目描述

在ACM_DIY群中,有一位叫做“傻崽”的同学由于在数论方面造诣很高,被称为数轮之神!对于任何数论问题,他都能瞬间秒杀!一天他在群里面问了一个神题: 对于给定的3个非负整数 A,B,K 求出满足 (1) X^A = B(mod 2*K + 1) (2) X 在范围[0, 2K] 内的X的个数!自然数论之神是可以瞬间秒杀此题的,那么你呢?

输入

第一行有一个正整数T,表示接下来的数据的组数( T <= 1000) 之后对于每组数据,给出了3个整数A,B,K (1 <= A, B <= 10^9, 1 <= K <= 5 * 10^8)

输出

输出一行,表示答案

样例输入

3
213 46290770 80175784
3 46290770 80175784
3333 46290770 80175784

样例输出

27
27
297
 

这题果然数论神题啊,写了三遍才A掉。做这题之前建议先做一下弱化版BZOJ1319
模数是$2*K+1$(我们设$P=2*K+1$),显然不一定是质数。一般碰到模数不是质数的情况我们会想能否将模数质因数分解成互质的几个数然后通过解决子问题来获得问题答案,这道题也是可以的,对于$P$将它质因数分解为$p_{1}^{a_{1}}*p_{2}^{a_{2}}*...*p_{k}^{a_{k}}$,然后分别解决$x^A\equiv B(mod\ p_{i}^{a_{i}})$的解的个数,再将每个子问题解个数相乘即可得到原问题解的个数。为什么?因为$x^A\equiv B(mod P)$即$x^A-n*P=B$,那么$x^A-n_{i}*p_{i}^{a_{i}}=B$(即将P的其他部分移到$n_{i}$中)。那么$x_{i}^A\equiv B(mod\ p_{i}^{a_{i}})$(不考虑$p_{i}^{a_{i}}$与$B$的大小关系),因为在原式中$x\in [0,P)$,而新式子中$x_{i}\in [0,p_{i}^{a_{i}})$,所以$x\equiv x_{i}(mod\ p_{i}^{a_{i}})$,对于这$k$个子问题每个求出一个解$x_{i}$,用中国剩余定理都能求出唯一的原式解$x$。
现在考虑如何解决$x_{i}^A \equiv B(mod\ p_{i}^{a_{i}})$,先考虑$B$与$p_{i}^{a_{i}}$的关系:
1、$gcd(B,p_{i}^{a_{i}})=p_{i}^{a_{i}}$,即$B$是$p_{i}^{a_{i}}$的倍数,那么$x_{i}^A$中至少要有$p_{i}^{a_{i}}$,$x_{i}$中至少要有$p_{i}^{\left \lceil \frac{a_{i}}{A} \right \rceil}$,所以$x_{i}$要是$p_{i}^{\left \lceil \frac{a_{i}}{A} \right \rceil}$的倍数,在$[0,p_{i}^{a_{i}})$中这样的数个数为$p_{i}^{a_{i}-\left \lceil \frac{a_{i}}{A} \right \rceil}$。
2、$gcd(B,p_{i}^{a_{i}})=1$,我们求出$p_{i}^{a_{i}}$的原根$g$(求原根的方法是求出$\phi(p_{i}^{a_{i}})$的所有质因子$q_{1}...q_{k}$,然后从$2$开始枚举原根,如果一个数$G$满足$\forall G^{\frac{\phi(p_{i}^{a_{i}})}{q_{j}}}mod\ p_{i}^{a_{i}}\neq1$即为原根)。将$B$和$x_{i}$分别求出指标$indB$和$indx_{i}$,原式就变成了$A*indx_{i}\equiv indB(mod\ \phi(p_{i}^{a_{i}}))$,那么根据扩展欧几里得可知,若$indB\%gcd(A,\phi(p_{i}^{a_{i}})=0$,那么解的个数就是$gcd(A,\phi(p_{i}^{a_{i}})$。(解一下不定方程即可得到此结论)
3、$gcd(B,p_{i}^{a_{i}})>1$,设$B=m*p_{i}^{z}$,显然$x^A$中$p_{i}$的幂次只能是$z$,那么$A$必须是$z$的约数。我们将$x_{i}^A,B,p_{i}^{a_{i}}$都除掉$p_{i}^{z}$这一部分,问题就转化成第二种情况了。但并不是这样就完事了,我们发现原式中$x_{i}\in[0,p_{i}^{a_{i}}),\frac{x_{i}}{p_{i}^{\frac{z}{A}}}\in[0,p_{i}^{a_{i}-\frac{z}{A}})$,而除完的式子中的$x_{i}$(即原式中的$\frac{x_{i}}{p_{i}^{\frac{z}{A}}}$)$\in[0,p_{i}^{a_{i}-z})$,值域减小了,对于同余方程我们将值域扩大$k$倍,解的个数就会扩大$k$倍,所以最后还要乘上$p_{i}^{z-\frac{z}{A}}$。

#include<set>
#include<map>
#include<queue>
#include<stack>
#include<cmath>
#include<cstdio>
#include<vector>
#include<bitset>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll long long 
#define INF 2e9
using namespace std;
ll A,B,K;
ll g,f,p;
ll phi;
int T;
int prime[100010];
int cnt;
map<ll,int>mp;
ll quick(ll x,ll y,ll mod)
{
	ll res=1ll;
	while(y)
	{
		if(y&1)
		{
			res=res*x%mod;
		}
		x=x*x%mod;
		y>>=1;
	}
	return res;
}
ll gcd(ll x,ll y)
{
	return y==0?x:gcd(y,x%y);
}
ll get_ori(ll mod,ll phi)
{
	ll n=phi;
	cnt=0;
	for(int i=2;1ll*i*i<=phi;i++)
	{
		if(phi%i==0)
		{
			prime[++cnt]=i;
			while(phi%i==0)
			{
				phi/=i;
			}
		}
	}
	if(phi!=1)
	{
		prime[++cnt]=phi;
	}
	for(int i=1;i<=n;i++)
	{
		bool flag=true;
		for(int j=1;j<=cnt;j++)
		{
			if(quick(i,n/prime[j],mod)==1)
			{
				flag=false;
				break;
			}
		}
		if(flag)
		{
			return i;
		}
	}
}
ll BSGS(ll g,ll mod,ll x,ll phi)
{
	ll n=ceil(sqrt(mod));
	mp.clear();
	ll sum=1ll;
	for(int i=1;i<=n;i++)
	{
		sum*=g,sum%=mod;
		mp[sum]=i;
	}
	ll num=1ll;
	for(int i=0;i<=n;i++)
	{
		ll inv=quick(num,phi-1,mod)*x%mod;
		if(mp[inv])
		{
			return i*n+mp[inv];
		}
		num*=sum,num%=mod;
	}
}
ll solve(ll p,ll k)
{
	ll mod=quick(p,k,INF);
	if(B%mod==0)
	{
		return quick(p,k-(k-1)/A-1,INF);
	}
	ll b=B;
	ll num=0;
	while(b%p==0)
	{
		num++;
		b/=p;
		mod/=p;
		k--;
	}
	if(num%A)
	{
		return 0;
	}
	ll phi=mod-(mod/p);
	g=get_ori(mod,phi);
	f=BSGS(g,mod,b,phi);
	ll d=gcd(A,phi);
	if(f%d)
	{
		return 0;
	}
	return d*quick(p,num-num/A,INF);
}
int main()
{
	scanf("%d",&T);
	while(T--)
	{
		scanf("%lld%lld%lld",&A,&B,&K);
		ll ans=1ll;
		p=2*K+1;
		for(int i=2;1ll*i*i<=p;i++)
		{
			if(p%i==0)
			{
				int num=0;
				while(p%i==0)
				{
					num++;
					p/=i;
				}
				ans*=solve(i,num);
			}
		}
		if(p!=1)
		{
			ans*=solve(p,1);
		}
		printf("%lld\n",ans);
	}
}
posted @ 2019-02-15 09:42  The_Virtuoso  阅读(431)  评论(0编辑  收藏  举报