/* $OpenXM: OpenXM_contrib2/asir2000/lib/gr,v 1.1 1999/12/03 07:39:11 noro Exp $ */
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)
{
	G = dp_gr_main(B,V,0,1,O);
	return G;
}

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

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

def hgr_mod(B,V,O,M)
{
	G = dp_gr_mod_main(B,V,1,M,O);
	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)
{
	TM = TE = TNF = 0;
	N = length(V); HM = hmlist(G0,V,O); ZD = zero_dim(HM,V,O);
	if ( !ZD )
		error("tolex : ideal is not zero-dimensional!");
	MB = dp_mbase(map(dp_ptod,HM,V));
	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];
		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 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 ) {
			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);
		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);
	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)
{
	DIM = length(MB);
	DV = newvect(DIM);
	for ( T = GM, SL = [], LCM = 1; T != []; T = cdr(T) ) {
		S = p_terms(car(T),V,2);
		dp_ord(O); RHS = termstomat(NF,map(dp_ptod,cdr(S),V),MB,M);
		dp_ord(0); 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);
		print(["DN",B[1]]);
	}
	return SL;
}

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!");

	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 R[0];
	}
}

/* subroutines */

def gennf(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(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);
		}
		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 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);
	}
}

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)];
	}
	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);
		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);
		IL = cons(I,IL);
	}
	L = dp_true_nf(IL,DP,DB,1);
	return [dp_dtop(L[0],V),L[1]];
}

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)
{
	for ( T = A; T != []; T = cdr(T) ) {
		for ( S = B, M = car(T), N = -M; S != []; S = cdr(S) )
			if ( car(S) == M || car(S) == N )
				break;
		if ( S == [] )
			break;
	}
	return T == [] ? 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);
	for ( I = 0, C = BB, X = 0, 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)
{
	P0 = P[0]; P1 = P[1]; P = [P0,P1];
	flush(P0); flush(P1);
	rpc(P0,"dp_gr_main",G,V,0,1,O);
	rpc(P1,"dp_gr_main",G,V,1,1,O);
	F = select(P);
	R = rpcrecv(F[0]); flush(P0); flush(P1);
	return 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;
}
end$