/* $OpenXM: OpenXM_contrib2/asir2018/lib/primdec_mod,v 1.2 2021/02/18 05:35:01 noro Exp $ */

extern First_Component,Minipoly_SBA,NDGR$
extern Hom,GBTime$
extern DIVLIST,INTIDEAL,ORIGINAL,ORIGINALDIMENSION,STOP,Trials,REM$
extern T_GRF,T_INT,T_PD,T_MP$
extern BuchbergerMinipoly,PartialDecompByLex,ParallelMinipoly$
extern B_Win,D_Win$
extern COMMONCHECK_SF,CID_SF$

if (!module_definedp("fff")) load("fff"); else $
if (!module_definedp("gr")) load("gr"); else $
module primdec_mod $
  /* Empty for now. It will be used in a future. */
endmodule $

/*==============================================*/
/*  prime decomposition of ideals over 		*/
/*  finite fields				*/
/*        part 1 				*/
/*	radical computation 			*/
/*  28/12/2002	for Basic Small Finite Field	*/
/*   4/1/2003	for Arbitary Small Finite Field	*/
/*==============================================*/

/*==============================================*/
/*  radical computation of ideals over 		*/
/*  finite fields by Matsumoto's algorithm	*/
/*						*/
/* The radical of I is computed as the 		*/
/* kernel of a power of Frobenius map.		*/
/* 						*/
/*  radical 					*/
/*	Input: a list P of polynomials		*/
/*	Output: a list P1 of polynomials 	*/
/*	     such that <P1>=rad(<P>)		*/
/*						*/
/*  frobeniuskernel_main{2}			*/
/*	Input: a list of polynomials P,		*/
/* 	       a list of variables VSet,	*/
/*	   and a list of other variables WSet	*/
/*	Output: a list of polynomials Q		*/
/*	   such that Q is the kernel of 	*/
/*         single Frobenius map: x--> x^q	*/
/*	   where q is the characteristic.	*/
/*	version 2 uses successive kernel	*/
/*	computation suggested by Matsumoto	*/
/*						*/
/*  coefficientfrobeniuskernel			*/
/*	Input: a list of polynomials P,		*/
/*	Output: a list of polynomials Q		*/
/*	   such that Q is the kernel of 	*/
/*         single coefficient Frobenius 	*/
/*	   map: a in GF(q) --> a^q		*/
/*	   where q is the characteristic.	*/
/*						*/
/*  elimination					*/
/*	Input: a list P of polynomials		*/
/*	       a list V of variables		*/
/*	Output: a list P1 of polynomials 	*/
/*  		such that P1={f\in P | 		*/
/*		  vars(f)\cap V =\emptyset}	*/
/*						*/
/*  checkequality				*/
/*	Input: a list P of polynomials		*/
/*	       a list Q of polynomials		*/
/*	Output: 1 if <P>=<Q>			*/
/*		0 otherwise			*/
/*						*/
/*						*/
/*  						*/
/*==============================================*/

def radicalideal(P,Choice,VSet)
	{

	GBTime=0;
	Ord=0;

	/* VSet=vars(P);*/
	WSet=makecounterpart(VSet);

T000=time()[0];
	P1=dp_gr_f_main(P,VSet,Hom,Ord);
T001=time()[0];
GBTime +=T001-T000;

	if ( Choice == 3 )
		{
		P2=frobeniuskernel_main3(P1,VSet,WSet);
		}
	else if ( Choice == 2 )
		{
		P2=frobeniuskernel_main2(P1,VSet,WSet);	
		}
	else if ( Choice == 4 )
		{
		P2=frobeniuskernel_main4(P1,VSet,WSet);	
		}
	else  
		{
		P2=frobeniuskernel_main(P1,VSet,WSet);	
		}

	M=checkequality(P1,P2,VSet);

	while ( M !=1 )
		{
		P1=P2;
		P2=frobeniuskernel_main2(P2,VSet,WSet);
		M=checkequality(P1,P2,VSet);
		}

	return P1;
	}


def frobeniuskernel_main(P,VSet,WSet)
	{
	NV=length(VSet);

	NewP=coefficientfrobeniuskernel(P);

	NW=length(NewP);

	TempP=NewP;

	for (K=NW-1;K>=0;K--)
		{
		Poly=TempP[K];
		for (J=0;J<NV;J++)
			{
			Poly=subst(Poly,VSet[J],WSet[J]);
			}
		NewP=cons(Poly,NewP);
		}

	XSet=append(VSet,WSet);
	NewOrder=[[0,length(VSet)],[0,length(WSet)]];

	Char=characteristic_ff();

	for (I=0;I<NV;I++)
		{
		AddPoly=VSet[I]^Char-WSet[I];
		NewP=cons(AddPoly,NewP);
		}

	G=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
	PW=elimination(G,WSet);

	NW=length(PW);
	ANS=[];
	for (I=NW-1;I>=0;I--)
		{
		Poly=PW[I];
		for (J=0;J<NV;J++)
			{
			Poly=subst(Poly,WSet[J],VSet[J]);
			}
		ANS=cons(Poly,ANS);
		}
	return ANS;
	}

def frobeniuskernel_main2(P,VSet,WSet)
	{
	NV=length(VSet);

	NewP=coefficientfrobeniuskernel(P);

	XSet=append(VSet,WSet);
	NewOrder=[[0,NV],[0,NV]];

	Char=characteristic_ff();

	for (I=0;I<NV;I++)
		{
		for (J=0;J<NV;J++)
			{
			if ( I == J )
				{
				AddPoly=VSet[J]^Char-WSet[J];
				}
			else
				{
				AddPoly=VSet[J]-WSet[J];
				}			
			NewP=cons(AddPoly,NewP);
			}

T000=time()[0];
		GP=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
T001=time()[0];
GBTime += T001-T000;

		PW=elimination(GP,WSet);
		NW=length(PW);
		NewP=[];
		for (K=NW-1;K>=0;K--)
			{
			Poly=PW[K];
			for (J=0;J<NV;J++)
				{
				Poly=subst(Poly,WSet[J],VSet[J]);
				}
			NewP=cons(Poly,NewP);
			}
		}
	return NewP;
	}


def frobeniuskernel_main4(P,VSet,WSet)
	{
	NV=length(VSet);

	NewP=coefficientfrobeniuskernel(P);

	XSet=append(VSet,WSet);
	NewOrder=[[0,NV],[0,NV]];

	Char=characteristic_ff();

	for (I=0;I<NV;I++)
		{
		NW=length(NewP);
		TempP=NewP;

		for (J=0;J<NV;J++)
			{
			if ( I == J )
				{
				AddPoly=VSet[J]^Char-WSet[J];
				}
			else
				{
				AddPoly=VSet[J]-WSet[J];
				}			
			NewP=cons(AddPoly,NewP);
			}

		/* for (K=NW-1;K>=0;K--)
			{
			Poly=TempP[K];
			for (J=0;J<NV;J++)
				{
				if ( J == I )
					{
					Poly=subst(Poly,VSet[J],WSet[J]);
					}
				else
					{
					Poly=subst(Poly,VSet[J],WSet[J]^Char);
					}
				}
			NewP=cons(Poly,NewP);
			}*/
T000=time()[0];
		GP=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
T001=time()[0];
GBTime += T001-T000;

		PW=elimination(GP,WSet);
		NW=length(PW);
		NewP=[];
		for (K=NW-1;K>=0;K--)
			{
			Poly=PW[K];
			for (J=0;J<NV;J++)
				{
				Poly=subst(Poly,WSet[J],VSet[J]);
				}
			NewP=cons(Poly,NewP);
			}
		}
	return NewP;
	}

def frobeniuskernel_main3(P,VSet,WSet)
	{
	NV=length(VSet);

	NewP=coefficientfrobeniuskernel(P);

	Char=characteristic_ff();

	for (I=0;I<NV;I++)
		{
		XWSet=[];
		for (J=0;J<NV;J++)
			{
			if ( I == J )
				{
				AddPoly=VSet[I]^Char-WSet[I];
				NewP=cons(AddPoly,NewP);
				XWSet=append(XWSet,[WSet[I]]);
				}
			else
				{
				XWSet =append(XWSet,[VSet[J]]);
				}			
			}

		XSet=cons(VSet[I],XWSet);

		NewOrder=[[0,1],[0,NV]];
T000=time()[0];
		GP=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
T001=time()[0];
GBTime += T001-T000;

		PW=elimination(GP,XWSet);
		NW=length(PW);
		NewP=[];
		for (K=NW-1;K>=0;K--)
			{
			Poly=PW[K];
			Poly=subst(Poly,WSet[I],VSet[I]);
			NewP=cons(Poly,NewP);
			}
		}
	return NewP;
	}

def coefficientfrobeniuskernel(P)
	{
	if ( setmod_ff()[1] == 0 )
		{
		return P;
		}

	NP=length(P);
	ANS=[];
	for (I=0;I<NP;I++)
		{
		TEST=P[I];
		Q=coefficientfrobeniuskernel_main(TEST);
		ANS=append(ANS,[Q]);
		}
	return ANS;
	}

def coefficientfrobeniuskernel_main(Poly)
	{
	Vars=vars(Poly);
	QP=dp_ptod(Poly,Vars);
	ANS=0;
	FOrd=extdeg_ff();
	Char=characteristic_ff();
	Pow=Char^(FOrd-1);

	while(QP !=0 )
		{
		Coef=dp_hc(QP);
		HT=dp_ht(QP);

		ANS=ANS+Coef^Pow*dp_dtop(HT,Vars);

		QP=dp_rest(QP);
		}
	return ANS;
	}


def elimination(G,V)
	{
	ANS=[];
	NG=length(G);

	for (I=NG-1;I>=0;I--)
		{
		VSet=vars(G[I]);
		DIFF=setminus(VSet,V);

		if ( DIFF ==[] ) 
			{
			ANS=cons(G[I],ANS);
			}
		}
	return ANS;
	}

def setminus(A,B)
	{
	NA=length(A);
	NB=length(B);

	ANS=[];

	for (I=0;I<NA;I++)
		{
		Flag =0;
		for (J=0;J<NB;J++)
			{
			if (A[I] == B[J])
				{
				Flag=1;
				break;
				}
			}

		if ( Flag == 0 )
			{
			ANS=append(ANS,[A[I]]);
			}
		}
	return ANS;
	}

def makecounterpart(V)
	{
	NV=length(V);

	A="tmp";

	ANS=[];
	for (I=NV-1;I>=0;I--)
		{
		T=strtov(A+rtostr(I));
		ANS=cons(T,ANS);
		}
	return ANS;
	}

def checkequality(P,Q,VSet)
	{
	QQ=dp_gr_f_main(Q,VSet,Hom,Ord);
	PP=dp_gr_f_main(P,VSet,Hom,Ord);

	NP=length(PP);
	NQ=length(QQ);

	VarsP=vars(P);
	VarsQ=vars(Q);
	if ( NP != NQ || VarsP !=VarsQ )
		{
		return 0;
		}

	for (I=0;I<NP;I++)
		{
		HCP=dp_hc(dp_ptod(PP[I],VSet));
		HCQ=dp_hc(dp_ptod(QQ[I],VSet));

		if ( HCP*QQ[I]-HCQ*PP[I] != 0 )
			{
			return 0;
			}
		}

	return 1;
	}

/*==============================================*/
/*  prime decomposition of ideals over          */
/*  finite fields                               */
/*        part 2                                */
/*      ideal dimension		                */
/*        3/1/2003                    		*/
/*==============================================*/

/*==============================================*/
/*  ideal dimension computation 	        */
/*                                              */
/* The dimension of I is computed as the       	*/
/* maximal size of MSI sets.          		*/
/*                                              */
/*  idealdimension                              */
/*      Input: a list P of polynomials          */
/*	       a list V of variables		*/
/*      Output: the dimension of I 		*/
/*           and the MSI with max size		*/
/*                                              */
/*  findmsi		                        */
/*      Input: a list of head monomials         */
/*         and a list of variables 		*/
/*      Output: a list of all MSI sets		*/
/*             and SI sets			*/
/*                                              */
/*  findmsi_main				*/
/*      Input: a list P of polynomials          */
/*             a list V of variables            */
/*      Output: a list P1 of polynomials        */
/*              such that P1={f\in P |          */
/*                vars(f)\cap V =\emptyset}     */
/*                                              */
/*  headtermset					*/
/*      Input: a list P of polynomials          */
/*             a list V of variables		*/
/*      Output: a list of head monomials 	*/
/*                                              */
/*                                              */
/*                                              */
/*==============================================*/

