/*
 * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
 * All rights reserved.
 * 
 * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
 * non-exclusive and royalty-free license to use, copy, modify and
 * redistribute, solely for non-commercial and non-profit purposes, the
 * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
 * conditions of this Agreement. For the avoidance of doubt, you acquire
 * only a limited right to use the SOFTWARE hereunder, and FLL or any
 * third party developer retains all rights, including but not limited to
 * copyrights, in and to the SOFTWARE.
 * 
 * (1) FLL does not grant you a license in any way for commercial
 * purposes. You may use the SOFTWARE only for non-commercial and
 * non-profit purposes only, such as academic, research and internal
 * business use.
 * (2) The SOFTWARE is protected by the Copyright Law of Japan and
 * international copyright treaties. If you make copies of the SOFTWARE,
 * with or without modification, as permitted hereunder, you shall affix
 * to all such copies of the SOFTWARE the above copyright notice.
 * (3) An explicit reference to this SOFTWARE and its copyright owner
 * shall be made on your publication or presentation in any form of the
 * results obtained by use of the SOFTWARE.
 * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
 * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
 * for such modification or the source code of the modified part of the
 * SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
 * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
 * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
 * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
 * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
 * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
 * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
 * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
 * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
 * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
 * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
 * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
 * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
 * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
 * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
 * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
 * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
 *
 * $OpenXM: OpenXM_contrib2/asir2018/lib/gr,v 1.1 2018/09/19 05:45:08 noro Exp $
*/

module gr $
  /* Empty for now.  It will be used in a future. */
endmodule $

extern INIT_COUNT,ITOR_FAIL$
extern REMOTE_MATRIX,REMOTE_NF,REMOTE_VARS$

#define MAX(a,b) ((a)>(b)?(a):(b))
#define HigherDim 0
#define ZeroDim   1
#define MiniPoly  2

/* toplevel functions for Groebner basis computation */

def gr(B,V,O)
{
	Ord=dp_ord();
	G = dp_gr_main(B,V,0,1,O);
	dp_ord(Ord);
	return G;
}

def hgr(B,V,O)
{
	Ord=dp_ord();
	G = dp_gr_main(B,V,1,1,O);
	dp_ord(Ord);
	return G;
}

def gr_mod(B,V,O,M)
{
	Ord=dp_ord();
	G = dp_gr_mod_main(B,V,0,M,O);
	dp_ord(Ord);
	return G;
}

def hgr_mod(B,V,O,M)
{
	Ord=dp_ord();
	G = dp_gr_mod_main(B,V,1,M,O);
	dp_ord(Ord);
	return G;
}

/* toplevel functions for change-of-ordering */

def lex_hensel(B,V,O,W,H)
{
	G = dp_gr_main(B,V,H,1,O);
	return tolex(G,V,O,W);
}

def lex_hensel_gsl(B,V,O,W,H)
{
	G = dp_gr_main(B,V,H,1,O);
	return tolex_gsl(G,V,O,W);
}

def gr_minipoly(B,V,O,P,V0,H)
{
	G = dp_gr_main(B,V,H,1,O);
	return minipoly(G,V,O,P,V0);
}

def lex_tl(B,V,O,W,H)
{
	G = dp_gr_main(B,V,H,1,O);
	return tolex_tl(G,V,O,W,H);
}

def tolex_tl(G0,V,O,W,H)
{
	N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
	for ( I = 0; ; I++ ) {
		M = lprime(I);
		if ( !valid_modulus(HM,M) )
			continue;
		if ( ZD ) {
			if ( G3 = dp_gr_main(G0,W,H,-M,3) )
				for ( J = 0; ; J++ )
					if ( G2 = dp_gr_main(G3,W,0,-lprime(J),2) )
						return G2;
		} else if ( G2 = dp_gr_main(G0,W,H,-M,2) )
			return G2;
	}
}

def tolex(G0,V,O,W)
{
	Procs = getopt(procs);

	TM = TE = TNF = 0;
	N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
	if ( ZD )
		MB = dp_mbase(map(dp_ptod,HM,V));
	else
		MB = 0;
	for ( J = 0; ; J++ ) {
		M = lprime(J);
		if ( !valid_modulus(HM,M) )
			continue;
		T0 = time()[0]; 
		if ( ZD ) {
			GM = tolexm(G0,V,O,W,M); 
			dp_ord(2);
			DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
			D = newvect(N); TL = [];
			do 
				TL = cons(dp_dtop(dp_vtoe(D),W),TL);
			while ( nextm(D,DL,N) );
		} else {
			HVN = "h";
			WN = map(rtostr,W);
			while ( member(HVN,WN) ) HVN += "h";
			HV = strtov(HVN);
			G0H = map(homogenize,G0,W,HV);
			GMH = nd_gr(G0H,append(W,[HV]),M,1);
			GMH=map(subst,GMH,HV,1);
			GM = nd_gr_postproc(GMH,W,M,2,0);
			dp_ord(2);
			for ( T = GM, S = 0; T != []; T = cdr(T) )
				for ( D = dp_ptod(car(T),V); D; D = dp_rest(D) )
					S += dp_ht(D);
			TL = dp_terms(S,V);
		}
		TM += time()[0] - T0;
		T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],ZD)[0];
		TNF += time()[0] - T0;
		T0 = time()[0];
		if ( type(Procs) != -1 )
			R = tolex_d_main(V,O,NF,GM,M,MB,Procs);
		else
			R = tolex_main(V,O,NF,GM,M,MB);
		TE += time()[0] - T0;
		if ( R ) {
			if ( dp_gr_print() )
				print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
			return R;
		}
	}
}

def tolex_gsl(G0,V,O,W)
{
	TM = TE = TNF = 0;
	N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
	MB = dp_mbase(map(dp_ptod,HM,V));
	if ( !ZD )
		error("tolex_gsl : ideal is not zero-dimensional!");
	for ( J = 0; ; J++ ) {
		M = lprime(J);
		if ( !valid_modulus(HM,M) )
			continue;
		T0 = time()[0]; GM = tolexm(G0,V,O,W,M); TM += time()[0] - T0;
		dp_ord(2);
		DL = map(dp_etov,map(dp_ht,map(dp_ptod,GM,W)));
		D = newvect(N); TL = [];
		do 
			TL = cons(dp_dtop(dp_vtoe(D),W),TL);
		while ( nextm(D,DL,N) );
		L = npos_check(DL); NPOSV = L[0]; DIM = L[1];
		if ( NPOSV >= 0 ) {
			if ( dp_gr_print() ) print("shape base");
			V0 = W[NPOSV];
			T0 = time()[0]; NFL = gennf(G0,TL,V,O,V0,1);
			TNF += time()[0] - T0;
			T0 = time()[0];
			R = tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB);
			TE += time()[0] - T0;
		} else {
			T0 = time()[0]; NF = gennf(G0,TL,V,O,W[N-1],1)[0];
			TNF += time()[0] - T0;
			T0 = time()[0];
			R = tolex_main(V,O,NF,GM,M,MB);
			TE += time()[0] - T0;
		}
		if ( R ) {
			if ( dp_gr_print() )
				print("mod="+rtostr(TM)+",nf="+rtostr(TNF)+",eq="+rtostr(TE));
			return R;
		}
	}
}

