basic.hpp

求解 形如 a*x^2+b*x+c=0 (mod p)的二次同余方程

//定义一些同余方程求解计算的基本函数

// 处理类型为Type的数据,其应该定义+,-,*,/,mod,power,&,|,~运算


#ifndef BASIC				//避免重复包含
#define BASIC
#define INF -1				//错误标志

#include <typeinfo.h>
#include "verylong.cpp"
// -------------------------------------------------------- //函数声明 

template <class Type>  Type mod(const Type &,const Type &);					//模运算 
template <class Type>  Type gcd(const Type &,const Type &);					//最大公约数,史泰因Stein算法
template <class Type>  
Type power(const Type &x,const Type &y,const Type &p);	//计算x^y (mod p)
template <class Type>  Type power(const Type &x,const Type &y);			//计算x^y
template <class Type>  Type ninv(const Type &a,const Type &n);			//欧几里德算法,求同余逆
template <class Type>
Type Euclid(const Type &aa,const Type &nn,Type *gcd); //欧几里德算法

// -------------------------------------------------------- //


template <class Type>	//模运算 
Type mod(const Type &a,const Type &b)	//模运算  
{
	if(b<=0)
		return INF;
	if(a>=0)
		return a%b;
	else
		return a%b+b;
}

template <class Type>  	//最大公约数,扩展欧几里德算法
Type gcd(const Type &a,const Type &b) //定义&运算
{
	Type r,*gcd=new Type;
	Euclid(a,b,gcd);
	r=*gcd;
	delete gcd;
	return r;
}

template <class Type>		//求 x^y mod p,p应为正数
Type power(const Type &x,const Type &yy,const Type &pp)
{
	Type fw;  //fw 求值
	long lengh=sizeof(Type)*8;
	Type xx=x;
    if(xx<0)
		xx+=pp;

	if(pp<=0)
	{
		cout<<"参数错误!";
		fw=INF;
	}
	else
	{
		if(yy==0)
			fw=1;
		else   
		{
			fw=1;
			for(long i=lengh-1;i>=0;i--)
			{
				fw=mod(fw*fw,pp);		//反复平方乘
				if(yy & (1<<i))			//依次取yy的每一位
					fw=mod(fw*xx,pp);	
			}
		}
	}
	return fw;
}

template <class Type>
Type power(const Type &xx,const Type &yy)
{
	Type fw;  //fw 求值
	long lengh=sizeof(Type)*8;

	if(yy==0)
		fw=1;
	else   
	{
		fw=1;
		for(long i=lengh-1;i>=0;i--)
		{
			fw=fw*fw;			//反复平方乘
			if(yy & (1<<i))		//依次取yy的每一位
				fw=fw*xx;	
		}
	}
	return fw;
}

//利用辗转相除法求乘法逆
// 求aa模nn的乘法逆
template <class Type> 
Type ninv(const Type &a,const Type &nn)
{
	Type fw;
	Type aa=mod(a,nn);
	Type *g=new Type;
	fw=Euclid(aa,nn,g);		//扩展欧几里得算法求逆及最大公约数 gcd(a,n) fw=a^(-1) mod n
	delete g;
	return fw;
}


//辗转相除法求乘法逆或最大公约数
template <class Type>
Type Euclid(const Type &a,const Type &nn,Type *gcd)
{
	Type fw,aa=a;
	if(nn<=0 || aa==0)
	{
		if(aa==0&&nn>0)
			*gcd=nn;
		fw=INF;
		*gcd=INF;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 
	}
	else
	{
	   if(aa<0)
		   aa=nn-aa;

	   Type rs2=aa,rs1=nn,rr=1;
	   Type us2=1,us1=0;
	   Type vs2=0,vs1=1;
	   Type qq,uu,vv,ww;
	   while(rr!=0)
	   {
		   qq=(Type)(rs2/rs1);
		   rr=rs2-qq*rs1;
		   if(rr!=0)
		   {	  
			   uu=us2-qq*us1;
			   vv=vs2-qq*vs1;
			   rs2=rs1;rs1=rr;
			   us2=us1;us1=uu;
			   vs2=vs1;vs1=vv;
		   }
     
	   }
	   ww=rs1;
	   if(ww==1)
		  fw=mod(uu,nn);
	   else
	   {
           fw=INF;
	   }
	   *gcd=ww;
	}	
	return fw;
}

//多项式的带余除法计算 yb=yq*ya+yr,系数mod p
//mm为yb的次数,nn为ya的次数
template <class Type>
Type polydiv(Type *yb,Type *ya,Type *yq,Type *yr,long mm,long nn,const Type &pp)
{
	if(mm<nn)
	{
		for(long ii=0;ii<=mm;ii++)
			yr[ii]=yb[ii];
	   
	    for(ii=mm+1;ii<=nn-1;ii++)
		    yr[ii]=0;
	}
	else 
	{ 
		long kk=mm-nn;
		Type v0=ninv(ya[nn],pp);
		yq[kk]=mod(v0*yb[mm],pp);

		Type v1,v2;
		long j0,j1;
	    for(long tt=1;tt<=kk;tt++)
		{ 
			j0=0;
			v1=0;
		    if(nn>tt)
			   j0=nn-tt;  //保证yb的下标不越界
		   
		    for(long jj=j0;jj<=nn-1;jj++)
			   v1=v1+ya[jj]*yq[mm-tt-jj];
		   
		    yq[kk-tt]=mod(v0*(yb[mm-tt]-v1),pp);
		}
		for(long ss=0;ss<=nn-1;ss++)
		{ 
			if(ss>=kk)
				j1=kk;    //保证yb的下标不越界
			else
				j1=ss;
		  
			v2=0;
			for(long jj=0;jj<=j1;jj++)
				 v2=v2+mod(yq[jj]*ya[ss-jj],pp);
			yr[ss]=mod(yb[ss]-v2,pp);
	   }
	}
	return 1;
}

template <class Type>
Type Mod(Type *x,long n,const Type &p)	//矩阵求模,矩阵以行优先存在向量x中,n为矩阵的阶
{
	for(long i=0;i<n;i++)
		for(long j=0;j<n;j++) 
			x[i*n+j]=mod(x[i*n+j],p);
	return 1;
}

template <class Type>
void Multi(Type *x,Type *y,long n,const Type &p,Type *t)	//矩阵的乘法,结果保存在t中
{
	Type *g=new Type[n*n];
	Mod(x,n,p);
	Mod(y,n,p);
	long i,j,k;
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			g[i*n+j]=0;
			for(k=0;k<n;k++)
				g[i*n+j]+=mod(x[i*n+k]*y[k*n+j],p);
		}
	}
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			t[i*n+j]=mod(g[i*n+j],p);

	delete []g;
	return;
}


//矩阵的模幕,n为阶,fw以行优先保存结果
template <class Type>
void matrixpower(Type *x,long n,const Type &yy,const Type &pp,Type *fw)
{
	if(pp<=0 || yy<0)	//参数错误
		return;

	long i,j;
	long length=sizeof(Type)*8;
	Type *tw=new Type[n*n];	//tw 为单位矩阵
	for (i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			if(i==j)
				tw[i*n+j]=1;
			else
				tw[i*n+j]=0;

			x[i*n+j]=mod(x[i*n+j],pp);
		}
	}

	if(yy!=0)
	{
		for(i=length-1;i>=0;i--)
		{
			Multi(tw,tw,n,pp,tw);
			if(yy & (1<<i))				//依次取yy的每一位
				Multi(tw,xx,n,pp,tw);	
		}
	}

	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			fw[i][j]=tw[i][j];
	}
	return;
}
#endif