def idealdimension(P,V,Ord)
	{
	GP=dp_gr_f_main(P,V,Hom,Ord);
	HeadTermset=headtermset(GP,V,Ord);
	MSIset=findmsi(HeadTermset,V);
	if ( MSIset == [] )
		{
		return [0,[]];
		}
	Index=findmaxindex(MSIset);
	MSI=MSIset[Index];
	Dimension=length(MSI);
	return [Dimension, MSI];
	}

def findmaxindex(S)
	{
	NS=length(S);
	Index=0;
	Defaultsize=length(S[0]);
	for (I=1;I<NS;I++)
		{
		Targetsize=length(S[I]);

		if ( Targetsize > Defaultsize )
			{
			Index=I;
			}
		Defaultsize= Targetsize;
		}
	return Index;
	}

def headtermset(P,V,Ord)
	{
	ANS=[];
	NP=length(P);
	for ( I=NP-1;I>=0;I--)
		{
		Head=headterm(P[I],V,Ord);
		ANS=cons(Head,ANS);
		}
	return ANS;
	}

def headterm(P,V,Ord)
	{
	dp_ord(Ord);
	Q=dp_ptod(P,V);
	Headdp=dp_hm(Q);
	Head=dp_dtop(Headdp,V);
	return Head;
	}

def findmsi_main(TermsetIndex,MSIsetIndex,Int,N)
	{
	ANS=[];
	BNS=[];
	TempMSIsetIndex=MSIsetIndex;

	if (Int == 0)
		{
		for ( I=0;I<N;I++)
			{
			TEST1=[I];
			Check=intersection(TermsetIndex,TEST1);
			if ( Check == 0 )
				{
				ANS=cons(TEST1,ANS);
				}
			}
		}

	while( Int !=0 && TempMSIsetIndex != [] )
		{
		TEST=TempMSIsetIndex[0];
		Flag=0;
		for ( I=TEST[0]+1;I<N;I++)
			{
			TEST1=cons(I,TEST);
			Check=intersection(TermsetIndex,TEST1);
			if ( Check == 0 )
				{
				Flag=1;
				ANS=cons(TEST1,ANS);
				}
			}
		if ( Flag == 0 )
			{
			BNS=cons(TEST,BNS);
			}
		TempMSIsetIndex=cdr(TempMSIsetIndex);
		}
	return [ANS,BNS];
	}
		
def findmsi(Termset,Vset)
	{
	ANS=[];
	MSIsetIndex=[];
	N=length(Vset);
	TermsetIndex=makeindex(Termset,Vset);

	for (I=0;I<N;I++)
		{
		Temp=findmsi_main(TermsetIndex,ANS,I,N);
		ANS=Temp[0];
		BNS=Temp[1];
		MSIsetIndex=append(BNS,MSIsetIndex);
		}

	MSIsetIndex=append(ANS,MSIsetIndex);
	MSIset=makevarset(MSIsetIndex,Vset);
	return MSIset;
	}

def makeindex(Termset,Vset)
	{
	ANS=[];

	N=length(Vset);

	TempTermset=Termset;
	while( TempTermset !=[] )
		{
		TEST=TempTermset[0];
		Index=[];
		for (I=0;I<N;I++)
			{
			Q=srem(TEST,Vset[I]);
			if ( Q == 0 )
				{
				Index=cons(I,Index);
				}
			}
		ANS=cons(Index,ANS);
		TempTermset=cdr(TempTermset);
		}
	return ANS;
	}

def makevarset(IndexSet,Vset)
	{
	ANS=[];
	NI=length(IndexSet);
		
	for (I=0;I<NI;I++)
		{
		TEST=IndexSet[I];

		TestV=makevar(TEST,Vset);
		ANS=append(ANS,[TestV]);
		}
	return ANS;
	}

def makevar(Index,Vset)
	{
	ANS=[];
	TempIndex=Index;

	while( TempIndex !=[] )
		{
		ANS=cons(Vset[TempIndex[0]],ANS);
		TempIndex=cdr(TempIndex);
		}

	return ANS;
	}

def intersection(IndexSet,TestIndex)
	{
	TempIndexSet=IndexSet;

	while(TempIndexSet !=[] )
		{
		Sample=TempIndexSet[0];

		Inclusion=testinclusion(Sample,TestIndex);

		if ( Inclusion == 1 )
			{
			return 1;
			}

		TempIndexSet=cdr(TempIndexSet);
		}
	return 0;
	}

def testinclusion(Sample,Test)
	{
	NS=length(Sample);
	NT=length(Test);
	
	for (I=0;I<NS;I++)
		{
		Flag=0;
		for (J=0;J<NT;J++)
			{
			if (Sample[I]==Test[J])
				{
				Flag=1;
				break;
				}
			}
		if ( Flag == 0 )
			{
			return 0;
			}
		}
	return 1;	
	}


/*======================================================*/
/*  prime decomposition of ideals over          	*/
/*  small finite fields                         	*/
/*        part 3                                	*/
/*      prime decomposition                     	*/
/*	  9/1/2003  1st Implementation			*/
/*	  11/1/2003  Divid 3 cases in zeroprime		*/
/*		decomposition		  (3a)		*/
/*	  12/1/2003 Redundant elimination (3e)		*/
/*	  13/1/2003 gr_fctr_sf version    (3f)		*/
/*	  14/1/2003 Flag on Early Termination		*/
/*	  15/1/2002 drl version				*/
/*	  17/1/2002 depth oriented search 		*/
/*	  17/1/2002 dimension oriented search 		*/
/*							*/
/*  global variables					*/
/*  Hom: used in dp_gr_f_main 				*/
/*  DIVLIST: the set of prime ideals			*/
/*  ORIGINAL: the GB of the original input ideal	*/
/*  ORIGINALDIMENSION: the dimension 			*/
/*  STOP: the flag on terminatation			*/
/*  Trials: a bound on the number of trials             */
/*  REM: the vector list of remaining ideals 		*/
/*							*/
/*======================================================*/

/*======================================================*/
/*  	prime decomposition				*/
/*                                              	*/
/*                                              	*/
/*  primedec(P,V,Ord,F) 		                */
/*      Input: a list P of polynomials          	*/
/*             a list V of variables            	*/
/*	       a term order Ord				*/
/*	       a flag F on Early Termination		*/
/*      Output: 0					*/
/*     DVILIST: the set of prime divisors of <P>	*/
/*  [subprocedure]	                        	*/
/*	 o --gr_fctr_sf    (mplex2)			*/
/*	 o --primedecomposition 			*/
/*  Remark: if F = 1, we empoly Early Termination	*/
/*							*/
/*							*/
/*  primedecomposition(P,V,Ord,C,F)                 	*/
/*      Input: a list P of polynomials          	*/
/*             a list V of variables            	*/
/*	       a term order Ord				*/
/*	       a level C of recursive call		*/
/*	       a flag F on Early Termination		*/
/*      Output: 0					*/
/*     DVILIST: the set of prime divisors of <P>^ec	*/
/*  [subprocedure]	                        	*/
/*	 o --idealdimension (part2)			*/
/*	 o --setminus	   (part1)			*/
/*	 o --zeroprimedecomposition 			*/
/*	 o --extcont					*/
/*	 o --checkadd2					*/
/*	 o --ideal_intersection_sfrat 			*/
/*	 o --radical_equality				*/
/*							*/
/*  zeroprimedecomposition(P,V,W)               	*/
/*      Input: a list P of polynomials          	*/
/*             	lists V,W of variables           	*/
/*	       	such that <P> is 0-dimensional		*/
/*	       	in K[V] and W is the original set	*/
/*	       	of variables				*/
/*      Output: the set of prime/primary 		*/
/*		divisors of <P>				*/
/*		(P is a GB w.r.t. DRL)			*/
/*  [subprocedure]	                        	*/
/*	 o --partical_decomp (mplex2)			*/
/*	 o --checkgeneric				*/
/*	 o --separableclosure				*/
/*	 o --zerogenericprimedecomposition		*/
/*	 o --zeroseparableprimedecomposition		*/
/*	 o --convertdivisor				*/
/*	 o --contraction				*/
/*	 o --radicalideal (part1)			*/
/*                                              	*/
/*  extcont(P,V,Ord)         				*/
/*      Input: a list P of polynomials			*/
/*             	a set V of variables    		*/
/*      Output: a polynomial f				*/
/*             	such that \sqrt<P>^c=\sqrt(P:f)		*/
/*	    Then \sqrt<P>=\sqrt(P:f)\cap \sqrt(P+<f>)	*/
/*                                              	*/
/*  separableclosure([P,MP],V,W)			*/
/*      Input: a pair [P,MP] of a list P of 		*/
/*	       	polynomials and a list MP of 		*/
/*	       	minimal polynomials of variables	*/
/*             	list V,W of variables          		*/
/*	       	such that <P> is 0-dimensional 		*/
/*	       	in K[V] and W is the original set	*/
/*	       	of variables				*/
/*      Output: [Q,E], a list Q of polynomials		*/
/*		 an exponent vector E			*/
/*             	 such that <Q> is the separable 	*/
/*	       	 closure of <P>				*/
/*  [subprocedure]	                        	*/
/*	 o --makecounterpart(part1)			*/
/*	 o --elimination(part1)				*/
/*	 o --checkseparable				*/
/*                                              	*/
/*  zeroseparableprimedecomposition(P,V,W)		*/
/*      Input: a list P of polynomials          	*/
/*             lists V,W of variables           	*/
/*	       	such that <P> is 0-dimensional		*/
/*	       	and all minimal polynomial 		*/
/*	       	are irreducible (i.e. <P> is an		*/
/*	       	intermediate ideal in the paper		*/
/*             	or the output of partial_decomp)	*/
/*	       	in K[V] and W is the original set	*/
/*	       	of variables				*/
/*      Output: the set of prime/primary 		*/
/*		 divisors of <P>			*/
/*  [subprocedure]	                        	*/
/*	 o --findgeneric				*/
/*	 o --elimination(part1)				*/
/*                                              	*/
/*  zerogenericprimedecomposition([P,MP],M,V,W) 	*/
/*      Input: a pair [P,MP] of a list P of 		*/
/*	       	polynomials and a list MP of 		*/
/*	       	minimal polynomials of variables	*/
/*	       	such that <P> is 0-dimensional		*/
/*	       	and M is a minimal poly of a 		*/
/*	       	variable in generic position 		*/
/*	       	in K[V] and W is the original 		*/
/*	       	set of variables			*/
/*      Output: the set of prime/primary 		*/
/*		divisors of <P>				*/
/*  [subprocedure]	                        	*/
/*                                              	*/
/*  convertdivisor(P,V,W,E)                     	*/
/*      Input: a list P of polynomials          	*/
/*             	list V,W of variables          		*/
/*	       	such that <P> is a primary/prime	*/
/*	       	0-dimensional ideal in K[V] and 	*/
/*	       	W is the original set			*/
/*	       	of variables				*/
/*	       	the exponent vector E			*/
/*      Output: the corresponding prime ideal		*/
/*             in K[W]		                	*/
/*  [subprocedure]	                        	*/
/*	 o --elimination(part1)				*/
/*                                              	*/
/*  contraction(P,V,W)	                        	*/
/*      Input: a list P of polynomials          	*/
/*             	list V,W of variables          		*/
/*	       	such that <P> is a primary/prime	*/
/*	       	0-dimensional ideal in K[V] and 	*/
/*	       	W is the original set			*/
/*	       	of variables				*/
/*      Output: the contraction <P>^c in K[W]		*/
/*                                              	*/
/*  findgeneric(P,V)	                        	*/
/*      Input: a list P of polynomials          	*/
/*             	a list V of variables          		*/
/*	       	such that <P> is a separable 		*/
/*	       	0-dimensional ideal in K[V]  		*/
/*      Output: [F,Min,X], a polynomial F in 		*/
/*    	       	generic position, its minimal 		*/ 
/*             	polynomial Min in a new variable 	*/
/*             	X                                	*/
/*  [subprocedure]	                        	*/
/*	 o --lineardimension				*/
/*	 o --makemagic					*/
/*	 o --minipoly_sf (mplex2)			*/
/*							*/
/*======================================================*/

def primedec_mod(P,VSet,Ord,Mod,Strategy)
{
    First_Component = getopt(first);
    if ( type(First_Component) == -1 ) First_Component = 0;
    Minipoly_SBA = getopt(sba);
    if ( type(Minipoly_SBA) == -1 ) Minipoly_SBA = 0;
    NDGR = getopt(ndgr);
    if ( type(NDGR) == -1 ) NDGR = 0;
	for ( Q = Mod, E = 1; Q < 2^14; Q *= Mod, E++ );
	Q /= Mod;
	E--;
	if ( !(E%2) ) {
		Q /= Mod;
		E--;
	}
	setmod_ff(Mod,E);
	Pq = map(simp_ff,P);
	primedec_sf(Pq,VSet,Ord,Strategy);
	R = convsf(DIVLIST,VSet,0,1);
	return R;
}