def termstomat(NF,TERMS,MB,MOD)
{
	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(MAT,MOD),DN];
}

def tolex_gsl_main(G0,V,O,W,NFL,NPOSV,GM,M,MB)
{
	NF = NFL[0]; PS = NFL[1]; GI = NFL[2];
	V0 = W[NPOSV]; N = length(W);
	DIM = length(MB);
	DV = newvect(DIM);
	TERMS = gather_terms(GM,W,M,NPOSV);
	Len = length(TERMS);
	dp_ord(O); RHS = termstomat(NF,map(dp_ptod,TERMS,V),MB,M);
	for ( T = GM; T != []; T = cdr(T) )
		if ( vars(car(T)) == [V0]	)
			break;
	dp_ord(0); NHT = nf_tab_gsl(dp_ptod(V0^deg(car(T),V0),V),NF);
	dptov(NHT[0],DV,MB);
	B = hen_ttob_gsl([DV,NHT[1]],RHS,TERMS,M);
	if ( !B )
		return 0;
	for ( I = 0, U = B[1]*V0^deg(car(T),V0); I < Len; I++ )
		U += B[0][I]*TERMS[I];
	DN0 = diff(U,V0);
	dp_ord(O); DN0NF = nf_tab_gsl(dp_ptod(DN0,V),NF);
	SL = [[V0,U,DN0]];
	for ( I = N-1, LCM = 1; I >= 0; I-- ) {
		if ( I == NPOSV )
			continue;
		V1 = W[I];
		dp_ord(O); L = nf(GI,DN0NF[0]*dp_ptod(-LCM*V1,V),DN0NF[1],PS);
		L = remove_cont(L);
		dptov(L[0],DV,MB);
		dp_ord(O); B = hen_ttob_gsl([DV,L[1]],RHS,TERMS,M);
		if ( !B )
			return 0;
		for ( K = 0, R = 0; K < Len; K++ )
			R += B[0][K]*TERMS[K];
		LCM *= B[1];
		SL = cons(cons(V1,[R,LCM]),SL);
		if ( dp_gr_print() )
			print(["DN",B[1]]);
	}
	return SL;
}

def hen_ttob_gsl(LHS,RHS,TERMS,M)
{
	LDN = LHS[1]; RDN = RHS[1]; LCM = ilcm(LDN,RDN);
	L1 = idiv(LCM,LDN); R1 = idiv(LCM,RDN);
	T0 = time()[0];
	S = henleq_gsl(RHS[0],LHS[0]*L1,M);
	if ( dp_gr_print() )
		print(["henleq_gsl",time()[0]-T0]);
	N = length(TERMS);
	return [S[0],S[1]*R1];
}

def	gather_terms(GM,W,M,NPOSV)
{
	N = length(W); V0 = W[NPOSV];
	for ( T = GM; T != []; T = cdr(T) ) {
		if ( vars(car(T)) == [V0] )
			break;
	}
	U = car(T); DU = diff(U,V0);
	R = tpoly(cdr(p_terms(U,W,2)));
	for ( I = 0; I < N; I++ ) {
		if ( I == NPOSV )
			continue;
		V1 = W[I];
		for ( T = GM; T != []; T = cdr(T) )
			if ( member(V1,vars(car(T))) )
				break;
		P = car(T);
		R += tpoly(p_terms(srem(DU*coef(P,0,V1),U,M),W,2));
	}
	return p_terms(R,W,2);
}

def tpoly(L)
{
	for ( R = 0; L != []; L = cdr(L) )
		R += car(L);
	return R;
}

def dptov(P,W,MB)
{
	N = size(W)[0];
	for ( I = 0; I < N; I++ )
		W[I] = 0;
	for ( I = 0, S = MB; P; P = dp_rest(P) ) {
		HM = dp_hm(P); C = dp_hc(HM); T = dp_ht(HM);
		for ( ; T != car(S); S = cdr(S), I++ );
		W[I] = C;
		I++; S = cdr(S);
	}
}

def tolex_main(V,O,NF,GM,M,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(NF,map(dp_ptod,cdr(S),V),MB,M);
		dp_ord(O); NHT = nf_tab_gsl(dp_ptod(LCM*car(S),V),NF);
		dptov(NHT[0],DV,MB);
		dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
		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];
		R = ptozp(U);
		SL = cons(R,SL);
		if ( dp_gr_print() )
			print(["DN",B[1]]);
	}
	return SL;
}

def tolex_d_main(V,O,NF,GM,M,MB,Procs)
{
	map(ox_reset,Procs);
	/* register data in servers */
	map(ox_cmo_rpc,Procs,"register_data_for_find_base",NF,V,O,MB,M);
	/* discard return value in stack */
	map(ox_pop_cmo,Procs);
	Free = Procs;
	Busy = [];
	T = GM;
	SL = [];
	while ( T != [] || Busy != []  ){
		if ( Free == [] || T == [] ) {
			/* someone is working; wait for data */
			Ready = ox_select(Busy);
			Busy = setminus(Busy,Ready);
			Free = append(Ready,Free);
			for ( ; Ready != []; Ready = cdr(Ready) )
				SL = cons(ox_get(car(Ready)),SL);
		} else {
			P = car(Free);
			Free = cdr(Free);
			Busy = cons(P,Busy);
			Template = car(T);
			T = cdr(T);
			ox_cmo_rpc(P,"find_base",Template);
			ox_push_cmd(P,262); /* 262 = OX_popCMO */
		}
	}
	return SL;
}

struct find_base_data { NF,V,O,MB,M,PosDim,DV }$
extern Find_base$

def register_data_for_find_base(NF,V,O,MB,M)
{
	Find_base = newstruct(find_base_data);
	Find_base->NF = NF;
	Find_base->V = V;
	Find_base->O = O;
	Find_base->M = M;
	Find_base->MB = MB;

	if ( MB ) {
		Find_base->PosDim = 0;
		DIM = length(MB);
		Find_base->DV = newvect(DIM);
	} else
		Find_base->PosDim = 1;
}

def find_base(S) {
	NF = Find_base->NF;
	V = Find_base->V;
	O = Find_base->O;
	MB = Find_base->MB;
	M = Find_base->M;
	PosDim = Find_base->PosDim;
	DV = Find_base->DV;

	S = p_terms(S,V,2);
	if ( PosDim ) {
		MB = gather_nf_terms(S,NF,V,O);
		DV = newvect(length(MB));
	}
	dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
	dp_ord(O); NHT = nf_tab_gsl(dp_ptod(car(S),V),NF);
	dptov(NHT[0],DV,MB);
	dp_ord(O); B = hen_ttob_gsl([DV,NHT[1]],RHS,cdr(S),M);
	if ( !B )
		return 0;
	Len = length(S);
	for ( U = B[1]*car(S), I = 1; I < Len; I++  )
		U += B[0][I-1]*S[I];
	R = ptozp(U);
	return R;
}

/*
 * NF = [Pairs,DN]
 *  Pairs = [[NF1,T1],[NF2,T2],...]
 */

