p is the field characteristic.


#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>
#include <conio.h>
#include <io.h>
#include <tchar.h> 

#include <math.h>
#include <process.h>
/*Now the great library by Victor Shoup www.shoup.net*/
#include <NTL/ZZ.h>
#include <NTL/vec_ZZ.h>
#include <NTL/mat_ZZ.h>
#include <NTL/GF2EX.h>
#include <NTL/lzz_pEX.h>
#include <NTL/GF2XFactoring.h>
#include <NTL/GF2EXFactoring.h>
#include <NTL/lzz_pXFactoring.h>
#include <NTL/lzz_pEXFactoring.h>
#include <NTL/mat_GF2.h>
#include <NTL/mat_GF2E.h>
#include <NTL/mat_lzz_p.h>
#include <NTL/mat_lzz_pE.h>


#if p>2
	#define GF zz_p
	#define GFE zz_pE
	#define GFX zz_pX
	#define GFEX zz_pEX
	#define GFBak zz_pBak
	#define GFEBak zz_pEBak
	#define GFEContext zz_pEContext
	#define vec_GF vec_zz_p
	#define vec_GFX vec_zz_pEX
	#define vec_GFE vec_zz_pE
	#define vec_GFEX vec_zz_pEX
	#define mat_GF mat_zz_p
	#define mat_GFX mat_zz_pX
	#define mat_GFE mat_zz_pE
	#define mat_GFEX mat_zz_pEX
	#define random_GFE (random_zz_pE())
	#define GFEXModulus zz_pEXModulus
	#define vec_pair_GFEX_long vec_pair_zz_pEX_long
#else
	#define GF GF2
	#define GFE GF2E
	#define GFX GF2X
	#define GFEX GF2EX
	#define GFBak GF2Bak
	#define GFEBak GF2EBak
	#define GFEContext GF2EContext
	#define vec_GF vec_GF2
	#define vec_GFE vec_GF2E
	#define vec_GFEX vec_GF2EX
	#define mat_GF mat_GF2
	#define mat_GFX mat_GF2X
	#define mat_GFE mat_GF2E
	#define mat_GFEX mat_GF2EX
	#define random_GFE (random_GF2E())
	#define GFEXModulus GF2EXModulus
	#define vec_pair_GFEX_long vec_pair_GF2EX_long
#endif







/*returns the number of roots, the correct n of GFE is determined, */
/*the fastest method I found for Hfe*/
s32 OneRoot(GFE &Root, GFEX monHfe)
{
			if(deg(monHfe)<2) return SolveLinearOnly(Root,monHfe);

			s32 i=0,ourn;
			ourn=GFE::degree();
			GFEX Big,Pgcd;
			vec_pair_GFEX_long factors;

			/*compute the X^(p^ourn)-X mod MonHfe*/
			GFEXModulus OptimizedHfe(monHfe);/*precomputation*/
			#if (p>2)/*normal version optimized with GF2EXModulus*/
				Big=GFEX(1,1);
				Big=Big%monHfe;/*needed only if deg(monHfe)=1	<=	deg('X')*/
				for(i=1;i<=ourn;i++)/*computing Big=Big^p mod MonHfe (ourn times) */
					Big=PowerMod(Big,p,OptimizedHfe);/**/
			#else/*faster than repeated squaring when d>n [Kaltofen & Shoup]*/
				Big=FrobeniusMap(OptimizedHfe);
			#endif
			Big-=GFEX(1,1);
			/*Show(Big));/*Show(monHfe));*/

			/*gcd*/
			/*Q: Shoup's GCD hang was probably due before ???*/
			Pgcd=GCD(Big,monHfe);
			
			/*Show(Pgcd));*/
			if (IsOne(Pgcd)) 
				return 0;/*No roots*/
			else
			{
				factors=SquareFreeDecomp(Pgcd);
				Root=FindRoot(factors[0].a);/*you could extract other roots too with factors[i]...*/
				return deg(Pgcd);
			};
}

/*the same BUT different, this one BLOCKS if it has no roots, WHY ?*/
s32 OneRoot2(GFE &Root, GFEX monHfe)
{
			if(deg(monHfe)<2) return SolveLinearOnly(Root,monHfe);

			s32 i=0,ourn;
			ourn=GFE::degree();
			vec_pair_GFEX_long factors;

				factors=SquareFreeDecomp(monHfe);
				if (factors.length()>0)
				{
					Root=FindRoot(factors[0].a);
					return 1987;
				}
				else
					return 0;
}

/*5 times slower than OneRoot*/
s32 OneRoot3(GFE &Root, GFEX monHfe)
{
			if(deg(monHfe)<2) return SolveLinearOnly(Root,monHfe);

			s32 i=0,ourn;
			ourn=GFE::degree();
			GFEX Big,Pgcd;
			vec_pair_GFEX_long factors;

			factors=SquareFreeDecomp(monHfe);
				if (factors.length()>0)
				{
					factors=NewDDF(factors[0].a,Big);
					if (factors[0].b==1)
					{
						Root=FindRoot(factors[0].a);
						return 1987;
					}
					else
						return 0;
				}
				else
					return 0;
}


/*fully independent of external p and n, provided (p?=2 is conserved) and sensible usage */
void SolvingHfeTiming(s32 Hfen, s32 Hfed)
{
	GFEBak bak;s32 PrevGFEextdeg=CurrGFEextdeg;	bak.save();
	BuildGFE(Hfen);
	
	/*random Hfe*/
			GFEX monHfe;
			monHfe=RandomSecretHfe(Hfen,Hfed,1);
			/*Show(monHfe));*/

	/*find one root*/
			GFE Root,Eval;
			s32 nbofroots;
			/**********************************/
			/**********************************/
			struct _timeb tstruct1,tstruct2;
			_ftime( &tstruct1 );
			/**********************************/
			/**********************************/
			nbofroots=OneRoot(Root,monHfe);
			/**********************************/
			/**********************************/
			_ftime( &tstruct2 );
			s32 Inms=1000*(tstruct2.time-tstruct1.time)+tstruct2.millitm-tstruct1.millitm;
			/**********************************/
			/**********************************/
			char buffer[256];
			FILE *speedfl=NULL;
			speedfl=Fopen("speed.log","a+");
			s32 pos=sprintf(buffer,"Hfe GF(%d) n=%d, d=%d, ",p,Hfen,Hfed);
			pos+=sprintf(buffer+pos,"%4.3f s, ",(double)Inms/1000);
			pos+=sprintf(buffer+pos,"%d roots",nbofroots);
			printf("%s\n",buffer);
			if (speedfl)
			{
				fprintf(speedfl,"%s\n",buffer);
				fclose(speedfl);
			};
	/*test*/
			/*Show(Root));
			Eval=eval(monHfe,Root);
			Show(Eval));*/
			if (nbofroots>0) 
				if( !IsZero(eval(monHfe,Root)))
					SayExit(6446,"False Solution to Hfe");



	BuildGFE(PrevGFEextdeg);bak.restore();/*GFE back again*/
}


void RepeatedSolvingHfeTiming(s32 Hfen, s32 Hfed)
{	
	for(s32 i=1;i<=5;i++)
		SolvingHfeTiming(Hfen,Hfed);
}