ercitongyu.hpp

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

//二次同余方程的求解

#include "basic.hpp"

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

template <class Type> 									// 求pp的一个二次非剩余
Type nqr(const Type &pp);
template <class Type> 
long sol2pp0(const Type &aa,const Type &pp,Type *eq);	// 求解方程 x^2 =a(mod p)
template <class Type>                                   // 求解方程 h2*x^2+h1*x+h0=0(mod p)
long sol2pp(const Type &h2,const Type &h1,const Type &h0,const Type &pp,Type *eq);

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

#ifndef INF 
#define INF = -1
#endif

template <class Type> 
Type nqr(const Type &pp)	// 求pp的一个二次非剩余
{
	Type s20,s24,s40,s44;
	Type fw;
	Type bb,hb,p0;

	//下面代码需要进一步了解
	s20=mod(pp,20);
	s24=mod(pp,24);
	s40=mod(pp,40);
	s44=mod(pp,44);
	if(mod(pp,5)==2 || mod(pp,5)==3)
		fw=5;
	else if(mod(pp,6)==3 || mod(pp,6)==5)
	    fw=pp-3;
	else if(mod(pp,12)==5 || mod(pp,12)==7)
	    fw=3;
	else if(s20==5 || s20==11 || s20==13 || s20==15 || s20==17 || s20==19)
		fw=pp-5;
	else if(s24==7 || s24==11 || s24==13 || s24==17)
		fw=6;
	else if(s40==5 || s40==7 || s40==11 || s40==15 || s40==17 || s40==19
			|| s40==21 || s40==23 || s40==25 || s40==29 || s40==33 || s40==35) 
	    fw=10;
	else if(s44==3 || s44==13 || s44==15 || s44==17 || s44==21 || s44==23
			|| s44==27  || s44==29 || s44==31 || s44==41)   
	    fw=11;
	else    
	{
		bb=1;
		hb=1;
		p0=(pp-1)/2;
	    while(hb!=pp-1 && bb<pp)
		{
			bb=bb+1;
			hb=power(bb,p0,pp);     
		}   
	   if(hb==pp-1)
		   fw=bb;      //满足 pp-1==power(fw,(pp-1)/2,pp),欧拉准则：pp 为奇素数 
	   else            //此时二次同余方程 x^2 = fw (mod pp)  无解
		   fw=0;
	}
	return fw;
}