def primedec_sf(P,VSet,Ord,Strategy)
	{
	/* DIVLIST = the set of all computed candidates for iredundant divisors */
	/* INTIDEAL = the intersection of all computed candidates	 	*/
	/* ORIGINAL = a Groebner basis of the input ideal <P> w.r.t. Ord	*/
	/* ORIGINALDIMENSION = the dimension of the input ideal <P>		*/
	/* STOP = the flag on termination					*/
	/* REM = the list of remaining 	ideals					*/

	T0 = time();
	T_GRF=0;
	T_INT=0;
	T_PD=0;
	T_MP=0;
	DIVLIST=[];
	INTIDEAL=[];
	STOP=0;
	B_Win=0;
	D_Win=0;

	ORIGINAL=dp_gr_f_main(P,VSet,Hom,Ord);

	Dimeset=idealdimension(ORIGINAL,VSet,Ord);
	ORIGINALDIMENSION=Dimeset[0]; 

	/* REM0=newvect(ORIGINALDIMENSION+1);*/
	REM=newvect(ORIGINALDIMENSION+1);

	for (I=0;I<ORIGINALDIMENSION+1;I++)
		{
		/* REM0[I]=[];*/
		REM[I]=[];
		}

	if ( dp_gr_print() ) {
		print("The dimension of the ideal is ",2);print(ORIGINALDIMENSION,2);
		print(". ");
	}

	if ( ORIGINALDIMENSION == 0 )
		{
		Strategy = 0;
		STOP = 0;
		}

 	ANS=gr_fctr_sf([ORIGINAL],VSet,Ord);
	NANS=length(ANS);
	if ( dp_gr_print() ) {
		print("There are ",2);print(NANS,2);print(" partial components. ");
    }
	for (I=0;I<NANS;I++)
		{
		TempI=ANS[I];
		DimI=idealdimension(TempI,VSet,Ord)[0];	

		REM[ORIGINALDIMENSION-DimI]=cons(TempI,REM[ORIGINALDIMENSION-DimI]);
		}

	REM[0]=ANS;

	INDEX=0;

	while  ( INDEX != -1 )
		{
		/* print("We deal with the ",2);print(I,2);print("-th partial component. ");*/

		TESTIDEAL=car(REM[INDEX]);
		REM[INDEX]=cdr(REM[INDEX]);
		primedecomposition(TESTIDEAL,VSet,Ord,INDEX,Strategy);

		if (STOP == 1 )
			{
			DIVLIST = prime_irred_sf_by_first(DIVLIST,VSet,0);
			DIVLIST = monic_sf_first(DIVLIST,VSet);
			if ( dp_gr_print() ) {
				print("We finish the computation. ");
				T_TOTAL = time()[3]-T0[3];
				print(["T_TOTAL",T_TOTAL,"T_GRF",T_GRF,"T_PD",T_PD,"T_MP",T_MP,"T_INT",T_INT,"B_Win",B_Win,"D_Win",D_Win]);
			}
			return 0;
			}

		INDEX=findindex(REM);
		}

	DIVLIST = prime_irred_sf_by_first(DIVLIST,VSet,0);
	DIVLIST = monic_sf_first(DIVLIST,VSet);
	T_TOTAL = time()[3]-T0[3];
	if ( dp_gr_print() ) {
		print(["T_TOTAL",T_TOTAL,"T_GRF",T_GRF,"T_PD",T_PD,"T_MP",T_MP,"T_INT",T_INT,"B_Win",B_Win,"D_Win",D_Win]);
	}
	return 0;
	}





def findindex(VECTORLIST)
	{
	NL=size(VECTORLIST)[0];

	for (I=0;I<NL;I++)
		{
		if ( VECTORLIST[I] == [] )
			{
			continue;
			}
		return I;
		}
	return -1;
	}

#if 0 
def checkadd2(ADD,TESTADDLIST,VSet)
	{
	/* This function will be used for eliminating redundant divisors	*/

	NTESTADDLIST=length(TESTADDLIST);
	for (I=0;I<NTESTADDLIST;I++)
		{
		TESTLIST=TESTADDLIST[I][0];

 		if ( setminus(TESTLIST,ADD) == [] )
			{
			return 1;
			}
		if ( inclusion_test(TESTLIST,ADD,VSet,Ord) == 1 )
			{
			return 1;
			}
		}
	return 0;
	}
#endif 

def checkadd2a(ADD,DIM,TESTADDLIST,VSet)
	{
	/* This function will be used for eliminating redundant divisors	*/

	NTESTADDLIST=length(TESTADDLIST);
	for (I=0;I<NTESTADDLIST;I++)
		{
		TESTLIST=TESTADDLIST[I][0];
		TESTDIM=TESTADDLIST[I][1];

		if (DIM > TESTDIM )
			{
			continue;
			}

 		if ( setminus(TESTLIST,ADD) == [] )
			{
			return 1;
			}
		if ( inclusion_test(TESTLIST,ADD,VSet,Ord) == 1 )
			{
			return 1;
			}
		}
	return 0;
	}

def checkadd3(ADD,TESTADDLIST,VSet)
	{
	/* This function will be used for eliminating redundant divisors	*/

	NTESTADDLIST=length(TESTADDLIST);
	for (I=0;I<NTESTADDLIST;I++)
		{
		TESTLIST=TESTADDLIST[I];

 		if ( setminus(TESTLIST,ADD) == [] ) 
			{
			return 1;
			}


		/* if ( inclusion_test(TESTLIST,ADD,VSet,0) == 1 )
			{
			return 1;
			}*/
		}
	return 0;
	}

def radical_equality(A,B,VSet,Ord)
	{
	/* This function will be used for checking the termination.	*/

	NA=length(A);

	Ord1=[[0,1],[0,length(VSet)]];    
	Vt=newt;

	for (I=0;I<NA;I++)
		{
		NewB=cons(1-newt*A[I],B);
	        GB = dp_gr_f_main(NewB,cons(Vt,VSet),0,Ord1);
		K=elimination(GB,VSet);

		if ( vars(K) !=[] )
			{
			return 0;
			}

		}

	return 1;
	}

def primedecomposition(P,VSet,Ord,COUNTER,Strategy)
	{
	/* DIVLIST = the set of all computed candidates for iredundant divisors */
	/* INTIDEAL = the intersection of all computed candidates	 	*/
	/* ORIGINAL = a Groebner basis of the input ideal <P> w.r.t. Ord	*/
	/* ORIGINALDIMENSION = the dimension of the input ideal <P>		*/

	RP=P;
	GP=dp_gr_f_main(RP,VSet,Hom,Ord);

	/* First we compute the MSI and the dimension of the 	*/
	/* ideal <P> in K[Vset], where K is the ground field.	*/

	Dimeset=idealdimension(GP,VSet,Ord);
	Dimension=Dimeset[0]; 
	MSI=Dimeset[1];

	if ( dp_gr_print() ) {
		print("The dimension of the ideal is ",2); print(Dimension,2);
		print(".");
	}
	TargetVSet=setminus(VSet,MSI);
	NewGP=dp_gr_f_main(GP,TargetVSet,Hom,Ord);

	if ( vars(NewGP) == [] )
		{
		return [];
		}

	/* Then the ideal is 0-dimension in K[TargetVSet].	*/

	if ( dp_gr_print() ) {
		print("We enter Zero-dimension Prime Decomposition. ",2);
	}

	QP=zeroprimedecomposition(NewGP,TargetVSet,VSet);

	ANS=[];
	NQP=length(QP);

	if ( dp_gr_print() ) {
		print("The number of the newly found component is ",2);
		print(NQP,2);print(". ",2);
	}
	for (I=0;I<NQP;I++)
		{
		ZPrimeideal=QP[I];
		Primedivisor=ZPrimeideal;
		CHECKADD=checkadd2a(Primedivisor,Dimension,DIVLIST,VSet);
		if (CHECKADD == 0 )
			{
			DIVLIST=append(DIVLIST,[[Primedivisor,Dimension]]);
            if ( First_Component ) return 0;
			if (Strategy != 1 )
				{
				/* NO-OPERATION */
				}
			else if ( COUNTER == 0 && Dimension == 0 && ORIGINALDIMENSION == 0 )
				{
				/* NO-OPERATION */
				}
			else if ( INTIDEAL == [] )
				{
				INTIDEAL=Primedivisor;
				}
			else 
				{
				INTIDEAL=ideal_intersection_sfrat(Primedivisor,INTIDEAL,VSet);
				}
			}
		}

	/* We compute prime decomposition for remaining ideal 	*/
	/* I+<ExtPoly> by recursive call.			*/

	if ( Strategy == 1 )
		{
		CHECK=radical_equality(INTIDEAL,ORIGINAL,VSet,Ord);
	
		if (CHECK==1)
			{
			if ( dp_gr_print() ) {
				print("We already obtain all divisor. ");
			}
			STOP = 1;
			return 0;
			}
		}


	if ( Dimension == 0 )
		{
		return 0;
		}

	Ord1 = [[0,length(TargetVSet)],[0,length(MSI)]];
	GP1 = dp_gr_f_main(GP,append(TargetVSet,MSI),Hom,Ord1);
	ExtpolyFactor=extcont_factor(GP1,TargetVSet,0);

	COUNTER=COUNTER+1;

	for ( Tmp = ExtpolyFactor; Tmp != []; Tmp = cdr(Tmp) ) 
		{
		AddPoly=car(Tmp);

		NewGP=cons(AddPoly,GP);

		GP1 = dp_gr_f_main(cons(AddPoly,GP),VSet,Hom,Ord);

		if ( Strategy == 1 )
			{
			CHECKADD=radical_equality(INTIDEAL,NewGP,VSet,Ord);

			if ( CHECKADD != 0 )
				{
				if ( dp_gr_print() ) {
					print("Avoid unnecessary computation. ",2);
				}
				continue;
				}
			}


		DimI=idealdimension(GP1,VSet,Ord)[0];

		if ( REM[ORIGINALDIMENSION-DimI] !=[] )
			{
			REM[ORIGINALDIMENSION-DimI]=cons(GP1,REM[ORIGINALDIMENSION-DimI]);
			}
		else
			{
			REM[ORIGINALDIMENSION-DimI]=[GP1];
			}

		/* primedecomposition(GP1,VSet,Ord,COUNTER,Strategy);
		if ( STOP == 1)
			{
			return 0;
			}*/
		}
	return 0;
	}

/* returns 1 if Poly is a subset of Id(GB) */
/* we assume GB is a gb wrt (V,Ord)        */

def inclusion_test(Poly,GB,V,Ord)
{
	Len = length(GB);
	PS = newvect(Len);
	dp_ord(Ord);
	for ( J = 0, T = GB; T != []; T = cdr(T), J++ )
		PS[J] = dp_ptod(car(T),V);
	for ( J = Len-1, Ind = []; J >= 0; J-- )
		Ind = cons(J,Ind);

	for ( T = Poly; T != []; T = cdr(T) )
		if ( nf_sfrat(Ind,dp_ptod(car(T),V),1,PS)[0] )
			return 0;
	return 1;
}