def gather_nf_terms(S,NF,V,O)
{
	R = 0;
	for ( T = S; T != []; T = cdr(T) ) {
		DT = dp_ptod(car(T),V);
		for ( U = NF[0]; U != []; U = cdr(U) )
			if ( car(U)[1] == DT ) {
				R += tpoly(dp_terms(car(U)[0],V));
				break;
			}
	}
	return map(dp_ptod,p_terms(R,V,O),V);
}

def reduce_dn(L)
{
	NM = L[0]; DN = L[1]; V = vars(NM);
	T = remove_cont([dp_ptod(NM,V),DN]);
	return [dp_dtop(T[0],V),T[1]];
}

/* a function for computation of  minimal polynomial */

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

	Pin = P;
	P = ptozp(P);
	CP = sdiv(P,Pin);
	G1 = cons(V0-P,G0);
	O1 = [[0,1],[O,length(V)]];
	V1 = cons(V0,V);
	W = append(V,[V0]);	

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

	for ( J = 0; ; J++ ) {
		M = lprime(J);
		if ( !valid_modulus(HM,M) )
			continue;
		MP = minipolym(G0,V,O,P,V0,M);
		for ( D = deg(MP,V0), TL = [], J = 0; J <= D; J++ )
			TL = cons(V0^J,TL);
		NF = gennf(G1,TL,V1,O1,V0,1)[0];
		R = tolex_main(V1,O1,NF,[MP],M,MB);
		return ptozp(subst(R[0],V0,CP*V0));
	}
}

/* subroutines */

def gennf(G,TL,V,O,V0,FLAG)
{
	F = dp_gr_flags();
	for ( T = F; T != []; T = cdr(T) ) {
		Key = car(T); T = cdr(T);
		if ( Key == "Demand" ) {
			Dir = car(T); break;
		}
	}
	if ( Dir )
		return gennf_demand(G,TL,V,O,V0,FLAG,Dir);
	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(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(dp_true_nf(GI,U*MB[I],PS,1))];
			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(L[0],UTAB);
				NF = [NF[0],dp_hc(L[1])*NF[1]*T];
			} else
				NF = nf(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
		} else
			NF = nf(GI,T,T,PS);
		NF = remove_cont(NF);
		H = cons(NF,H);
		LCM = ilcm(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,H,LCM),LCM],PS,GI];
}

def gennf_demand(G,TL,V,O,V0,FLAG,Dir)
{
	N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len);
	NTL = length(TL);
	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);

	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(dp_true_nf(GI,U*MB[I],PS,1))];
			if ( dp_gr_print() )
				print(".",2);
		}
		if ( dp_gr_print() )
			print("");
		TTAB = time()[0]-T0;
	}

	T0 = time()[0];
	for ( LCM = 1, Index = 0, H = []; DTL != []; Index++ ) {
		if ( dp_gr_print() )
			print(".",2);
		T = car(DTL); DTL = cdr(DTL);
		if ( L = search_redble(T,H) ) {
			L = nf_load(Dir,L[0]);
			DD = dp_subd(T,L[1]);
			if ( USE_TAB && (DD == U) ) {
				NF = nf_tab(L[0],UTAB);
				NF = [NF[0],dp_hc(L[1])*NF[1]*T];
			} else
				NF = nf(GI,L[0]*dp_subd(T,L[1]),dp_hc(L[1])*T,PS);
		} else
			NF = nf(GI,T,T,PS);
		NF = remove_cont(NF);
		nf_save(NF,Dir,Index);
		H = cons([Index,NF[1]],H);
		LCM = ilcm(LCM,dp_hc(NF[1]));
	}
	TNF = time()[0]-T0;
	if ( dp_gr_print() )
		print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")");
	
	for ( I = 0; I < NTL; I++ ) {
		NF = nf_load(Dir,I);
		NF = adj_dn(NF,LCM);
		nf_save(NF,Dir,I);
	}
	for ( H = [], I = NTL-1; I >= 0; I-- )
		H = cons(nf_load(Dir,I),H);
	return [[H,LCM],PS,GI];
}

def nf_load(Dir,I)
{
	return bload(Dir+"/nf"+rtostr(I));
}

def nf_save(NF,Dir,I)
{
	bsave(NF,Dir+"/nf"+rtostr(I));
}

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

def hen_ttob(T,NF,LHS,V,MOD)
{
	if ( length(T) == 1 )
		return car(T);
	T0 = time()[0]; M = etom(leq_nf(T,NF,LHS,V)); TE = time()[0] - T0;
	T0 = time()[0]; U = henleq(M,MOD); TH = time()[0] - T0;
	if ( dp_gr_print() ) {
		print("(etom="+rtostr(TE)+" hen="+rtostr(TH)+")");
	}
	return U ? vtop(T,U,LHS) : 0;
}

def vtop(S,L,GSL)
{
	U = L[0]; H = L[1];
	if ( GSL ) {
		for ( A = 0, I = 0; S != []; S = cdr(S), I++ )
			A += U[I]*car(S);
		return [A,H];
	} else {
		for ( A = H*car(S), S = cdr(S), I = 0; S != []; S = cdr(S), I++ )
			A += U[I]*car(S);
		return ptozp(A);
	}
}

/* broken */

def leq_nf(TL,NF,LHS,V)
{
	TLen = length(NF);
	T = newvect(TLen); M = newvect(TLen);
	for ( I = 0; I < TLen; I++ ) {
		T[I] = dp_ht(NF[I][1]);
		M[I] = dp_hc(NF[I][1]);
	}
	Len = length(TL); INDEX = newvect(Len); COEF = newvect(Len);
	for ( L = TL, J = 0; L != []; L = cdr(L), J++ ) {
		D = dp_ptod(car(L),V);
		for ( I = 0; I < TLen; I++ )
			if ( D == T[I] )
				break;
		INDEX[J] = I; COEF[J] = strtov("u"+rtostr(J));
	}
	if ( !LHS ) {
		COEF[0] = 1; NM = 0; DN = 1;
	} else {
		NM = LHS[0]; DN = LHS[1];
	}
	for ( J = 0, S = -NM; J < Len; J++ ) {
		DNJ = M[INDEX[J]];
		GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD;
		S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J];
		DN *= CS;
	}
	for ( D = S, E = []; D; D = dp_rest(D) )
		E = cons(dp_hc(D),E);
	BOUND = LHS ? 0 : 1;
	for ( I = Len - 1, W = []; I >= BOUND; I-- )	
			W = cons(COEF[I],W);
	return [E,W];	
}

def nf_tab(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 = igcd(DN,D); DN1 = idiv(DN,G); D1 = idiv(D,G);
		NM = D1*NM + DN1*dp_hc(F)*N; DN *= D1;
	}
	return [NM,DN];
}

def nf_tab_gsl(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([R,DN]);
}

def redble(D1,D2,N)
{
	for ( I = 0; I < N; I++ )
		if ( D1[I] > D2[I] )
			break;
	return I == N ? 1 : 0;
}

def tolexm(G,V,O,W,M)
{
	N = length(V); Len = length(G);
	dp_ord(O); setmod(M); PS = newvect(Len);
	for ( I = 0, T = G; T != []; T = cdr(T), I++ )
		PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
	for ( I = Len-1, HL = []; I >= 0; I-- )
		HL = cons(dp_ht(PS[I]),HL);
	G2 = tolexm_main(PS,HL,V,W,M,ZeroDim);
	L = map(dp_dtop,G2,V);
	return L;
}