//用构造法求解二次同余方程x^2=a(mod p),出错返回INF
template <class Type> 
long sol2pp0(const Type &a,const Type &pp,Type *eq)
{
	Type aa=mod(a,pp);
	if((pp<=0) || (mod(pp,2)==0))
	{
		eq[0]=-1;
		return INF;
	}
	else if(mod(aa,pp)==0) // x^2 = 0 (mod p)
	{
		eq[0]=2;
		eq[1]=0;
		eq[2]=0;
	}
	else  //mod(aa,pp)!=0
	{
	   if(power(aa,(pp-1)/2,pp)==pp-1)   //由欧拉准则 (a/p)=a^((p-1)/2)，无解
	   {
		   eq[0]=0;
		   return 1;
	   }
	   else   
	   {
		   Type kk=1;
		   Type ll=(pp-1)/2;
		   while(mod(ll,2)==0)   //参看<离散数学> 李盘林 李丽双   高等教育出版社 p173
		   {
			   kk=kk+1;
			   ll=ll/2;
		   }                         //2^kk*ll==pp-1
  
		   if(kk==1)				// p mod 4 = 3, x = a^((p+1)/4)
		   {
			   Type hh=(pp+1)/4;        
			   Type fw=power(aa,hh,pp);
			   eq[0]=2;
			   eq[1]=fw;
			   eq[2]=pp-fw;
		   }
		   else if(kk==2)			// p mod 8 = 5
		   {
			   Type hh=(pp+3)/8;  //用作指数  
			   Type fw;			  //保存一个解
			   Type fw0=power(aa,hh,pp);
			   Type hh0=(pp-1)/4;
			   Type gg=power(aa,hh0,pp); //用作判断
			   if(gg==1)			// p mod 8 = 5 && a^((p-1)/4)=1,x=a^((p+3)/8)
				   fw=fw0;
			   else   //gg==-1		// p mod 8 = 5 && a^((p-1)/4)=-1,x=2^((p-1)/4)*a^((p+3)/8)
			   {
				   Type gg1=power(2,hh0,pp);   // p = (+-)3 mod 8时,2是p的二次非剩余
				   fw=mod(fw0*gg1,pp);    // (2^((p-1)/4))^2=-1
			   }                         
			   eq[0]=2;
			   eq[1]=fw;
			   eq[2]=pp-fw;
		   }
		   else if(kk==3 && power(aa,(pp-1)/4,pp)==1)
		   {
			   Type hh=(pp+7)/16;               //pp=16*k+7,k为正整数
			   Type fw;			 //保存一个解
			   Type fw0=power(aa,hh,pp);  
			   Type hh0=(pp-1)/8;
			   Type gg=power(aa,hh0,pp); 
			   if(gg==1)       //x=a^((p+7)/16)为其一个解
				   fw=fw0;
			   else
			   {
				   Type bb=nqr(pp);  //b为p的一个二次非剩余
				   if(bb!=0)    //(b^((p-1)/4))^2=-1,(a^((p+7)/16))^2=-a
				   {
					   Type gg1=power(bb,hh0*2,pp); 
				   	   fw=mod(fw0*gg1,pp);  //x=(a^((p+7)/16))*(b^((p-1)/4))
				   }
				   else
				   {
					   eq[0]=-1;   //错误  
					   return INF;
				   }
			   }
			   eq[0]=2;
		 	   eq[1]=fw;
			   eq[2]=pp-fw;
		   }
		   else //kk>=4 || kk==3 && power(aa,(pp-1)/4,pp)!=1
		   {
			   Type bb=nqr(pp);   //b为p的一个二次非剩余
			   if(bb==0)
			   {
				   eq[0]=-1;		//错误
				   return INF;
			   }
			   else  
			   {
				   Type cc=power(aa,ll,pp);
				   Type gg=power(bb,ll,pp);
				   Type ff=mod(gg*gg,pp); 
				   Type ff0=1;
				   Type tt=0;
				   Type temp=power(2,kk-1);
				   while(mod(ff0-cc,pp)!=0 && tt<=temp)
				   {
					   ff0=mod(ff*ff0,pp);
					   tt=tt+1;
				   }
				   //存在t<=2^k, b^(2lt)=a^l,其中l满足 p = 2^k*l+1
				   if(mod(ff0-cc,pp)!=0)
				   {
					   eq[0]=-1;		//不存在上述t,错误!
					   return INF;
				   }
				   else  //离散数学 p174
				   {
					   Type gg1=ninv(gg,pp);			//g1=(b^l)^(-1)
					   Type cc1=power(aa,(ll+1)/2,pp);  //c1=a^((l+1)/2)
					   Type fw0=power(gg1,tt,pp);		//fw0=(b^l)^(-t)
					   Type fw=mod(fw0*cc1,pp);			//x=(b^(-lt))*a^((l+1)/2),注意b^(2lt)=a^l
					   eq[0]=2;
					   eq[1]=fw;
					   eq[2]=pp-fw;
				   } 
			  }
		  }
	  }
	}
	return 1;
}

/*
% 求解同余方程 f(x)=h2*x^2+h1*x+h0==0(mod pp) 或者等价方式 x^2+va*x+vb==0(mod pp)
% 结果用数组 eq 表示，eq[0]存放解的个数，eq[0]=0即无解，eq[0]=-1 非正常退出
% eq[0]=1 时 eq[1]存放唯一解,eq[0]=2 时 eq[1],eq[2]存放两个解。
% 要调用子函数 ninv : 求整数模逆 ;power 求整数模幂 ;nqr：二次非剩余候选数
*/
template <class Type>
long sol2pp(const Type &hh2,const Type &hh1,const Type &hh0,const Type &pp,Type *eq)
{
	Type h2=mod(hh2,pp),h1=mod(hh1,pp),h0=mod(hh0,pp);
	if(pp<=0 || mod(pp,2)==0)
	{
		eq[0]=-1;
		return INF;
	}
	else if(mod(h2,pp)==0)
	{
		if(mod(h1,pp)==0)
		{
			eq[0]=-1;
			return INF;
		}
		else
		{
			eq[0]=1;eq[1]=mod(-ninv(h1,pp)*h0,pp);  // 化为h1*x+h0=0(mod p)
		}
	}
	else   //mod(h2,pp)!=0
	{
		Type v0=ninv(h2,pp);
		Type va=mod(v0*h1,pp);
		Type vb=mod(v0*h0,pp);						//化为 x^2+va*x+vb=0(mod p)
		Type dt=mod(va*va-4*vb,pp);					//判别式=va^2-4*1*vb
		Type iv2=ninv(Type(2),pp);   
/*		if(mod(dt,pp)==0)
		{
			eq[0]=2;eq[1]=mod(iv2*(-va),pp);eq[2]=eq[1];   //x=-va/2
		}
		else if(power(dt,(pp-1)/2,pp)==pp-1)	//令y=2*x+va,化为y^2=dt (mod p)
		{
			eq[0]=0;
		}*/
//		else 
		{
			sol2pp0(dt,pp,eq);					//求解 y^2=dt(mod p)
			if(eq[0]>0)
			for(long kk=1;kk<=eq[0];kk++)
				eq[kk]=mod(iv2*(eq[kk]-va),pp);	// x=(y-va)/2
		}
	} 
	return 1;
}