def zeroprimedecomposition(P,TargetVSet,VSet)
	{
	/* We compute prime decomposition for an ideal <P>	*/
	/* such that <P> is 0-dimensional in K[TargetVSet].	*/

	PD=partial_decomp(P,TargetVSet);

	/* each ideal in PD is an intermediate ideal,		*/
	/* i.e., all minimal polynomials are irreducible.	*/
	/* Thus, each ideal is very likely a prime ideal.	*/
	/* PD=[[P,MP],...], where P is a GB w.r.t. DRL and 	*/
	/* MP is the set of all minimal polynomials.		*/

	NPD=length(PD);
	ANS=[];

	for (I=0;I<NPD;I++)
		{
		COMP=PD[I];

		/* check if the ideal has a variable in generic position*/
	
		CHECK=checkgeneric(COMP,TargetVSet,VSet);

		/* If there is a variable in generic position, then we apply a special 	*/
		/* procedure (zerogenericprimedecomposition). Otherwise we apply a 	*/
		/* general procedure (zeroseparableprimedecomposition).			*/	
		/*									*/
		/* CHECK=[1,M], where M is the minimal polynomail of a variable x	*/
		/* such that x is in generic position, if there is such a variable x.	*/
		/* Otherwise, CHECK=0.							*/

		if (CHECK != 0 )
			{
			if ( TargetVSet != VSet )
				{
				PDiv=contraction(COMP[0],TargetVSet,VSet)[0];
				}
			else
				{
				PDiv=COMP[0];
				}

			ZDecomp=[PDiv];
		
			if ( dp_gr_print() ) {
				print("An intermediate ideal is of generic type. "); 
			}
			}
		else
			{
			if ( dp_gr_print() ) {
				print("An intermediate ideal is not of generic type. ",2); 
			}

			/* We compute the separable closure of <P> by using minimal polynomails.*/
			/* separableclosure outputs 						*/
			/* [a GB Q of separable closure, the exponent vector].			*/
			/* If <P> is already separable, then the exponent vector is 0.		*/

			Sep=separableclosure(COMP,TargetVSet,VSet);

			/* Then, we compute prime divisors of the separable closure by using	*/
			/* generic position.							*/
			/* If Sep[1]=0, COMP[0] is already separable ideal, and hence,		*/
			/* we do not apply radicalideal to its divisors.			*/
	
			if ( Sep[1] != 0 )
				{
				if ( dp_gr_print() ) {
					print("The ideal is inseparable. ",2); 
				}
				CHECK2=checkgeneric2(Sep[2]);
				}	
			else
				{
				if ( dp_gr_print() ) {
					print("The ideal is already separable. ",2); 
				}
				}	

			if ( Sep[1] !=0 && CHECK2 == 1 )
				{
				if ( dp_gr_print() ) {
					print("The separable closure is of generic type. ",2); 
					print("So, the intermediate ideal is prime or primary. ",2); 
				}
				PDiv=convertdivisor(Sep[0],TargetVSet,VSet,Sep[1]);				
				if ( TargetVSet != VSet )
					{
					PDiv=contraction(PDiv,TargetVSet,VSet)[0]; 
					}
				Divisor=radicalideal(PDiv,2,VSet);
				ZDecomp=[Divisor];
				}
			else
				{
#if 0
				ZSDecomp=zeroseparableprimedecomposition(Sep[0],TargetVSet,VSet);
#else
				ZSDecomp=zerosepdec(Sep[0],TargetVSet,VSet,Ord);
#endif
				/* We convert a prime separable ideal in K[TargetVSet] to its 		*/
				/* corresponding prime ideal in K[VSet].				*/

				NComp=length(ZSDecomp);
				ZDecomp=[];

				for (J=0;J<NComp;J++)
					{
					SDiv=ZSDecomp[J];

				/* First we convert a prime separable ideal in K[TargetVSet] to its	*/
				/* corresponding prime/primary ideal in K[TargetVSet] by computing the	*/
				/* image of Frobenius map.						*/

					PDiv=convertdivisor(SDiv,TargetVSet,VSet,Sep[1]);

				/* Then we compute the true prime divisor by contraction and radical	*/
				/* computation.								*/
				
					if ( TargetVSet != VSet )
						{
						PDiv=contraction(PDiv,TargetVSet,VSet)[0]; 
						}

					if (Sep[1] != 0 )
						{
						Divisor=radicalideal(PDiv,2,VSet);
						}
					else
						{
						Divisor=PDiv;
						}

	
					ZDecomp=append(ZDecomp,[Divisor]);
					}
				}
			}
						
		ANS=append(ANS, ZDecomp);
		}
	return map(monic_hc,ANS,VSet);
	}

def monic_hc(Poly,V)
{
	if ( type(Poly) == 4 )
		map(monic_hc,Poly,V);
	else {
		T = dp_ptod(Poly,V);
		return dp_dtop(T/dp_hc(T),V);
	}
}

def zeroseparableprimedecomposition(P,TargetVSet,VSet)
	{
	/* We compute prime decomposition for a separable	*/
	/* 0-dimensional ideal <P> in K[TargetVSet] by using	*/
	/* generic position.					*/

	ANS=[];
	Ord=0;

	NVSet=length(TargetVSet);

	/* The P is already a GB w.r.t. DRL. So, the following	*/
	/* must be replaced with  NewP=P. 			*/

	/* NewGP=dp_gr_f_main(P,TargetVSet,Hom,Ord); */
	NewGP = P;

	/* We search a polynomial f in generic position.	*/
	/* Generic=[f, minimal polynomial of f in newt, newt], 	*/
	/* where newt (X) is a newly introduced variable.	*/

	if ( dp_gr_print() ) {
		print("We search for a linear sum of variables in generic position. ",2);
	}
	Generic=findgeneric(NewGP,TargetVSet,VSet);

	X=Generic[2]; /* newly introduced variable		*/
	Factors=sffctr(Generic[1]);
	NFactors=length(Factors);
	/* irreducible case */
	if ( NFactors == 2 && Factors[1][1] == 1 )
		return [NewGP];
#if 0
	for (J=1;J<NFactors;J++)
		{
		AddPoly=Factors[J][0];
		AddPoly2=X-Generic[0]; 

		TestGP=append(NewGP,[AddPoly, AddPoly2]);
		VS=cons(X,TargetVSet);
		Q=dp_gr_f_main(TestGP,VS,Hom,[[0,1],[Ord,NVSet]]); 
		QR=elimination(Q,VSet);
		ANS=append(ANS,[QR]);
		}
#else
	/* noro */
	dp_ord([[0,1],[Ord,NVSet]]); 
	XVSet = cons(X,TargetVSet);
	PS = newvect(length(NewGP)+1,map(dp_ptod,cons(X-Generic[0],NewGP),XVSet));
	for ( I = length(NewGP), Ind = [];I >= 0; I-- )
		Ind = cons(I,Ind);
	Ind = reverse(Ind);
	YSet = setminus(VSet,TargetVSet);
	NYSet = length(YSet);
	ElimVSet = append(TargetVSet,YSet);
	ElimOrd = [[Ord,NVSet],[0,NYSet]];
	for ( J = 1; J < NFactors; J++ ) {
		dp_ord([[0,1],[Ord,NVSet]]); 
		Factor = dp_dtop(nf_sfrat(Ind,dp_ptod(Factors[J][0],XVSet),1,PS)[0],XVSet);
#if 0
		Q = dp_gr_f_main(cons(Factor,NewGP),ElimVSet,Hom,ElimOrd);
#else
        if ( NDGR ) {
		  Q0 = nd_gr(cons(Factor,NewGP),ElimVSet,-1,0);
		  Q = nd_gr(Q0,ElimVSet,-1,ElimOrd);
		  Q = nd_gr(Q,TargetVSet,-1,Ord);
          dp_ord(0);
        } else {
		  Q0 = dp_gr_f_main(cons(Factor,NewGP),ElimVSet,Hom,0);
		  Q = dp_gr_f_main(Q0,ElimVSet,Hom,ElimOrd);
		  Q = dp_gr_f_main(Q,TargetVSet,Hom,Ord);
        }
#endif
		ANS = cons(Q,ANS);
	}
#endif
	return ANS;
	}

/* partial decomposition by sums of variables 
   -> zeroseparableprimedecomposition */

#if 0
def zerosepdec(P,W,V,Ord)
{
	/* P is GB wrt (W,Ord) */	
	dp_ord(Ord);
	DIM = length(dp_mbase(map(dp_ptod,P,W)));
	/* preprocessing */
	N = length(W);
	C = newvect(N);
	C[0] = 1;
	Vt = ttttt;
	do {			
		for ( I = 0, LF = 0; I < N; I++ )
			if ( C[I] ) LF += W[I];
		LF = simp_ff(LF);
		MP = minipoly_sf(P,W,Ord,LF,Vt);
		Factors = map(first,cdr(sffctr(MP)));
		if ( deg(MP,Vt) == DIM ) {
			if ( length(Factors) == 1 )
				return [P];
			else
				return zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));		
		} else if ( length(Factors) > 1 ) {
			R =  zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));		
			S = [];
			for ( TR = R; TR != []; TR = cdr(TR) )
				S = append(zerosepdec(car(TR),W,V,Ord),S);	
			return S;
		}
	} while ( !nextchoice(C) );
	/* we could not find any useful combination */
	R = zeroseparableprimedecomposition(P,W,V);
	return R;
}
#else

/* we only search for useful poly among v[i]+v[j] */
def zerosepdec(P,W,V,Ord)
{
	/* P is GB wrt (W,Ord) */	
	dp_ord(Ord);
	DIM = length(dp_mbase(map(dp_ptod,P,W)));
	/* preprocessing */
	N = length(W);
	C = newvect(N);
	C[0] = 1;
	Vt = ttttt;
	for ( I1 = 0; I1 < N; I1++ )
		for ( I2 = I1+1; I2 < N; I2++ ) {
			LF = simp_ff(W[I1]+W[I2]);
			MP = minipoly_sf(P,W,Ord,LF,Vt);
			Factors = map(first,cdr(sffctr(MP)));
			if ( deg(MP,Vt) == DIM ) {
				if ( length(Factors) == 1 )
					return [P];
				else
					return zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));		
			} else if ( length(Factors) > 1 ) {
				R =  zerosepdec_main(P,W,V,Ord,map(subst,Factors,Vt,LF));		
				S = [];
				for ( TR = R; TR != []; TR = cdr(TR) )
					S = append(zerosepdec(car(TR),W,V,Ord),S);	
				return S;
			}
		}
	/* we could not find any useful combination */
	R = zeroseparableprimedecomposition(P,W,V);
	return R;
}
#endif

def zerosepdec_main(P,W,V,Ord,Factors)
{
	dp_ord(Ord); 
	N = length(P);
	NFactors = length(Factors);
	PS = newvect(length(P),map(dp_ptod,P,W));
	for ( I = length(P)-1, Ind = [];I >= 0; I-- )
		Ind = cons(I,Ind);
	Y= setminus(V,W);
	NW = length(W);
	NY = length(Y);
	ElimV = append(W,Y);
	ElimOrd = [[Ord,NW],[0,NY]];
	ANS = [];
	for ( J = 0; J < NFactors; J++ ) {
		Factor = dp_dtop(nf_sfrat(Ind,dp_ptod(Factors[J],W),1,PS)[0],W);
        if ( NDGR ) {
		  Q = nd_gr(cons(Factor,P),ElimV,-1,ElimOrd);
		  Q = nd_gr(Q,W,-1,Ord);
          dp_ord(0);
        } else {
		  Q = dp_gr_f_main(cons(Factor,P),ElimV,Hom,ElimOrd);
		  Q = dp_gr_f_main(Q,W,Hom,Ord);
        }
		ANS = cons(Q,ANS);
        if ( First_Component ) break;
	}
	return ANS;
}

def nextchoice(C)
{
	N = size(C)[0];
	for ( I = N-1, Carry = 1; I >= 0; I-- ) {
		if ( C[I] ) {
			C[I] = !Carry;
			Carry = !C[I];
		} else {
			C[I] = Carry;
			Carry = 0;
		}
	}
	return Carry;
}

def separableclosure(CP,TargetVSet,VSet)
	{
	/* We compute a separable ideal <Q> from <P> by 	*/
	/* computing inverse Frobenius map.			*/

	Ord=0;

	NVSet=length(TargetVSet);
	IndexVector=newvect(NVSet);

	/* First we check if <P> is already separable.		*/
	/* CHECK == 1 if <P> is already separable,		*/
	/* Otherwise, CHECK is the list of pairs of		*/
	/* the separable closure sc(m) of the minimal polynomial*/
	/* m and the exponent e, i.e. m=sc(m)(t^(q^e)).		*/

	CHECK=checkseparable(CP,TargetVSet,VSet);

	if ( CHECK == 1 )
		{
		if ( dp_gr_print() ) {
			print("This is already a separable ideal.", 2);
		}
		return [CP[0],0];
		}

	if ( dp_gr_print() ) {
		print("This is not a separable ideal, so we make its separable closure.", 2);
	}
	WSet=makecounterpart(TargetVSet);
	Char=characteristic_ff();

	NewP=CP[0];
	EXPVECTOR=newvect(NVSet);
	CHECKEXPVECTOR=newvect(NVSet);

	for (I=0;I<NVSet;I++)
		{
		EXPVECTOR[I]=CHECK[NVSet-I-1][0];
		CHECKEXPVECTOR[I]=deg(CHECK[NVSet-I-1][1],TargetVSet[I]);

		POW=Char^CHECK[NVSet-I-1][0];

		AddPoly=TargetVSet[I]^POW-WSet[I];
		/*SepClosure=subst(CHECK[I][0],TargetVSet[I],WSet[I]);*/

		NewP=cons(AddPoly,NewP);
		}

	NewOrder=[[Ord,NVSet],[Ord,NVSet]];
	XSet=append(TargetVSet,WSet);
	YSet=append(VSet,WSet);

	NewG=dp_gr_f_main(NewP,XSet,Hom,NewOrder);

	NWSet=setminus(YSet,TargetVSet);
	
	Q=elimination(NewG,NWSet);

	NQ=length(Q);

	ANS=[];
	for (I=NQ-1;I>=0;I--)
		{
		Poly=Q[I];
		for (J=0;J<NVSet;J++)
			{
			Poly=subst(Poly,WSet[J],TargetVSet[J]);
			}
		ANS=cons(Poly,ANS);
		}
	return [ANS,EXPVECTOR,CHECKEXPVECTOR];
	}

