HYSBZ 3283 运算器 exBSGS+exLucas

http://www.caiyiwen.tech/article/21.html

• 运算器
  • 这是一道模板题，通过此题即可获得两个重要模板——exBSGS和exLucas。
  • exBSGS和exLucas分别用以解决以下2、3问题：
  • 
  • 网上找到的模板并不稳定，经过反复调试最终确定模板如下代码。

代码：

#pragma GCC optimize(2)
#include<cstdio>
#include<cmath>
#include<cstring>
using namespace std;
namespace exLucas
{
	const int N=30,M=(int)1e6+10;
	int pri[M],p;
	char npri[M];
	long long aa[N],mm[N];
	void exgcd(long long a,long long b,long long&x,long long&y)
	{
		if(b==0)
		{
			x=1,y=0;
		}
		else
		{
			exgcd(b,a%b,y,x),y-=a/b*x;
		}
	}
	long long inv(long long x,long long mo)
	{
		long long nx, ny;
		exgcd(x, mo, nx, ny);
		return (nx%mo+mo)%mo;
	}
	long long qpow(long long x,long long y,long long mo)
	{
		long long ret=1;
		while(y)
		{
			if(y&1ll)
			{
				ret=ret*x%mo;
			}
			x=x*x%mo, y>>=1ll;
		}
		return ret;
	}
	long long f(long long x,long long mo,long long p)
	{
		if(x<=1)
		{
			return 1;
		}
		long long ret=1;
		for(int i=2;i<=p;i++)
		{
			if(i%mo)
			{
				ret=ret*i%p;
			}
		}
		ret=qpow(ret,x/p,p);
		for(int i=2;i<=x%p;i++)
		{
			if(i%mo)
			{
				ret=ret*i%p;
			}
		}
		return ret*f(x/mo,mo,p)%p;
	}
	inline void init()
	{
		npri[1]=1;
		for(int i=2;i<M;i++)
		{
			if(!npri[i])
			{
				pri[p++]=i;
			}
			for(int j=0;j<p&&i*pri[j]<M;j++)
			{
				npri[i*pri[j]]=1;
				if(i%pri[j]==0)
				{
					break;
				}
			}
		}
	}
	long long exlucas(long long n,long long m,long long mo)
	{
		int t=0;
		long long rmo=mo;
		for(int i=0;i<p&&mo>1;i++)
		{
			if(mo%pri[i]==0)
			{
				mm[++t]=1;
				while(mo%pri[i]==0)
				{
					mm[t]*=pri[i],mo/=pri[i];
				}
				aa[t]=f(n,pri[i],mm[t])*inv(f(m,pri[i],mm[t]),mm[t])%mm[t]*inv(f(n-m,pri[i],mm[t]),mm[t])%mm[t];
				int k=0;
				for(int j=n/pri[i];j;j/=pri[i])
				{
					k+=j;
				}
				for(int j=m/pri[i];j;j/=pri[i])
				{
					k-=j;
				}
				for(int j=(n-m)/pri[i];j;j/=pri[i])
				{
					k-=j;
				}
				aa[t]=aa[t]*qpow(pri[i],k,mm[t])%mm[t];
			}
		}
		mo=rmo;
		long long ret=0;
		for(int i=1;i<=t;i++)
		{
			(ret+=aa[i]*mo/mm[i]%mo*inv(mo/mm[i], mm[i]))%=mo;
		}
		return ret;
	}
}
namespace exBSGS
{
	const int NN=1e6, M=1e6+10;
	int h[NN], nexp[M], sp=1;
	long long a[M], b[M];
	inline void ins(int x,long long a1,long long b1)
	{
		nexp[sp]=h[x],h[x]=sp,a[sp]=a1,b[sp]=b1,sp++;
	}
	void clr()
	{
		memset(h,0,sizeof(h)),sp=1;
	}
	void exgcd(long long a,long long b,long long&x,long long&y)
	{
		if(b==0)
		{
			x=1,y=0;
		}
		else
		{
			exgcd(b,a%b,y,x),y-=a/b*x;
		}
	}
	long long inv(long long x,long long mo)
	{
		long long nx, ny;
		exgcd(x, mo, nx, ny);
		return (nx%mo+mo)%mo;
	}
	long long gcd(long long a,long long b)
	{
		if(a<b)
		{
			long long temp=a;
			a=b;
			b=temp;
		}
		return b?gcd(b,a%b):a;
	}
	long long exbsgs(long long y,long long z,long long p)
	{
		if(p==1||z==1)
		{
			return 0;
		}
		long long d=gcd(y,p);
		if(d>1)
		{
			if(z%d)
			{
				return -0x3f3f3f3f;
			}
			z/=d,p/=d,z=z*inv(y/d,p)%p;
			return exbsgs(y,z,p)+1;
		}
		int unit=ceil(sqrt(p));
		clr();
		long long yz=1;
		for(int i=0;i<unit;i++,yz=yz*y%p)
		{
			ins(yz*z%p%NN,yz*z%p,i);
		}
		long long nx=1;
		for(int i=0;i<=unit;i++,nx=nx*yz%p)
		{
			for(int u=h[nx%NN];u;u=nexp[u])
			{
				if(a[u]==nx&&i*unit-b[u]>=0)
				{
					return i*unit-b[u];
				}
			}
		}
		return -0x3f3f3f3f;
	}
}
int main()
{
	exLucas::init();
	int n,op,y,z,p;
	scanf("%d",&n);
	while(n--)
	{
		scanf("%d%d%d%d",&op,&y,&z,&p);
		if(op==1)
		{
			printf("%lld\n",exLucas::qpow(y,z,p));
		}
		else if(op==2)
		{
			long long tmp=exBSGS::exbsgs(y,z,p);
			if(tmp<0)
			{
				printf("Math Error\n");
			}
			else
			{
				printf("%lld\n",tmp);
			}
		}
		else
		{
			printf("%lld\n",exLucas::exlucas(z,y,p));
		}
	}
	return 0;
}