def tolexm_main(PS,HL,V,W,M,FLAG)
{
	N = length(W); D = newvect(N); Len = size(PS)[0];
	for ( I = Len - 1, GI = []; I >= 0; I-- )
		GI = cons(I,GI);
	MB = dp_mbase(HL); DIM = length(MB);
	U = dp_mod(dp_ptod(W[N-1],V),M,[]);
	UTAB = newvect(DIM);
	for ( I = 0; I < DIM; I++ ) {
		if ( dp_gr_print() )
			print(".",2);
		UTAB[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
	}
	if ( dp_gr_print() )
		print("");
	T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
	H = G = [[T,T]];
	DL = []; G2 = [];
	TNF = 0;
	while ( 1 ) {
		if ( dp_gr_print() )
			print(".",2);
		S = nextm(D,DL,N);
		if ( !S )
			break;
		T = dp_mod(dp_ptod(dp_dtop(dp_vtoe(D),W),V),M,[]);
		T0 = time()[0];
		if ( L = search_redble(T,H) ) {
			DD = dp_mod(dp_subd(T,L[1]),M,[]);
			if ( DD == U )
				NT = dp_nf_tab_mod(L[0],UTAB,M);
			else
				NT = dp_nf_mod(GI,L[0]*DD,PS,1,M);
		} else
			NT = dp_nf_mod(GI,T,PS,1,M);
		TNF += time()[0] - T0;
		H = cons([NT,T],H);
		T0 = time()[0];
		L = dp_lnf_mod([NT,T],G,M); N1 = L[0]; N2 = L[1];
		TLNF += time()[0] - T0;
		if ( !N1 ) {
			G2 = cons(N2,G2);
			if ( FLAG == MiniPoly )
				break;
			D1 = newvect(N);
			for ( I = 0; I < N; I++ )
				D1[I] = D[I];
			DL = cons(D1,DL);
		} else
			G = insert(G,L);
	}
	if ( dp_gr_print() )
		print("tolexm(nfm="+rtostr(TNF)+" lnfm="+rtostr(TLNF)+")");
	return G2;
}

def minipolym(G,V,O,P,V0,M)
{
	N = length(V); Len = length(G);
	dp_ord(O); setmod(M); PS = newvect(Len);
	for ( I = 0, T = G; T != []; T = cdr(T), I++ )
		PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
	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_mod(dp_ptod(P,V),M,[]);
	for ( I = 0; I < DIM; I++ )
		UT[I] = [MB[I],dp_nf_mod(GI,U*dp_mod(MB[I],M,[]),PS,1,M)];
	T = dp_mod(<<0>>,M,[]); TT = dp_mod(dp_ptod(1,V),M,[]);
	G = H = [[TT,T]]; TNF = TLNF = 0;
	for ( I = 1; ; I++ ) {
		T = dp_mod(<<I>>,M,[]);
		T0 = time()[0]; NT = dp_nf_tab_mod(H[0][0],UT,M); TNF += time()[0] - T0;
		H = cons([NT,T],H);
		T0 = time()[0]; L = dp_lnf_mod([NT,T],G,M); TLNF += time()[0] - T0;
		if ( !L[0] ) {
			if ( dp_gr_print() ) print(["nfm",TNF,"lnfm",TLNF]);
			return dp_dtop(L[1],[V0]);
		} else
			G = insert(G,L);
	}
}

def nextm(D,DL,N)
{
	for ( I = N-1; I >= 0; ) {		
		D[I]++;
		for ( T = DL; T != []; T = cdr(T) )
			if ( car(T) == D )
				return 1;
			else if ( redble(car(T),D,N) )
				break;
		if ( T != [] ) {
			for ( J = N-1; J >= I; J-- )
				D[J] = 0;
			I--;
		} else
			break;
	}
	if ( I < 0 )
		return 0;
	else
		return 1;
}

def search_redble(T,G)
{
	for ( ; G != []; G = cdr(G) )
		if ( dp_redble(T,car(G)[1]) )
			return car(G);
	return 0;
}

def insert(G,A)
{
	if ( G == [] )
		return [A];
	else if ( dp_ht(car(A)) > dp_ht(car(car(G))) )
		return cons(A,G);
	else
		return cons(car(G),insert(cdr(G),A));
}

#if 0
def etom(L) {
	E = L[0]; W = L[1];
	LE = length(E); LW = length(W);
	M = newmat(LE,LW+1);
	for(J=0;J<LE;J++) {
		for ( T = E[J]; T && (type(T) == 2); )
			for ( V = var(T), I = 0; I < LW; I++ )
				if ( V == W[I] ) {
					M[J][I] = coef(T,1,V);
					T = coef(T,0,V);
				}
		M[J][LW] = T;
	}
	return M;
}
#endif

def etom(L) {
	E = L[0]; W = L[1];
	LE = length(E); LW = length(W);
	M = newmat(LE,LW+1);
	for(J=0;J<LE;J++) {
		for ( I = 0, T = E[J]; I < LW; I++ ) {
			M[J][I] = coef(T,1,W[I]); T = coef(T,0,W[I]);
		}
		M[J][LW] = T;
	}
	return M;
}

def calcb_old(M) {
	N = 2*M;
	T = gr_sqrt(N);
	if ( T^2 <= N && N < (T+1)^2 )
		return idiv(T,2);
	else
		error("afo");
}

def calcb_special(PK,P,K) { /* PK = P^K */
	N = 2*PK;
	T = sqrt_special(N,2,P,K);
	if ( T^2 <= N && N < (T+1)^2 )
		return idiv(T,2);
	else
	error("afo");
}

def sqrt_special(A,C,M,K) { /* A = C*M^K */
	L = idiv(K,2); B = M^L;
	if ( K % 2 )
		C *= M;
	D = 2^K; X = idiv((gr_sqrt(C*D^2)+1)*B,D)+1;
	while ( 1 )
		if ( (Y = X^2) <= A )
			return X;
		else
			X = idiv(A + Y,2*X);
}

def gr_sqrt(A) {
	for ( J = 0, T = A; T >= 2^27; J++ ) {
		T = idiv(T,2^27)+1;
	}
	for ( I = 0; T >= 2; I++ ) {
		S = idiv(T,2);
		if ( T = S+S )
			T = S;
		else
			T = S+1;
	}
	X = (2^27)^idiv(J,2)*2^idiv(I,2);
	while ( 1 ) {
		if ( (Y=X^2) < A )
			X += X;
		else if ( Y == A )
			return X;
		else
			break;
	}
	while ( 1 )
		if ( (Y = X^2) <= A )
			return X;
		else
			X = idiv(A + Y,2*X);
}

#define ABS(a) ((a)>=0?(a):(-a))

def inttorat_asir(C,M,B)
{
	if ( M < 0 )
		M = -M;
	C %= M;
	if ( C < 0 )
		C += M;
	U1 = 0; U2 = M; V1 = 1; V2 = C;
	while ( V2 >= B ) {
		L = iqr(U2,V2); Q = L[0]; R2 = L[1];
		R1 = U1 - Q*V1;
		U1 = V1; U2 = V2;
		V1 = R1; V2 = R2;
	}
	if ( ABS(V1) >= B )
		return 0;
	else
	if ( V1 < 0 )
		return [-V2,-V1];
	else
		return [V2,V1];
}

def intvtoratv(V,M,B) {
	if ( !B )
		B = 1;
	N = size(V)[0];
	W = newvect(N);
	if ( ITOR_FAIL >= 0 ) {
		if ( V[ITOR_FAIL] ) {
			T = inttorat(V[ITOR_FAIL],M,B);
			if ( !T ) {
				if ( dp_gr_print() ) {
					print("F",2);
				}
				return 0;
			}
		}
	}
	for ( I = 0, DN = 1; I < N; I++ )
		if ( V[I] ) {
			T = inttorat((V[I]*DN) % M,M,B);
			if ( !T ) {
				ITOR_FAIL = I;
				if ( dp_gr_print() ) {
#if 0
					print("intvtoratv : failed at I = ",0); print(ITOR_FAIL);
#endif
					print("F",2);
				}
				return 0;
			} else {
				for( J = 0; J < I; J++ )
					W[J] *= T[1];
				W[I] = T[0]; DN *= T[1];
			}
		}
	return [W,DN];
}

def nf(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 = igcd(dp_hc(G),dp_hc(R));
				CG = idiv(dp_hc(R),GCD); CR = idiv(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 remove_cont(L)
{
	if ( type(L[1]) == 1 ) {
		T = remove_cont([L[0],L[1]*<<0>>]);
		return [T[0],dp_hc(T[1])];
	} else if ( !L[0] )
		return [0,dp_ptozp(L[1])];
	else if ( !L[1] )
		return [dp_ptozp(L[0]),0];
	else {
		A0 = dp_ptozp(L[0]); A1 = dp_ptozp(L[1]);
		C0 = idiv(dp_hc(L[0]),dp_hc(A0)); C1 = idiv(dp_hc(L[1]),dp_hc(A1));
		GCD = igcd(C0,C1); M0 = idiv(C0,GCD); M1 = idiv(C1,GCD);
		return [M0*A0,M1*A1];
	}
}

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

def union1(A,E)
{
	if ( A == [] )
		return [E];
	else if ( car(A) == E )
		return A;
	else
		return cons(car(A),union1(cdr(A),E));
}

def setminus(A,B) {
	for ( T = reverse(A), R = []; T != []; T = cdr(T) ) {
		for ( S = B, M = car(T); S != []; S = cdr(S) )
			if ( car(S) == M )
				break;
		if ( S == [] )
			R = cons(M,R);
	}
	return R;
}

def member(A,L) {
	for ( ; L != []; L = cdr(L) )
		if ( A == car(L) )
			return 1;
	return 0;
}

/* several functions for computation of normal forms etc. */

def p_nf(P,B,V,O) {
	dp_ord(O); DP = dp_ptod(P,V);
	N = length(B); DB = newvect(N);
	for ( I = N-1, IL = []; I >= 0; I-- ) {
		DB[I] = dp_ptod(B[I],V);
		if ( DB[I] ) IL = cons(I,IL);
	}
	return dp_dtop(dp_nf(IL,DP,DB,1),V);
}

def p_true_nf(P,B,V,O) {
	dp_ord(O); DP = dp_ptod(P,V);
	N = length(B); DB = newvect(N);
	for ( I = N-1, IL = []; I >= 0; I-- ) {
		DB[I] = dp_ptod(B[I],V);
		if ( DB[I] ) IL = cons(I,IL);
	}
	L = dp_true_nf(IL,DP,DB,1);
	return [dp_dtop(L[0],V),L[1]];
}

def p_nf_mod(P,B,V,O,Mod) {
	setmod(Mod);
	dp_ord(O); DP = dp_mod(dp_ptod(P,V),Mod,[]);
	N = length(B); DB = newvect(N);
	for ( I = N-1, IL = []; I >= 0; I-- ) {
		DB[I] = dp_mod(dp_ptod(B[I],V),Mod,[]);
		IL = cons(I,IL);
	}
	return dp_dtop(dp_nf_mod(IL,DP,DB,1,Mod),V);
}

def p_terms(D,V,O)
{
	dp_ord(O);
	for ( L = [], T = dp_ptod(D,V); T; T = dp_rest(T) )
		L = cons(dp_dtop(dp_ht(T),V),L);
	return reverse(L);
}

def dp_terms(D,V)
{
	for ( L = [], T = D; T; T = dp_rest(T) )
		L = cons(dp_dtop(dp_ht(T),V),L);
	return reverse(L);
}

def gb_comp(A,B)
{
	LA = length(A);
	LB = length(B);
	if ( LA != LB )
		return 0;
	A = newvect(LA,A);
	B = newvect(LB,B);
	for ( I = 0; I < LA; I++ )
		A[I] *= headsgn(A[I]);
	for ( I = 0; I < LB; I++ )
		B[I] *= headsgn(B[I]);
	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;
}

def zero_dim(G,V,O) {
	dp_ord(O);
	HL = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
	for ( L = []; HL != []; HL = cdr(HL) )
		if ( length(vars(car(HL))) == 1 )
			L = cons(car(HL),L);
	return length(vars(L)) == length(V) ? 1 : 0;
}

def hmlist(G,V,O) {
	dp_ord(O);
	return map(dp_dtop,map(dp_hm,map(dp_ptod,G,V)),V);
}

def valid_modulus(HL,M) {
	V = vars(HL);
	for ( T = HL; T != []; T = cdr(T) )
		if ( !dp_mod(dp_ptod(car(T),V),M,[]) )
			break;
	return T == [] ? 1 : 0;
}

def npos_check(DL) {
	N = size(car(DL))[0];
	if ( length(DL) != N )
		return [-1,0];
	D = newvect(N);
	for ( I = 0; I < N; I++ ) {
		for ( J = 0; J < N; J++ )
			D[J] = 0;
		D[I] = 1;
		for ( T = DL; T != []; T = cdr(T) )
			if ( D == car(T) )
				break;
		if ( T != [] )
			DL = setminus(DL,[car(T)]);
	}
	if ( length(DL) != 1 )
		return [-1,0];
	U = car(DL);
	for ( I = 0, J = 0, I0 = -1; I < N; I++ )
		if ( U[I] ) {
			I0 = I; J++;
		}
	if ( J != 1 )
		return [-1,0];
	else
		return [I0,U[I0]];
}

def mult_mat(L,TAB,MB)
{
	A = L[0]; DN0 = L[1];
	for ( NM = 0, DN = 1, I = 0; A; A = dp_rest(A) ) {
		H = dp_ht(A);
		for ( ; MB[I] != H; I++ );
		NM1 = TAB[I][0]; DN1 = TAB[I][1]; I++;
		GCD = igcd(DN,DN1); C = DN1/GCD; C1 = DN/GCD;
		NM = C*NM + C1*dp_hc(A)*NM1;
		DN *= C;
	}
	Z=remove_cont([NM,DN*DN0]);
	return Z;
}

def sepm(MAT)
{
	S = size(MAT); N = S[0]; M = S[1]-1;
	A = newmat(N,M); B = newvect(N);
	for ( I = 0; I < N; I++ )
		for ( J = 0, T1 = MAT[I], T2 = A[I]; J < M; J++ )
			T2[J] = T1[J];
	for ( I = 0; I < N; I++ )
		B[I] = MAT[I][M];
	return [A,B];
}

def henleq(M,MOD)
{
	SIZE = size(M); ROW = SIZE[0]; COL = SIZE[1];
	W = newvect(COL);
	L = sepm(M); A = L[0]; B = L[1];
	COUNT = INIT_COUNT?INIT_COUNT:idiv(max_mag(M),54);
	if ( !COUNT )
		COUNT = 1;

	TINV = TC = TR = TS = TM = TDIV = 0;

	T0 = time()[0];
	L = geninvm_swap(A,MOD); INV = L[0]; INDEX = L[1];
	TS += time()[0] - T0;
	
	COL1 = COL - 1;
	AA = newmat(COL1,COL1); BB = newvect(COL1);
	for ( I = 0; I < COL1; I++ ) {
		for ( J = 0, T = AA[I], S = A[INDEX[I]]; J < COL1; J++ )
			T[J] = S[J];
		BB[I] = B[INDEX[I]];
	}
	if ( COL1 != ROW ) {
		RESTA = newmat(ROW-COL1,COL1); RESTB = newvect(ROW-COL1);
		for ( ; I < ROW; I++ ) {
			for ( J = 0, T = RESTA[I-COL1], S = A[INDEX[I]]; J < COL1; J++ )
				T[J] = S[J];
			RESTB[I-COL1] = B[INDEX[I]];	
		}
	} else
		RESTA = RESTB = 0;

	MOD2 = idiv(MOD,2);
	for ( I = 0, C = BB, X = 0, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
		I++, PK *= MOD ) {
		if ( COUNT == CCC ) {
			CCC = 0;
			T0 = time()[0]; 
			ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
			TR += time()[0]-T0;
			if ( ND ) {
				T0 = time()[0];
				F = ND[0]; LCM = ND[1]; T = AA*F+LCM*BB;
				TM += time()[0]-T0;
				if ( zerovector(T) ) {
					T0 = time()[0]; T = RESTA*F+LCM*RESTB; TM += time()[0]-T0;
					if ( zerovector(T) ) {
#if 0
						if ( dp_gr_print() ) print(["init",TS,"pinv",TINV,"c",TC,"div",TDIV,"rat",TR,"mul",TM]);
#endif
						if ( dp_gr_print() ) print("end",2);
						return [F,LCM];
					} else
						return 0;
				}
			} else {
#if 0
				if ( dp_gr_print() ) print(I);
#endif
			}
		} else {
#if 0
			if ( dp_gr_print() ) print([I,TINV,TC,TDIV]);
#endif
			if ( dp_gr_print() ) print(".",2);
			CCC++;
		}
		T0 = time()[0];
		XT = sremainder(INV*sremainder(-C,MOD),MOD);
		XT = map(adj_sgn,XT,MOD,MOD2);
		TINV += time()[0] - T0;
		X += XT*PK;
		T0 = time()[0];
		C += mul_mat_vect_int(AA,XT);
		TC += time()[0] - T0;
		T0 = time()[0]; C = map(idiv,C,MOD); TDIV += time()[0] - T0;
	}
}

def henleq_prep(A,MOD)
{
	SIZE = size(A); ROW = SIZE[0]; COL = SIZE[1];
	L = geninvm_swap(A,MOD); 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 henleq_gsl(L,B,MOD)
{
	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 = INIT_COUNT?INIT_COUNT:idiv(MAX(max_mag(A),max_mag_vect(B)),54);
	if ( !COUNT )
		COUNT = 1;
	MOD2 = idiv(MOD,2);
	X = newvect(size(AA)[0]);
	for ( I = 0, C = BB, PK = 1, CCC = 0, ITOR_FAIL = -1; ;
		I++, PK *= MOD ) {
		if ( zerovector(C) )
			if ( zerovector(RESTA*X+RESTB) ) {
				if ( dp_gr_print() ) print("end",0);
				return [X,1];
			} else
				return 0;
		else if ( COUNT == CCC ) {
			CCC = 0;
			ND = intvtoratv(X,PK,ishift(calcb_special(PK,MOD,I),32));
			if ( ND ) {
				F = ND[0]; LCM = ND[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++;
		}
		XT = sremainder(INV*sremainder(-C,MOD),MOD);
		XT = map(adj_sgn,XT,MOD,MOD2);
		X += XT*PK;
		C += mul_mat_vect_int(AA,XT);
		C = map(idiv,C,MOD);
	}
}

def adj_sgn(A,M,M2)
{
	return A > M2 ? A-M : A;
}

def zerovector(C)
{
	if ( !C )
		return 1;
	for ( I = size(C)[0]-1; I >= 0 && !C[I]; I-- );
	if ( I < 0 )
		return 1;
	else
		return 0;
}

def solvem(INV,COMP,V,MOD)
{
	T = COMP*V;
	N = size(T)[0];
	for ( I = 0; I < N; I++ )
		if ( T[I] % MOD )
			return 0;
	return modvect(INV*V,MOD);
}

def modmat(A,MOD)
{
	if ( !A )
		return 0;
	S = size(A); N = S[0]; M = S[1];
	MAT = newmat(N,M);
	for ( I = 0, NZ = 0; I < N; I++ )
		for ( J = 0, T1 = A[I], T2 = MAT[I]; J < M; J++ ) {
			T2[J] = T1[J] % MOD;
			NZ = NZ || T2[J];
		}
	return NZ?MAT:0;
}

def modvect(A,MOD)
{
	if ( !A )
		return 0;
	N = size(A)[0];
	VECT = newvect(N);
	for ( I = 0, NZ = 0; I < N; I++ ) {
		VECT[I] = A[I] % MOD;
		NZ = NZ || VECT[I];
	}
	return NZ?VECT:0;
}

def qrmat(A,MOD)
{
	if ( !A )
		return [0,0];
	S = size(A); N = S[0]; M = S[1];
	Q = newmat(N,M); R = newmat(N,M);
	for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ )
		for ( J = 0, TA = A[I], TQ = Q[I], TR = R[I]; J < M; J++ ) {
			L = iqr(TA[J],MOD); TQ[J] = L[0]; TR[J] = L[1];
			NZQ = NZQ || TQ[J]; NZR = NZR || TR[J];
		}
	return [NZQ?Q:0,NZR?R:0];
}

def qrvect(A,MOD)
{
	if ( !A )
		return [0,0];
	N = size(A)[0];
	Q = newvect(N); R = newvect(N);
	for ( I = 0, NZQ = 0, NZR = 0; I < N; I++ ) {
		L = iqr(A[I],MOD); Q[I] = L[0]; R[I] = L[1];
		NZQ = NZQ || Q[I]; NZR = NZR || R[I];
	}
	return [NZQ?Q:0,NZR?R:0];
}

def max_mag(M)
{
	R = size(M)[0];
	U = 1;
	for ( I = 0; I < R; I++ ) {
		A = max_mag_vect(M[I]);
		U = MAX(A,U);
	}
	return U;
}

def max_mag_vect(V)
{
	R = size(V)[0];
	U = 1;
	for ( I = 0; I < R; I++ ) {
		A = dp_mag(V[I]*<<0>>);
		U = MAX(A,U);
	}
	return U;
}

def gsl_check(B,V,S)
{
	N = length(V);
	U = S[N-1]; M = U[1]; D = U[2];
	W = setminus(V,[var(M)]);
	H = uc(); VH = append(W,[H]);
	for ( T = B; T != []; T = cdr(T) ) {
		A = car(T);
		AH = dp_dtop(dp_homo(dp_ptod(A,W)),VH);
		for ( I = 0, Z = S; I < N-1; I++, Z = cdr(Z) ) {
			L = car(Z); AH = ptozp(subst(AH,L[0],L[1]/L[2]));
		}
		AH = ptozp(subst(AH,H,D));
		R = srem(AH,M);
		if ( dp_gr_print() )
			if ( !R )
				print([A,"ok"]);
			else
				print([A,"bad"]);
		if ( R )
			break;
	}
	return T == [] ? 1 : 0;
}

def vs_dim(G,V,O)
{
	HM = hmlist(G,V,O); ZD = zero_dim(HM,V,O);
	if ( ZD ) {
		MB = dp_mbase(map(dp_ptod,HM,V));
		return length(MB);
	} else
		error("vs_dim : ideal is not zero-dimensional!");
}

def dgr(G,V,O)
{
	P = getopt(proc);
	if ( type(P) == -1 )
		return gr(G,V,O);
	P0 = P[0]; P1 = P[1]; P = [P0,P1];
	map(ox_reset,P);
	ox_cmo_rpc(P0,"dp_gr_main",G,V,0,1,O);
	ox_cmo_rpc(P1,"dp_gr_main",G,V,1,1,O);
	map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
	F = ox_select(P);
	R = ox_get(F[0]);
	if ( F[0] == P0 ) {
		Win = "nonhomo";
		Lose = P1;
	} else {
		Win = "homo";
		Lose = P0;
	}
	ox_reset(Lose);
	return [Win,R];
}

/* competitive Gbase computation : F4 vs. Bucbberger */
/* P : process list */

def dgrf4mod(G,V,M,O)
{
	P = getopt(proc);
	if ( type(P) == -1 )
		return dp_f4_mod_main(G,V,M,O);
	P0 = P[0]; P1 = P[1]; P = [P0,P1];
	map(ox_reset,P);
	ox_cmo_rpc(P0,"dp_f4_mod_main",G,V,M,O);
	ox_cmo_rpc(P1,"dp_gr_mod_main",G,V,0,M,O);
	map(ox_push_cmd,P,262); /* 262 = OX_popCMO */
	F = ox_select(P);
	R = ox_get(F[0]);
	if ( F[0] == P0 ) {
		Win = "F4";
		Lose = P1;
	} else {
		Win = "Buchberger";
		Lose = P0;
	}
	ox_reset(Lose);
	return [Win,R];
}

/* functions for rpc */

def register_matrix(M)
{
	REMOTE_MATRIX = M; return 0;
}

def register_nfv(L)
{
	REMOTE_NF = L[0]; REMOTE_VARS = L[1]; return 0;
}

def r_ttob(T,M)
{
	return hen_ttob(T,REMOTE_NF,0,REMOTE_VARS,M);
}

def r_ttob_gsl(L,M)
{
	return cons(L[2],hen_ttob(L[0],REMOTE_NF,L[1],REMOTE_VARS,M));
}

def get_matrix()
{
	REMOTE_MATRIX;
}

extern NFArray$

/*
 * HL = [[c,i,m,d],...]
 * if c != 0
 *   g = 0
 *   g = (c*g + m*gi)/d
 *   ...
 *   finally compare g with NF
 *   if g == NF then NFArray[NFIndex] = g
 *
 * if c = 0 then HL consists of single history [0,i,0,0],
 * which means that dehomogenization of NFArray[i] should be 
 * eqall to NF.
 */

def check_trace(NF,NFIndex,HL)
{
	if ( !car(HL)[0] ) {
		/* dehomogenization */
		DH = dp_dehomo(NFArray[car(HL)[1]]);
		if ( NF == DH ) {
			realloc_NFArray(NFIndex);
			NFArray[NFIndex] = NF;
			return 0;
		} else
			error("check_trace(dehomo)");
	}

	for ( G = 0, T = HL; T != []; T = cdr(T) ) {
		H = car(T);

		Coeff = H[0];
		Index = H[1];
		Monomial = H[2];
		Denominator = H[3];

		Reducer = NFArray[Index];
		G = (Coeff*G+Monomial*Reducer)/Denominator;
	}
	if ( NF == G ) {
		realloc_NFArray(NFIndex);
		NFArray[NFIndex] = NF;
		return 0;
	} else
		error("check_trace");
}

/*
 * Trace = [Input,[[j1,[[c,i,m,d],...]],[j2,[[...],...]],...]]
 * if c != 0
 *   g = 0
 *   g = (c*g + m*gi)/d
 *   ...
 *   finally fj = g
 */

def show_trace(Trace,V)
{
	Input = Trace[0];
	for ( I = 0, T = Input; T != []; T = cdr(T), I++ ) {
		print("F"+rtostr(I)+"=",0);
		print(dp_dtop(car(T),V));
	}
	Trace = cdr(Trace);
	for ( T = Trace; T != []; T = cdr(T) ) {
		HL = car(T);
		J = car(HL); HL = HL[1];
		L = length(HL);
		print("F"+rtostr(J)+"=",0);	
		for ( I = 0; I < L; I++ ) print("(",0);
		for ( First = 1, S = HL; S != []; S = cdr(S) ) {
			H = car(S);
	
			Coeff = H[0];
			Index = H[1];
			Monomial = H[2];
			Denominator = H[3];
			if ( First ) {
				if ( Monomial != 1 ) {
					print("(",0);
					print(type(Monomial)==9?dp_dtop(Monomial,V):Monomial,0);
					print(")*",0);
				}
				print("F"+rtostr(Index)+")",0);	
			} else {
				if ( Coeff != 1 ) {
					print("*(",0); print(Coeff,0); print(")",0);
				}
				print("+",0);
				if ( Monomial != 1 ) {
					print("(",0);
					print(type(Monomial)==9?dp_dtop(Monomial,V):Monomial,0);
					print(")*",0);
				}
				print("F"+rtostr(Index)+")",0);	
				if ( Denominator != 1 ) {
					print("/",0); print(Denominator,0);
				}
			}
			if ( First ) First = 0;
		}
		print("");
	}
}

def generating_relation(Trace,V)
{
	Trace = cdr(Trace);
	Tab = [];
	for ( T = Trace; T != []; T = cdr(T) ) {
		HL = car(T);
		J = car(HL); HL = HL[1];
		L = length(HL);
		LHS = strtov("f"+rtostr(J));
		Dn = 1;
		for ( First = 1, S = HL; S != []; S = cdr(S) ) {
			H = car(S);
	
			Coeff = H[0];
			Index = H[1];
			Monomial = type(H[2])==9?dp_dtop(H[2],V):H[2];
			Denominator = H[3];
			F = strtov("f"+rtostr(Index));
			for ( Z = Tab; Z != []; Z = cdr(Z) )
				if ( Z[0][0] == F ) break;
			if ( Z != [] ) Value = Z[0][1];
			else Value = [F,1];
			if ( First ) {
				RHS = Monomial*Value[0];
				Dn *= Value[1];
			} else {
				RHS = RHS*Coeff*Value[1]+Dn*Value[0]*Monomial;
				Dn = Value[1]*Dn*Denominator;
			}
			VVVV = tttttttt;
			P = ptozp(Dn*VVVV+RHS);
			RHS = coef(P,0,VVVV);
			Dn = coef(P,1,VVVV);
			if ( First ) First = 0;
		}
		Tab = cons([LHS,[RHS,Dn]],Tab);
	}
	return Tab;
}

def generating_relation_mod(Trace,V,M)
{
	Trace = cdr(Trace);
	Tab = [];
	for ( T = Trace; T != []; T = cdr(T) ) {
		HL = car(T);
		J = car(HL); HL = HL[1];
		L = length(HL);
		LHS = strtov("f"+rtostr(J));
		Dn = 1;
		for ( First = 1, S = HL; S != []; S = cdr(S) ) {
			H = car(S);
	
			Coeff = H[0];
			Index = H[1];
			Monomial = type(H[2])==9?dp_dtop(H[2],V):H[2];
			F = strtov("f"+rtostr(Index));
			for ( Z = Tab; Z != []; Z = cdr(Z) )
				if ( Z[0][0] == F ) break;
			if ( Z != [] ) Value = Z[0][1];
			else Value = F;
			if ( First ) {
				RHS = (Monomial*Value)%M;
			} else {
				RHS = ((RHS*Coeff+Value*Monomial)*inv(H[3],M))%M;
			}
			if ( First ) First = 0;
		}
		Tab = cons([LHS,RHS],Tab);
	}
	return Tab;
}
/*
 * realloc NFArray so that it can hold * an element as NFArray[Ind].
 */

def realloc_NFArray(Ind)
{
	if ( Ind == size(NFArray)[0] ) {
		New = newvect(Ind + 100);
		for ( I = 0; I < Ind; I++ )
			New[I] = NFArray[I];
		NFArray = New;
	}
}

/*
 * create NFArray and initialize it by List.
 */

def register_input(List)
{
	Len = length(List);
	NFArray = newvect(Len+100,List);
}

/*
	tracetogen(): preliminary version

	dp_gr_main() returns  [GB,GBIndex,Trace].
	GB : groebner basis
	GBIndex : IndexList (corresponding to Trace)
	Trace : [InputList,Trace0,Trace1,...]
	TraceI : [Index,TraceList]
	TraceList : [[Coef,Index,Monomial,Denominator],...]
	Poly <- 0
	Poly <- (Coef*Poly+Monomial*PolyList[Index])/Denominator
*/

def tracetogen(G)
{
	GB = G[0]; GBIndex = G[1]; Trace = G[2];

	InputList = Trace[0];
	Trace = cdr(Trace);

	/* number of initial basis */
	Nini = length(InputList);

	/* number of generated basis */
	Ngen = length(Trace);
	
	N = Nini + Ngen;

	/* stores traces */
	Tr = vector(N);

	/* stores coeffs */
	Coef = vector(N);

	/* XXX create dp_ptod(1,V) */
	HT = dp_ht(InputList[0]);
	One = dp_subd(HT,HT);

	for ( I = 0; I < Nini; I++ ) {
		Tr[I] = [1,I,One,1];
		C = vector(Nini);
		C[I] = One;
		Coef[I] = C; 
	}
	for ( ; I < N; I++ )
		Tr[I] = Trace[I-Nini][1];

	for ( T = GBIndex; T != []; T = cdr(T) )
		compute_coef_by_trace(car(T),Tr,Coef);
	return Coef;
}

def compute_coef_by_trace(I,Tr,Coef)
{
	if ( Coef[I] )
		return;

	/* XXX */
	Nini = size(Coef[0])[0];

	/* initialize coef vector */
	CI = vector(Nini);

	for ( T = Tr[I]; T != []; T = cdr(T) ) {
		/*	Trace = [Coef,Index,Monomial,Denominator] */
		Trace = car(T);
		C = Trace[0];
		Ind = Trace[1];
		Mon = Trace[2];
		Den = Trace[3];
		if ( !Coef[Ind] )
			compute_coef_by_trace(Ind,Tr,Coef);

		/* XXX */
		CT = newvect(Nini);
		for ( J = 0; J < Nini; J++ )
			CT[J] = (C*CI[J]+Mon*Coef[Ind][J])/Den;
		CI = CT;
	}
	Coef[I] = CI;
}

extern Gbcheck_DP,Gbcheck_IL$

def register_data_for_gbcheck(DPL)
{
	for ( IL = [], I = length(DPL)-1; I >= 0; I-- )
		IL = cons(I,IL);
	Gbcheck_DP = newvect(length(DPL),DPL);
	Gbcheck_IL = IL;
}

def sp_nf_for_gbcheck(Pair)
{
	SP = dp_sp(Gbcheck_DP[Pair[0]],Gbcheck_DP[Pair[1]]);
	return dp_nf(Gbcheck_IL,SP,Gbcheck_DP,1);
}

def gbcheck(B,V,O)
{
	dp_ord(O);
	D = map(dp_ptod,B,V);	
	L = dp_gr_checklist(D,length(V));
	DP = L[0]; Plist = L[1];
	for ( IL = [], I = size(DP)[0]-1; I >= 0; I-- )
		IL = cons(I,IL);
	Procs = getopt(proc);
	if ( type(Procs) == 4 ) {
		map(ox_reset,Procs);
		/* register DP in servers */
		map(ox_cmo_rpc,Procs,"register_data_for_gbcheck",vtol(DP));
		/* discard return value in stack */
		map(ox_pop_cmo,Procs);
		Free = Procs;
		Busy = [];
		T = Plist;
		while ( T != [] || Busy != []  ){
			if ( Free == [] || T == [] ) {
				/* someone is working; wait for data */
				Ready = ox_select(Busy);
				Busy = setminus(Busy,Ready);
				Free = append(Ready,Free);
				for ( ; Ready != []; Ready = cdr(Ready) ) {
					if ( ox_get(car(Ready)) ) {
						map(ox_reset,Procs);
						return 0;
					}
				}
			} else {
				P = car(Free);
				Free = cdr(Free);
				Busy = cons(P,Busy);
				Pair = car(T);
				T = cdr(T);
				ox_cmo_rpc(P,"sp_nf_for_gbcheck",Pair);
				ox_push_cmd(P,262); /* 262 = OX_popCMO */
			}
		}
		map(ox_reset,Procs);
		return 1;
	} else {
		for ( T = Plist; T != []; T = cdr(T) ) {
			Pair = T[0];
			SP = dp_sp(DP[Pair[0]],DP[Pair[1]]);
			if ( dp_nf(IL,SP,DP,1) )
				return 0;
		}
		return 1;
	}
}
end$