def convertdivisor(P,TargetVSet,VSet,ExVector)
	{
	/* We compute a corresponding ideal <Q> from <P> by 	*/
	/* computing Frobenius map.				*/

	if (ExVector == 0 )
		{
		return P;
		}

	NVSet=length(TargetVSet);
	WSet=makecounterpart(TargetVSet);
	Char=characteristic_ff();
	Ord=0;

	NewP=P;

	for (I=0;I<NVSet;I++)
		{
		POW=Char^ExVector[I];

		AddPoly=TargetVSet[I]-WSet[I]^POW;
		/*SepClosure=subst(CHECK[I][0],TargetVSet[I],WSet[I]);*/

		NewP=cons(AddPoly,NewP);
		}

	NewOrder=[[Ord,NVSet],[Ord,NVSet]];

	XSet=append(TargetVSet,WSet);
	YSet=append(VSet,WSet);

	NewG=dp_gr_f_main(NewP,XSet,Hom,NewOrder);
	
	NWSet=setminus(YSet,TargetVSet);

	Q=elimination(NewG,NWSet);

	NQ=length(Q);

	ANS=[];
	for (I=NQ-1;I>=0;I--)
		{
		Poly=Q[I];
		for (J=0;J<NVSet;J++)
			{
			Poly=subst(Poly,WSet[J],TargetVSet[J]);
			}
		ANS=cons(Poly,ANS);
		}
	return ANS;
	}


def findgeneric(P,TargetVSet,VSet)
	{
	/* The number of trials is set as Trails. */
	NTargetVSet=length(TargetVSet);
	MNumber=0;
	/* Trials=100;*/

	if ( Trials == 0 )
		Trials = 2;

	Lineardimension=lineardimension(P,TargetVSet);
	NewVar=newt;
	Count=0;  
	Ord=0;

#if 0
	while(Count < Trials )
		{
		Candidate=0;
		MNumber=MNumber+1;

		MAGIC=makemagic(TargetVSet,MNumber);

		for (I=0;I<NTargetVSet;I++)
			{
			Candidate=Candidate+MAGIC[I]*TargetVSet[I];
			}
		MinPoly=minipoly_sf(P,TargetVSet,Ord,Candidate,NewVar);
		Deg=deg(MinPoly,NewVar);
		if ( Deg == Lineardimension )
			{
			return [Candidate,MinPoly,NewVar];
			}
		Count=Count+1;
		}
#else
	YSet = setminus(VSet,TargetVSet);
	NYSet = length(YSet);
	Eval = newvect(NYSet);
	Q1 = field_order_ff()-1;
	HM = hmlist(P,TargetVSet,Ord);
	for ( Count = 0; Count < Trials; Count++ ) {
		Candidate = ptosfp(2)*TargetVSet[0];
		Candidate = 0;
		for ( I = 0; I < NTargetVSet; I++ )
			Candidate += random_ff()*TargetVSet[I];
		do {
			for ( I = 0; I < NYSet; I++ )
				Eval[I] = random()%Q1;
		} while ( !valid_modulus_sfrat(HM,YSet,Eval) );
		P0 = map(eval_sfrat,P,YSet,Eval);
		MinPoly0 = minipoly_sf(P0,TargetVSet,Ord,Candidate,NewVar);
		Deg = deg(MinPoly0,NewVar);
		if ( Deg == Lineardimension ) {
			MinPoly=minipoly_sf(P,TargetVSet,Ord,Candidate,NewVar);
			if ( deg(MinPoly,NewVar) != Deg )
				error("findgeneric : cannot happen");
			return [Candidate,MinPoly,NewVar];
		}
	}
#endif
	if ( dp_gr_print() ) {
		print("Extend the ground field. ",2);
	}
	error();
	}

def makemagic(VSet,MNumber)
	{
	NVSet=length(VSet);
	MAGIC=[1];
	for (I=1;I<NVSet;I++)
		{
		U=ptosfp(MNumber^I);
		MAGIC=append(MAGIC,[U]);
		}
	return MAGIC;
	}

#if 0
def zerogenericprimedecomposition(CP,MinP,TargetVSet,VSet)
	{
	P=CP[0];

	FACTORS=sffctr(MinP);
	NFACTORS=length(FACTORS);

	if ( NFACTORS == 2 )
		{
		return [P];
		}
	
	ANS=[];

	for (I=1;I<NFACTORS;I++)
		{
		AddPoly=NFACTORS[I][0];
		NewP=cos(AddPoly,P);
		NewG=dp_gr_f_main(NewP,TargetVSet,Hom,Ord);
		ANS=cons(NewG,ANS);
		}
	return ANS;
	}
#endif 


def lineardimension(P,VSet)
	{
	Ord=0;

	PP=dp_gr_f_main(P,VSet,Hom,Ord);
	dp_ord(Ord);
	Dimension=length(dp_mbase(map(dp_ptod,PP,VSet)));
	return Dimension;
	}

def checkgeneric(CP,TargetVSet,VSet)
	{
	Dimension=lineardimension(CP[0],TargetVSet);
	NVSet=length(TargetVSet);
	Flag=0;

	for (I=0;I<NVSet;I++)
		{
		MinP=CP[1][I];

		if ( deg(MinP,TargetVSet[NVSet-I-1]) == Dimension )
			{
			return [1,MinP];
			}
		}
	return 0;
	}

def checkseparable(CP,TargetVSet,VSet)
	{
	NVSet=length(TargetVSet);
	EXPVector=newvect(NVSet);
	Flag=1;

	for (I=0;I<NVSet;I++)
		{
		MinP=CP[1][I];
		CHECK=checkseparablepoly(MinP,TargetVSet[NVSet-I-1]);
		if ( CHECK[0] != 0 )
			{
			Flag = 0; 
			}
		EXPVector[I]=CHECK;
		}
	if ( Flag == 0 )
		{
		return EXPVector;
		}
	return 1;
	}

def checkseparablepoly(F,V)
	{
	P = characteristic_ff();
	E = 0;
	Vt = newt;
	while ( 1 ) {
		if ( diff(F,V) != 0 )
			return [E,F];
		T = srem(F,V^P-Vt,V);
		F = subst(T,Vt,V);
		E++;
		}
	}

def extcont(P,V,Ord)
	{
	/* We assume P is a Groebner Basis w.r.t. Ord. 		*/
	/* We compute a polynomial f such that 			*/
	/* \sqrt{<P>}=\sqrt{<P>+<f>}\cap \sqrt{(<P>)^{ec}}	*/			
	/* by B.W. Proposition 8.96.				*/
	/* Remark: The square free part of f is also OK.	*/

	dp_ord(Ord); 
	HC = map(dp_hc,map(dp_ptod,P,V));
	LCM = car(HC);
	for ( T = cdr(HC); T != []; T = cdr(T) )
		LCM = lcm_sfrat(LCM,car(T));
	return LCM;
	}

def first(A)
{
	return A[0];
}

def set_union(A,B)
{
	for ( T = A; T != []; T = cdr(T) )
		B = insert_element(car(T),B);
	return B;
}

def insert_element(E,A)
{
	for ( T = A; T != []; T = cdr(T) )
		if ( E == car(T) )
			break;
	if ( T == [] )
		return cons(E,A);
	else
		return A;
}

def extcont_factor(P,V,Ord)
	{
	/* We assume P is a Groebner Basis w.r.t. Ord. 		*/
	/* We compute a polynomial f such that 			*/
	/* \sqrt{<P>}=\sqrt{<P>+<f>}\cap \sqrt{(<P>)^{ec}}	*/			
	/* by B.W. Proposition 8.96.				*/
	/* Remark: The square free part of f is also OK.	*/

	dp_ord(Ord); 
	HC = map(dp_hc,map(dp_ptod,P,V));
	F = [];
	for ( T = HC; T != []; T = cdr(T) ) {
		F1 = map(first,cdr(sffctr(car(T))));
		F = set_union(F1,F);
	}
	return F;
	}

def contraction(P,V,W)
	{
	/* We assume P is a Groebner Basis w.r.t. Ord.		*/
	/* We compute the contraction <P>^{c} by 		*/
	/*  <P>^c=(<G>:f^\infty)				*/			
	/* by B.W. Proposition 8.91.				*/
	/* This procedure is called by zeroprimedecomposition.	*/
	/* So, P is supposed to be a GB w.r.t. DRL.		*/

	Ord0 = dp_ord();
	Ord=0;
	YSet=setminus(W,V);

	Ord1 = [[Ord,length(V)],[0,length(YSet)]];
	W1 = append(V,YSet);
	GP1 = dp_gr_f_main(P,W1,Hom,Ord1);

	Factor = extcont_factor(GP1,V,Ord);
	for ( F = 1, T = Factor; T != []; T = cdr(T) )
		F *= car(T);
	Vt = newt;	

	G = dp_gr_f_main(cons(1-Vt*F,P),cons(Vt,W),0,[[0,1],[Ord,length(W)]]);
	R = [];
	for ( T = G; T != []; T = cdr(T) )
		if ( !member(Vt,vars(car(T))) )
			R = cons(car(T),R);
	dp_ord(Ord0);
	return [R,F];
	}

def checkgeneric2(LIST)
	{
	NList=size(LIST)[0];

	FLAG=0;

	for (I=0;I<NList;I++)
		{
		if ( LIST[I] > 1 )
			{
			FLAG=FLAG+1;
			}
		}

	if (FLAG < 2 )
		{
		return 1;
		}
	return 0;
	}

#if 0 
def checkseparablepoly(P,V)
	{
	TestP=P;
	CHECK=diff(TestP,V); 
	Count=0;

	while ( CHECK !=0 )
		{
		if  ( deg(TestP,V) != 0 )
			{
			break;
			}
		TestP=pdivide(TestP);
		CHECK=diff(TestP,V);
		Count=Count+1;
		}

	return [TestP,Count];
	}

def pdivide(F,V)
	{
	Char=characteristic_ff();
	TestP=P;

	Deg=ideg(TestP,V);
	DegP=idiv(Deg,Char); 

	if ( irem(Deg,Char) != 0 )
		{
		error;
		}
	ANS=0; 

	for (I=O;I<DegP;I++)
		{
		TempDeg=I*Char;
		ANS=ANS+coeff(TestP,V,TempDeg)*X^I;
		}
	return ANS;
	}
#endif 


def convsf(PP,VSet,Ord,Flag)
	{
	/* Flag = 0 or 1 				*/
	/* If Flag = 1, we compute the intersection 	*/
	/* of the Galois orbit.				*/

	CHECK=checkgaloisorbit(PP,VSet,Ord,Flag);

	NewPP=CHECK[0];
	
	ANS=[];

	NPP=length(NewPP);

	for (I=0;I<NPP;I++)
		{
		NewComp=convertsmallfield(CHECK[Flag][I],VSet,Ord);		
		ANS=cons(NewComp,ANS);
		}
	return ANS;
	}

def convertsmallfield(PP,VSet,Ord)
	{
	dp_ord(Ord);
	NVSet=length(VSet);
	Char=characteristic_ff();
	ExtDeg=extdeg_ff();

	NewV=pgpgpgpgpgpgpg;	
	MPP=map(monic_hc,PP,VSet);
	MPP=map(sfptopsfp,MPP,NewV);

	DefPoly=setmod_ff()[1];
	/* GF(p) case */
	if ( !DefPoly )
		return MPP;

	MinPoly=subst(DefPoly,var(DefPoly),NewV);
	XSet=cons(NewV,VSet);

	Ord1=[[0,1],[Ord,NVSet]];

	/* setmod_ff(Char,1);*/

	NewP=dp_gr_mod_main(cons(MinPoly,MPP),XSet,0,Char,Ord1);

	ANS=elimination(NewP,VSet);
	
	/* setmod_ff(Char,ExtDeg);*/

	return ANS;
	}

def checkgaloisorbit(PP,VSet,Ord,Flag)
	{
	NPP=length(PP);
	TmpPP=PP;
	ExtDeg=extdeg_ff();
	
	ANS=[];
	BNS=[];

	while (TmpPP != [] )
		{
		TestP=car(TmpPP);
		Dim=TestP[1];
		TmpPP=cdr(TmpPP);
		NewP=TestP[0];

		ANS=cons(TestP[0],ANS);
		
		for (J=1;J<ExtDeg;J++)
			{
			G1=map(sf_galois_action,TestP[0],J);		

			if ( setminus(G1,TestP[0]) == [] )
				{
				break;
				}
			if (Flag == 1 )
				{
				NewP=ideal_intersection_sfrat(NewP,G1,VSet);
				}
			TmpPP=deletecomponent(TmpPP,[G1,Dim]);
			}
		BNS=cons(NewP,BNS);

		}
	return [ANS,BNS];
	}

def deletecomponent(PP,G)
	{
	TmpPP=PP;

	while( TmpPP !=[] )
		{
		Test=car(TmpPP)[0];

		if ( setminus(Test,G[0]) == [] )
			{
			ANS=setminus(PP,[G]);
			return ANS;
			}

		TmpPP=cdr(TmpPP);
		}
	error();
	}


def pfctr(F)
{
	if ( type(F) == 1 )
		return [[F,1]];
	P = characteristic_ff();
	E = extdeg_ff();
	F = sfptop(F);
	check_coef(F,P);
	D = modfctr(F,P);
	setmod_ff(P,E);
	return map(mapsf_first,D);
}

def check_coef(F,P)
{
	V = vars(F);
	D = dp_ptod(F,V);
	for ( T = D; T; T = dp_rest(T) )
		if ( dp_hc(T) >= P )
			error("invalid coef");
}

def mapsf_first(D)
{
	return [ptosfp(D[0]),D[1]];
}

def partial_decomp(B,V)
{
	T0 = time();
	if ( ParallelMinipoly ) {
		map(ox_cmo_rpc,ParallelMinipoly,"setmod_ff",characteristic_ff(),extdeg_ff());
		map(ox_pop_cmo,ParallelMinipoly);
	}
	B = map(simp_ff,B);
	B = dp_gr_f_main(B,V,0,0);
	R = partial_decomp0(B,V,length(V)-1);
	if ( PartialDecompByLex ) {
		R0 = [];
		for ( TR = R; TR != []; TR = cdr(TR) ) {
			T = car(TR);
			S = dp_gr_f_main(T[0],V,0,0);
			R0 = cons([S,T[1]],R0);
		}
		R = reverse(R0);
	}
	T_PD += time()[3]-T0[3];
	return R;
}

def setintersection(A,B)
{ 
	for ( L = []; A != []; A = cdr(A) )
		if ( member(car(A),B) )
			L = cons(car(A),L);
	return L;   
}

/* returns [[Plist,Mlist],...] */

def partial_decomp0(B,V,I)
{
	N = length(V);
	if ( I < 0 )
		return [[B,[]]];
	Ord = PartialDecompByLex ? [[0,I],[2,N-I]] : 0;
#if 0
	if ( setminus(vars(B),V) == [] )
		B = fglm_sf_0dim(B,V,Ord0,Ord);
	else
#endif
		B = dp_gr_f_main(B,V,0,Ord);
	if ( type(B[0]) == 1 )
		return [];
	if ( !zero_dim(B,V,Ord) )
		error("non zero-dimensional ideal");
	/* XXX */
	Vt = ttttt;
	VI = V[I];
	if ( PartialDecompByLex ) {
		for ( J = 0, W = V; J < I; J++ )
			W = cdr(W);
		VW = setminus(V,W);
		for ( Bw = [], T = B; T != []; T = cdr(T) )
			if ( setintersection(vars(car(T)),VW) == [] )
				Bw = cons(car(T),Bw);
		MI = minipoly_sf(Bw,W,2,VI,Vt);
	} else
		MI = minipoly_sf(B,V,0,VI,Vt);
	MIF = sffctr(MI);
	/* if MI is irreducible, then process the next variable */
	if ( length(MIF) == 2 && MIF[1][1] == 1 ) {
		L = partial_decomp0(B,V,I-1);
		R = [];
		for ( S = L; S != []; S = cdr(S) ) {
			L = car(S);
			R = cons([L[0],cons(subst(MI,Vt,VI),L[1])],R);
		}
		return R;
	}

	/* for nf computation */
	Ord1 = PartialDecompByLex ? 2 : [[0,1],[0,N]];
	Len = length(B);
	PS = newvect(Len+1);
	dp_ord(Ord1);
	XV = cons(Vt,V);
	for ( J = 0, T = cons(Vt-VI,B); T != []; T = cdr(T), J++ )
		PS[J] = dp_ptod(car(T),XV);
	for ( J = 0, GI = []; J <= Len; J++ )
		GI = cons(J,GI);
	if ( PartialDecompByLex )
		GI = reverse(GI);
	
	R = [];
	for ( T = MIF; T != []; T = cdr(T) ) {
		Mt = car(car(T));
		if ( type(Mt) == 1 )
			continue;
		dp_ord(Ord1);
		NfMt = dp_dtop(nf_sfrat(GI,dp_ptod(Mt,XV),1,PS)[0],XV);
		if ( NfMt )
			B1 = dp_gr_f_main(cons(NfMt,B),V,0,Ord);
		else
			B1 = B;
		Rt = partial_decomp0(B1,V,I-1);
		for ( S = Rt; S != []; S = cdr(S) ) {
			L = car(S);
			R = cons([L[0],cons(subst(Mt,Vt,VI),L[1])],R);
		}
	}
	return R;
}


/* G is a gb wrt (V,O) */

def minipoly_sf(G,V,O,F,V0)
{
	T0 = time();
	Vc = cons(V0,setminus(vars(G),V));
	if ( ParallelMinipoly ) {
		Proc0 = ParallelMinipoly[0];
		Proc1 = ParallelMinipoly[1];
		if ( length(Vc) == 1 ) {
			ox_rpc(Proc0,"minipoly_sf_by_buchberger",G,V,O,F,V0,1);
			ox_rpc(Proc1,"minipoly_sf_0dim",G,V,O,F,V0,1);
			map(ox_get,ParallelMinipoly);
			/* 258=SM_popSerializedLocalObject */
			map(ox_push_cmd,ParallelMinipoly,258);
			F = ox_select(ParallelMinipoly);
			MP = ox_get(F[0]);
			if ( F[0] == Proc0 ) {
				if ( length(F) == 1 )
					B_Win++;
				else
					ox_get(Proc1);
			} else {
				if ( length(F) == 1 )
					D_Win++;
				else
					ox_get(Proc0);
			}
			ox_reset(Proc0);
			ox_reset(Proc1);
		} else if ( length(Vc) == 2 ) {
			ox_rpc(Proc0,"minipoly_sf_by_buchberger",G,V,O,F,V0,1);
			ox_rpc(Proc1,"minipoly_sfrat",G,V,O,F,V0,1);
			map(ox_get,ParallelMinipoly);
			/* 258=SM_popSerializedLocalObject */
			map(ox_push_cmd,ParallelMinipoly,258);
			F = ox_select(ParallelMinipoly);
			MP = ox_get(F[0]);
			if ( F[0] == Proc0 ) {
				if ( length(F) == 1 )
					B_Win++;
				else
					ox_get(Proc1);
			} else {
				if ( length(F) == 1 )
					D_Win++;
				else
					ox_get(Proc0);
			}
			ox_reset(Proc0);
			ox_reset(Proc1);
		} else
			MP = minipoly_sf_by_buchberger(G,V,O,F,V0,0);
	} else if ( BuchbergerMinipoly )
		MP = minipoly_sf_by_buchberger(G,V,O,F,V0,0);
	else {
		if ( length(Vc) == 1 )
			MP = minipoly_sf_0dim(G,V,O,F,V0,0);
		else if ( length(Vc) == 2 )
			MP = minipoly_sfrat(G,V,O,F,V0,0);
		else
			MP = minipoly_sf_by_buchberger(G,V,O,F,V0,0);
	}
	T_MP += time()[3]-T0[3];
	return MP;
}

def minipoly_sf_by_buchberger(G,V,O,F,V0,Server)
{
	if ( Server )
		ox_sync(0);
	Vc = cons(V0,setminus(vars(G),V));
	Gf = cons(simp_ff(V0-F),G);
	Vf = append(V,Vc);
	Gelim = dp_gr_f_main(Gf,Vf,1,[[0,length(V)],[0,length(Vc)]]);
	for ( Gc = [], T = Gelim; T != []; T = cdr(T) ) {
		Vt = setminus(vars(car(T)),Vc);
		if ( Vt == [] )
			Gc = cons(car(T),Gc);
	}
	Gt = dp_gr_f_main(Gc,Vc,1,[[0,1],[0,length(Vc)-1]]);	
	Pmin = car(Gt); Dmin = deg(Pmin,V0); 
	for ( T = cdr(Gt); T != []; T = cdr(T) ) {
		Dt = deg(car(T),V0);
		if ( Dt < Dmin ) {
			Pmin = car(T); Dmin = Dt; 
		}
	}
	Cont = sfcont(Pmin,V0);
	return sdiv(Pmin,Cont);
}

def minipoly_sf_0dim(G,V,O,F,V0,Server)
{
	if ( Server )
		ox_sync(0);

    if ( Minipoly_SBA ) {
       G1 = cons(simp_ff(V0-F),G);
       V1 = append(V,[V0]);
       G2 = nd_sba(G1,V1,-1,[[0,length(V)],[0,1]]|sba_modord=[[],0]);
       dp_ord(0);
       for ( T = G2; T != []; T = cdr(T) )
         if ( vars(car(T)) == [V0] ) break;
       if ( T == [] )
         error("minipoly does not exit");
       return car(T);
    }

	N = length(V);
	Len = length(G);
	dp_ord(O);
	PS = newvect(Len);
	for ( I = 0, T = G; T != []; T = cdr(T), I++ )
		PS[I] = dp_ptod(car(T),V);
	for ( I = Len-1, HL = []; I >= 0; I-- )
		HL = cons(dp_ht(PS[I]),HL);
	for ( I = Len - 1, GI = []; I >= 0; I-- )
		GI = cons(I,GI);
	MB = dp_mbase(HL); DIM = length(MB); UT = newvect(DIM);
	U = dp_ptod(simp_ff(F),V);
	U = dp_nf_f(GI,U,PS,1);
	for ( I = 0; I < DIM; I++ )
		UT[I] = [MB[I],dp_nf_f(GI,U*MB[I],PS,1)];

	T = dp_ptod(simp_ff(1),[V0]);
	TT = dp_ptod(simp_ff(1),V);
	G = H = [[TT,T]];

	for ( I = 1; ; I++ ) {
		if ( dp_gr_print() )
			print(".",2);
		T = dp_ptod(simp_ff(V0^I),[V0]);
		TT = dp_nf_tab_f(H[0][0],UT);
		H = cons([TT,T],H);
		L = dp_lnf_f([TT,T],G);
		if ( !L[0] ) {
			if ( dp_gr_print() )
				print("");
			return dp_dtop(L[1],[V0]); /* XXX */
		} else
			G = insert(G,L);
	}
}

#if 1
def minipoly_sf_rat(G,V,F,V0)
{
	Vc = setminus(vars(G),V);
	Gf = cons(V0-F,G);
	Vf = append(V,[V0]);
	G3 = dp_gr_f_main(map(simp_ff,Gf),Vf,0,3);
	for ( T = G3; T != []; T = cdr(T) ) {
		Vt = setminus(vars(car(T)),Vc);
		if ( Vt == [V0] )
			return car(T);
	}
}

def minipoly_mod(G,V,Mod,F,V0)
{
	Vc = setminus(vars(G),V);
	Gf = cons(V0-F,G);
	Vf = append(V,[V0]);
	G3 = dp_gr_mod_main(Gf,Vf,1,Mod,3);
	for ( T = G3; T != []; T = cdr(T) ) {
		Vt = setminus(vars(car(T)),Vc);
		if ( Vt == [V0] )
			return car(T);
	}
}
#endif

/* find N/D s.t. F = N/D mod V^(2M), deg(N),deg(D)<M */

def pade(F,M)
{
	V = var(F);
	A = newvect(M);
	D = 0;
	for ( I = 0; I < M; I++ ) {
		A[I] = strtov("b"+rtostr(I));
		D += A[I]*V^I;
	}
	T = F*D;
	M2 = M*2;
	R = [];
	for ( I = M; I < M2; I++ )
		R = cons(coef(T,I,V),R);
}

def minipoly_sfrat(G0,V,O,P,V0,Server)
{
	if ( Server )
		ox_sync(0);
	if ( !zero_dim(hmlist(G0,V,O),V,O) )
		error("minipoly_sfrat : ideal is not zero-dimensional!");

	G1 = cons(V0-P,G0);
	if ( type(O) <= 1 )
		O1 = [[0,1],[O,length(V)]];
	else
		O1 = cons([0,1],O);
	V1 = cons(V0,V);
	W = append(V,[V0]);	
	Vc = setminus(vars(G0),V);

	N = length(V1);
	dp_ord(O1);
	HM = hmlist(G1,V1,O1);
	MB = dp_mbase(map(dp_ptod,HM,V1));
	dp_ord(O);

	Nc = length(Vc);
	Eval = newvect(Nc);

	/* heristic lower bound of deg(MP) */
	Q1 = field_order_ff()-1;
	do {
		for ( I = 0; I < Nc; I++ )
			Eval[I] = random()%Q1;
	} while ( !valid_modulus_sfrat(HM,Vc,Eval) );
	G0E = map(eval_sfrat,G0,Vc,Eval);
	P0 = eval_sfrat(P,Vc,Eval);
	MP = minipoly_sf(G0E,V,O,P0,V0);
	DMP = deg(MP,V0);

	for ( I = 0; I < Nc; I++ )
		Eval[I] = 0;
	while ( 1 ) {
		if ( valid_modulus_sfrat(HM,Vc,Eval) ) {
			G0E = map(eval_sfrat,G0,Vc,Eval);
			P0 = eval_sfrat(P,Vc,Eval);
			MP = minipoly_sf(G0E,V,O,P0,V0);
			D = deg(MP,V0);
			if ( D >= DMP ) {
				DMP = D;
				for ( TL = [], MPT = 0, J = 0; J <= D; J++ ) {
					TL = cons(V0^J,TL);
					MPT += car(TL);
				}
				NF = gennf_sfrat(G1,TL,V1,O1,V0,1)[0];
				R = tolex_sfrat_main(V1,O1,NF,[MPT],Vc,Eval,MB);
				if ( R )
					return sdiv(R[0],sfcont(R[0],V0));
			}
		}
		next_eval_sfrat(Eval);
	}
}

def next_eval_sfrat(Eval)
{
	N = size(Eval)[0];
	for ( I = N-1; I >= 0; I-- )
		if ( Eval[I] ) break;
	if ( I < 0 ) Eval[N-1] = 1;
	else if ( I == 0 ) {
		T = Eval[0]; Eval[0] = 0; Eval[N-1] = T+1;
	} else {
		Eval[I-1]++; T = Eval[I];
		for ( J = I; J < N-1; J++ )
			Eval[J] = 0;
		Eval[N-1] = T-1;
	}
}

def eval_sfrat(F,V,Eval)
{
	for ( I = 0; V != []; V = cdr(V), I++ )
		F = subst(F,car(V),ptosfp(Eval[I]));
	return F;
}

def valid_modulus_sfrat(HL,Vc,Eval) {
	for ( T = HL; T != []; T = cdr(T) ) {
		C = car(T);
		for ( S = Vc, I = 0; S != []; S = cdr(S), I++ )
			C = subst(C,car(S),ptosfp(Eval[I]));
		if ( !C )
			break;
	}
	return T == [] ? 1 : 0;
}

def gennf_sfrat(G,TL,V,O,V0,FLAG)
{
	N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
	for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) {
		PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL);
	}
	for ( I = 0, DTL = []; TL != []; TL = cdr(TL) )
		DTL = cons(dp_ptod(car(TL),V),DTL);
	for ( I = Len - 1, GI = []; I >= 0; I-- )
		GI = cons(I,GI);
	T = car(DTL); DTL = cdr(DTL);
	H = [nf_sfrat(GI,T,T,PS)];

	USE_TAB = (FLAG != 0);
	if ( USE_TAB ) {
		T0 = time()[0];
		MB = dp_mbase(HL); DIM = length(MB);
		U = dp_ptod(V0,V);
		UTAB = newvect(DIM);
		for ( I = 0; I < DIM; I++ ) {
			UTAB[I] = [MB[I],remove_cont_sfrat(nf_sfrat(GI,U*MB[I],1,PS))];
			if ( dp_gr_print() )
				print(".",2);
		}
		if ( dp_gr_print() )
			print("");
		TTAB = time()[0]-T0;
	}

	T0 = time()[0];
	for ( LCM = 1; DTL != []; ) {
		if ( dp_gr_print() )
			print(".",2);
		T = car(DTL); DTL = cdr(DTL);
		if ( L = search_redble(T,H) ) {
			DD = dp_subd(T,L[1]);
			if ( USE_TAB && (DD == U) ) {
				NF = nf_tab_sfrat(L[0],UTAB);
				NF = [NF[0],dp_hc(L[1])*NF[1]*T];
			} else
				NF = nf_sfrat(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
		} else
			NF = nf_sfrat(GI,T,T,PS);
		NF = remove_cont_sfrat(NF);
		H = cons(NF,H);
		LCM = lcm_sfrat(LCM,dp_hc(NF[1]));
	}
	TNF = time()[0]-T0;
	if ( dp_gr_print() )
		print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
	return [[map(adj_dn_sfrat,H,LCM),LCM],PS,GI];
}

def lcm_sfrat(A,B)
{
	G = sfgcd(A,B);
	return sdiv(A,G)*B;
}

#define REDCONT(f) ((f)/sfcont(f))

def remove_cont_sfrat(L)
{
	if ( type(L[1]) <= 2 ) {
		T = remove_cont_sfrat([L[0],L[1]*<<0>>]);
		return [T[0],dp_hc(T[1])];
	} else if ( !L[0] )
		return [0,REDCONT(L[1])];
	else if ( !L[1] )
		return [REDCONT(L[0]),0];
	else {
		A0 = REDCONT(L[0]); A1 = REDCONT(L[1]);
		C0 = sdiv(dp_hc(L[0]),dp_hc(A0)); C1 = sdiv(dp_hc(L[1]),dp_hc(A1));
		GCD = sfgcd(C0,C1); M0 = sdiv(C0,GCD); M1 = sdiv(C1,GCD);
		return [M0*A0,M1*A1];
	}
}

def nf_sfrat(B,G,M,PS)
{
	for ( D = 0; G; ) {
		for ( U = 0, L = B; L != []; L = cdr(L) ) {
			if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
				GCD = sfgcd(dp_hc(G),dp_hc(R));
				CG = sdiv(dp_hc(R),GCD); CR = sdiv(dp_hc(G),GCD);
				U = CG*G-dp_subd(G,R)*CR*R;
				if ( !U )
					return [D,M];
				D *= CG; M *= CG;
				break;
			}
		}
		if ( U )
			G = U;
		else {
			D += dp_hm(G); G = dp_rest(G);
		}
	}
	return [D,M];
}

def nf_tab_sfrat(F,TAB)
{
	for ( NM = 0, DN = 1, I = 0; F; F = dp_rest(F) ) {
		T = dp_ht(F);
		for ( ; TAB[I][0] != T; I++);
		NF = TAB[I][1]; N = NF[0]; D = NF[1];
		G = sfgcd(DN,D); DN1 = sdiv(DN,G); D1 = sdiv(D,G);
		NM = D1*NM + DN1*dp_hc(F)*N; DN *= D1;
	}
	return [NM,DN];
}

def adj_dn_sfrat(P,D)
{
	return [(sdiv(D,dp_hc(P[1])))*P[0],dp_ht(P[1])];
}

def tolex_sfrat_main(V,O,NF,GM,Vc,Eval,MB)
{
	if ( MB ) {
		PosDim = 0;
		DIM = length(MB);
		DV = newvect(DIM);
	} else
		PosDim = 1;
	for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
		S = p_terms(car(T),V,2);
		if ( PosDim ) {
			MB = gather_nf_terms(S,NF,V,O);
			DV = newvect(length(MB));
		}
		dp_ord(O); RHS = termstomat_sfrat(NF,map(dp_ptod,cdr(S),V),MB,Vc,Eval);
		dp_ord(O); NHT = nf_tab_gsl_sfrat(dp_ptod(LCM*car(S),V),NF);
		dptov(NHT[0],DV,MB);
		dp_ord(O); B = hen_ttob_gsl_sfrat([DV,NHT[1]],RHS,cdr(S),Vc,Eval);
		if ( !B )
			return 0;
		Len = length(S);
		LCM *= B[1];
		for ( U = LCM*car(S), I = 1; I < Len; I++  )
			U += B[0][I-1]*S[I];
		SL = cons(U,SL);
		if ( dp_gr_print() )
			print(["DN",B[1]]);
	}
	return SL;
}

def termstomat_sfrat(NF,TERMS,MB,Vc,Eval)
{
	DN = NF[1];
	NF = NF[0];
	N = length(MB);
	M = length(TERMS);
	MAT = newmat(N,M);
	W = newvect(N);
	Len = length(NF);
	for ( I = 0; I < M; I++ ) {
		T = TERMS[I];
		for ( K = 0; K < Len; K++ )
			if ( T == NF[K][1] )
				break;
		dptov(NF[K][0],W,MB);
		for ( J = 0; J < N; J++ )
			MAT[J][I] = W[J];
	}
	return [henleq_prep_sfrat(MAT,Vc,Eval),DN];
}

def henleq_prep_sfrat(A,Vc,Eval)
{
	SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
	L = geninv_sf_swap(map(eval_sfrat,A,Vc,Eval)); INV = L[0]; INDEX = L[1];
	AA = newmat(COL,COL);
	for ( I = 0; I < COL; I++ )
		for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL; J++ )
			T[J] = S[J];
	if ( COL != ROW ) {
		RESTA = newmat(ROW-COL,COL);
		for ( ; I < ROW; I++ )
			for ( J = 0, T = RESTA[I-COL], S = A[INDEX[I]]; J < COL; J++ )
				T[J] = S[J];
	} else
		RESTA = 0;
	return [[A,AA,RESTA],L];
}

def hen_ttob_gsl_sfrat(LHS,RHS,TERMS,Vc,Eval)
{
	LDN = LHS[1]; RDN = RHS[1]; LCM = lcm_sfrat(LDN,RDN);
	L1 = sdiv(LCM,LDN); R1 = sdiv(LCM,RDN);
	T0 = time()[0];
	S = henleq_gsl_sfrat(RHS[0],LHS[0]*L1,Vc,Eval);
	if ( dp_gr_print() )
		print(["henleq_gsl_sfrat",time()[0]-T0]);
	if ( !S )
		return 0;
	N = length(TERMS);
	return [S[0],S[1]*R1];
}

def nf_tab_gsl_sfrat(A,NF)
{
	DN = NF[1];
	NF = NF[0];
	TLen = length(NF);
	for ( R = 0; A; A = dp_rest(A) ) {
		HM = dp_hm(A); C = dp_hc(HM); T = dp_ht(HM);
		for ( I = 0; I < TLen; I++ )
			if ( NF[I][1] == T )
				break;
		R += C*NF[I][0];
	}
	return remove_cont_sfrat([R,DN]);
}

def henleq_gsl_sfrat(L,B,Vc,Eval)
{
	Nc = length(Vc);
	if ( Nc > 1 )
		return henleq_gsl_sfrat_higher(L,B,Vc,Eval);

	V0 = Vc[0]; E0 = ptosfp(Eval[0]);

	AL = L[0]; INVL = L[1];
	A = AL[0]; AA = AL[1]; RESTA = AL[2];
	INV = INVL[0]; INDEX = INVL[1];

	SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
	BB = newvect(COL);
	for ( I = 0; I < COL; I++ )
		BB[I] = B[INDEX[I]];
	if ( COL != ROW ) {
		RESTB = newvect(ROW-COL);
		for ( ; I < ROW; I++ )
			RESTB[I-COL] = B[INDEX[I]];	
	} else
		RESTB = 0;

	COUNT = 2;
	if ( !COUNT )
		COUNT = 2;
	X = newvect(size(AA)[0]);
	for ( I = 0, C = BB, CCC = 0, ITOR_FAIL = -1; ; I++ ) {
		if ( zerovector(C) ) {
			X = map(subst,X,V0,V0-E0);
			if ( zerovector(RESTA*X+RESTB) ) {
				if ( dp_gr_print() ) print("end",0);
				return [X,simp_ff(1)];
			} else
				return 0;
		} else if ( COUNT == CCC ) {
			CCC = 0;
			ND = polyvtoratv(X,V0,I);
			if ( ND ) {
				F = map(subst,ND[0],V0,V0-E0);
				LCM = subst(ND[1],V0,V0-E0);
				T = AA*F+LCM*BB;
				if ( zerovector(T) ) {
					T = RESTA*F+LCM*RESTB;
					if ( zerovector(T) ) {
						if ( dp_gr_print() ) print("end",0);
						return [F,LCM];
					} else
						return 0;
				}
			} else {
			}
		} else {
			if ( dp_gr_print() ) print(".",2);
			CCC++;
		}
		XT = -INV*map(subst,C,V0,E0);
		X += XT*V0^I;
		C += AA*XT;
		C = map(sdiv,C,V0-E0);
	}
}

def henleq_gsl_sfrat_higher(L,B,Vc,Eval)
{
	Nc = length(Vc);
	E = map(ptosfp,Eval);

	AL = L[0]; INVL = L[1];
	A = AL[0]; AA0 = AL[1]; RESTA = AL[2];
	INV = INVL[0]; INDEX = INVL[1];

	SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];

	AA = map(mshift,AA0,Vc,E,1);
	BB = newvect(COL);
	for ( I = 0; I < COL; I++ )
		BB[I] = mshift(B[INDEX[I]],Vc,E,1);

	if ( COL != ROW ) {
		RESTB = newvect(ROW-COL);
		for ( ; I < ROW; I++ )
			RESTB[I-COL] = B[INDEX[I]];	
	} else
		RESTB = 0;

	COUNT = 2;
	if ( !COUNT )
		COUNT = 2;
	X = newvect(size(AA)[0]);
	for ( I = 0, R = BB, CCC = 0, ITOR_FAIL = -1; ; I++ ) {
		if ( zerovector(R) ) {
			X = map(mshift,X,Vc,E,-1);
			if ( zerovector(RESTA*X+RESTB) ) {
				if ( dp_gr_print() ) print("end",0);
				return [X,simp_ff(1)];
			} else
				return 0;
		} else if ( COUNT == CCC ) {
			CCC = 0;
			ND = polyvtoratv_higher(X,Vc,I);
			if ( ND ) {
				F = map(mshift,ND[0],Vc,E,-1);
				LCM = mshift(ND[1],Vc,E,-1);
				T = AA*F+LCM*BB;
				if ( zerovector(T) ) {
					T = RESTA*F+LCM*RESTB;
					if ( zerovector(T) ) {
						if ( dp_gr_print() ) print("end",0);
						return [F,LCM];
					} else
						return 0;
				}
			} else {
			}
		} else {
			if ( dp_gr_print() ) print(".",2);
			CCC++;
		}
		RK = map(mtrunc,R,I+1);
		XT = -INV*RK;
		X += XT;
		R += AA*XT;
	}
}

/* V -> V+Sgn*E */

def mshift(F,V,E,Sgn)
{
	N = length(V);
	for ( I = 0; I < N; I++ )
		F = subst(F,V[I],V[I]+Sgn*E[I]);
	return F;
}

/* truncate terms whose degree are higher than or equal to D */

def mtrunc(F,D)
{
	if ( type(F) <= 1 )
		return F;
	else {
		R = 0;
		V = var(F);
		while ( F ) {
			K = deg(F,V);
			C = coef(F,K,V);
			if ( K < D )
				R += mtrunc(C,D-K)*V^K;
			F -= C*V^K;
		}
		return R;
	}
}

/* for 1-dim case */

def polyvtoratv(Vect,V,K)
{
	N = size(Vect)[0];
	R = newvect(N);
	K2 = idiv(K,2);
	Mat = newmat(K2,K2);
	for ( I = 0; I < N; I++ ) {
		T = polytorat(Vect[I],V,Mat,K2);
		if ( T )
			R[I] = T;
		else
			return 0;
	}
	LCM = R[0][1];
	for ( I = 1; I < N; I++ )
		LCM = lcm_sfrat(LCM,R[I][1]);
	for ( I = 0; I < N; I++ )
		R[I] = R[I][0]*sdiv(LCM,R[I][1]);
	return [R,LCM];
}

/* for higher dim case */

def polyvtoratv_higher(Vect,Vc,K)
{
	N = size(Vect)[0];
	R = newvect(N);
	K2 = idiv(K,2);
	for ( I = 0; I < N; I++ ) {
		T = polytorat_higher(Vect[I],Vc,K2);
		if ( T )
			R[I] = T;
		else
			return 0;
	}
	LCM = R[0][1];
	for ( I = 1; I < N; I++ )
		LCM = lcm_sfrat(LCM,R[I][1]);
	for ( I = 0; I < N; I++ )
		R[I] = R[I][0]*sdiv(LCM,R[I][1]);
	return [R,LCM];
}

/* find F = N/D mod V^(2K), deg(N), deg(D) < K */
def polytorat_gcd(F,V,K)
{
	if ( deg(F,V) < K )
		return [F,simp_ff(1)];
	F1 = Mod^(K*2); F2 = F;
	B1 = 0; B2 = 1;
	while ( 1 ) {
		Q = sdiv(F1,F2); 
		F3 = F1-F2*Q;
		B3 = B1-B2*Q;
		if ( deg(F3,V) < K ) {
			if ( deg(B3,V) < K )
				return [F3,B3];
			else
				return 0;
		}
		F1 = F2; F2 = F3;
		B1 = B2; B2 = B3;
	}
}

/* 
 * for 1-dim case
 *
 * solve
 *
 * [fk      ... f1][  a0  ]   
 * [f(k+1)  ... f2][  a1  ] = 0
 * [        ...   ][ ...  ]
 * [f(2k-1) ... fk][a(k-1)]
 */

def polytorat(F,V,Mat,K)
{
	if ( deg(F,V) < K )
		return [F,simp_ff(1)];
	for ( I = 0; I < K; I++ )
		for ( J = 0; J < K; J++ )
			Mat[I][J] = coef(F,I+K-J);
	S = nullspace_ff(Mat);
	MT = S[0]; IND = S[1];
	for ( I = 0; I < K; I++ )
		if ( IND[I] )
			break;
	if ( I == K )
		return 0;
	D = null_to_poly_ff(MT,IND,V)[0];
	N = trunc(F*D,K-1);
	return [N,D];
}

/* find N,D s.t. tdeg(N), tdeg(D) < K, F = N/D mod <V0,...,VM>^(2K) */

def polytorat_higher(F,V,K)
{
	if ( K < 2 ) return 0;
	if ( homogeneous_deg(F) < K )
		return [F,simp_ff(1)];
	D = create_icpoly(V,K);
	C = extract_coef(D*F,V,K,2*K);
	Vc = vars(C);
	G = dp_gr_f_main(C,Vc,0,2);	
	PS = newvect(length(G),map(dp_ptod,G,Vc));
	for ( I = length(G)-1, Ind = [];I >= 0; I-- )
		Ind = cons(I,Ind);
	D = dp_dtop(nf_sfrat(Ind,dp_ptod(D,Vc),1,PS)[0],Vc);
	if ( !D )
		return 0;
	Vp = setminus(vars(D),V);
	D = subst(D,car(Vp),1);
	for ( T = cdr(Vp); T != []; T = cdr(T) )
		D = subst(D,car(T),0);
	N = mtrunc(F*D,K);
	return [N,D];
}

def create_icpoly(V,K)
{
	if ( V == [] )
		return uc();
	R = 0;
	for ( I = 0; I < K; I++ )
		R += create_icpoly(cdr(V),K-I)*V[0]^I;
	return R;
}

/* extract terms whose degrees are in [From,To) */

def extract_coef(F,V,From,To)
{
	if ( V == [] ) {
		if ( From == 0 )
			return [F];
		else 
			return [];
	}
	R = [];
	V0 = V[0];
	for ( I = 0; I < To; I++ ) {
		C = coef(F,I,V0);
		if ( C ) {
			C = extract_coef(C,cdr(V),From>=I?From-I:0,To-I);
			R = append(C,R);
		}
	}
	return R;
}

def saturation_sfrat(F,G,V,Ord)
{
	Vt = ttttt;
	G0 = dp_gr_f_main(cons(Vt*F-1,G),cons(Vt,V),0,[[0,1],[Ord,length(V)]]);
	return elimination(G0,V);
}

def ideal_list_intersection_sfrat(L,V)
{
	R = car(L);
	for ( TL = cdr(L); TL != []; TL = cdr(TL) )
		R = ideal_intersection_sfrat(R,car(TL),V);
	return R;
}

def ideal_intersection_sfrat(A,B,V)
{
	T0 = time();
	Vt = ttttt;
	C = [];
	for ( T = A; T != []; T = cdr(T) )
		C = cons(car(T)*Vt,C);
	for ( T = B; T != []; T = cdr(T) )
		C = cons(car(T)*(1-Vt),C);
	Ord = [[0,1],[0,length(V)]];	
#if 0
	G0 = dp_gr_f_main(C,cons(Vt,V),0,0);
	G = dp_gr_f_main(G0,cons(Vt,V),0,Ord);
#else
	G = dp_gr_f_main(C,cons(Vt,V),0,Ord);
#endif
	T_INT += time()[3]-T0[3];
	return elimination(G,V);
}

/* Shimoyama's gr_fctr */

def idealsqfr_sf(G)
{
	for(LL=[],I=length(G)-1;I>=0;I--) {
		for(A=1,L=sfsqfr(G[I]),J=1;J<length(L);J++)
			A*=L[J][0];
		LL=cons(A,LL);
	}
	return LL;
}

def ideal_uniq(L) /* sub procedure of welldec and normposdec */
{
	for (R = [],I = 0,N=length(L); I < N; I++) {
		if ( R == [] )
			R = append(R,[L[I]]);
		else {
			for (J = 0; J < length(R); J++)
				if ( gb_comp_old(L[I],R[J]) )
					break;
			if ( J == length(R) )
				R = append(R,[L[I]]);
		}
	}
	return R;
}

def ideal_uniq_by_first(L) /* sub procedure of welldec and normposdec */
{
	for (R = [],I = 0,N=length(L); I < N; I++) {
		if ( R == [] )
			R = append(R,[L[I]]);
		else {
			for (J = 0; J < length(R); J++)
				if ( gb_comp_old(L[I][0],R[J][0]) )
					break;
			if ( J == length(R) )
				R = append(R,[L[I]]);
		}
	}
	return R;
}

def prime_irred_sf(TP,VL,Ord)
{
	TP = ideal_uniq(TP);
	N = length(TP);
	for (P = [], I = 0; I < N; I++) {
		for (J = 0; J < N; J++)
			if ( I != J && inclusion_test(TP[J],TP[I],VL,Ord) )
				break;
		if (J == N)
			P = append(P,[TP[I]]);
	}
	return P;
}

def prime_irred_sf_by_first(TP,VL,Ord)
{
	TP = ideal_uniq_by_first(TP);
	N = length(TP);
	for (P = [], I = 0; I < N; I++) {
		for (J = 0; J < N; J++)
			if ( I != J && inclusion_test(car(TP[J]),car(TP[I]),VL,Ord) )
				break;
		if (J == N)
			P = append(P,[TP[I]]);
	}
	return P;
}

def monic_sf_first(L,V)
{
	for ( R = [], T = L; T != []; T = cdr(T) )
		R = cons([monic_sf(car(T)[0],V),car(T)[1]],R);
	return reverse(R);
}

def monic_sf(P,V)
{
	if ( type(P) == 4 )
		return map(monic_sf,P,V);
	else {
		D = dp_ptod(P,V);
		return dp_dtop(D/dp_hc(D),V);
	}
}

def gr_fctr_sf(FL,VL,Ord)
{
	T0 = time();
	COMMONCHECK_SF = 1;
	for (TP = [],I = 0; I<length(FL); I++ ) {
		F = FL[I];
		SF = idealsqfr_sf(F);
		if ( !gb_comp_old(F,SF) ) 
			F = dp_gr_f_main(SF,VL,0,Ord);
		CID_SF=[1];
		SP = gr_fctr_sub_sf(F,VL,Ord);
		TP = append(TP,SP);
	}
	TP = prime_irred_sf(TP,VL,Ord);
	T_GRF += time()[3]-T0[3];
	return TP;
}

def gr_fctr_sub_sf(G,VL,Ord)
{
	if ( length(G) == 1 && type(G[0]) == 1 )
		return [G];
	RL = [];
	for (I = 0; I < length(G); I++) {
		FL = sffctr(G[I]); L = length(FL); N = FL[1][1];
		if (L > 2 || N > 1) {
			TLL = [];
			for (J = 1; J < L; J++) {
				W = cons(FL[J][0],G);
				NG = dp_gr_f_main(W,VL,0,Ord);
				TNG = idealsqfr_sf(NG);
				if ( !gb_comp_old(NG,TNG) )
					NG = dp_gr_f_main(TNG,VL,0,Ord);
				if ( !inclusion_test(CID_SF,NG,VL,Ord) ) {
					DG = gr_fctr_sub_sf(NG,VL,Ord);
					RL = append(RL,DG);
					if ( J <= L-2 
						&& !inclusion_test(CID_SF,NG,VL,Ord)
						&& COMMONCHECK_SF ) {
						CID_SF=ideal_intersection_sfrat(CID_SF,NG,VL);
					}
				}
			}
			break;
		}
	}
	if (I == length(G))
		RL = append([G],RL);
	return RL;
}

def gb_comp_old(A,B)
{
	LA = length(A);
	LB = length(B);
	if ( LA != LB )
		return 0;
	A = newvect(LA,A);
	B = newvect(LB,B);
	A1 = qsort(A);
	B1 = qsort(B);
	for ( I = 0; I < LA; I++ )
		if ( A1[I] != B1[I] && A1[I] != -B1[I] )
			break;
	return I == LA ? 1 : 0;
}
end$