===================================================================
RCS file: /home/cvs/OpenXM/src/asir-contrib/testing/noro/Attic/pd.rr,v
retrieving revision 1.6
retrieving revision 1.7
diff -u -p -r1.6 -r1.7
--- OpenXM/src/asir-contrib/testing/noro/Attic/pd.rr	2010/05/21 06:45:06	1.6
+++ OpenXM/src/asir-contrib/testing/noro/Attic/pd.rr	2010/06/02 04:25:46	1.7
@@ -12,15 +12,15 @@ localf partial_decomp, partial_decomp0, zprimacomp, zp
 localf fast_gb, incremental_gb, elim_gb, ldim, make_mod_subst$
 localf rsgn, find_npos, gen_minipoly, indepset$
 localf maxindep, contraction, ideal_list_intersection, ideal_intersection$
-localf radical_membership, quick_radical_membership, modular_radical_membership$
+localf radical_membership, modular_radical_membership$
 localf radical_membership_rep, ideal_product, saturation$
 localf sat, satind, sat_ind, colon$
 localf ideal_colon, ideal_sat, ideal_inclusion, qd_simp_comp, qd_remove_redundant_comp$
-localf pd_remove_redundant_comp, ppart, sq$
-localf lcfactor, compute_deg0, compute_deg, member$
+localf pd_remove_redundant_comp, ppart, sq, gen_fctr, gen_nf, gen_gb_comp$
+localf gen_mptop, lcfactor, compute_deg0, compute_deg, member$
 localf elimination, setintersection, setminus, sep_list$
 localf first_element, comp_tdeg, tdeg, comp_by_ord, comp_by_second$
-localf gbcheck,f4,sathomo,qdcheck$
+localf gbcheck,f4,sathomo,qd_check$
 
 SatHomo=0$
 GBCheck=1$
@@ -71,11 +71,12 @@ def kill_procs()
 
 def qd_check(B,V,QD)
 {
-	G = nd_gr(B,V,0,0);
-	Iso = ideal_list_intersection(map(first_element,QD[0]),V,0);
-	Emb = ideal_list_intersection(map(first_element,QD[1]),V,0);
-	GG = ideal_intersection(Iso,Emb,V,0);
-	return gb_comp(G,GG);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	G = nd_gr(B,V,Mod,0);
+	Iso = ideal_list_intersection(map(first_element,QD[0]),V,0|mod=Mod);
+	Emb = ideal_list_intersection(map(first_element,QD[1]),V,0|mod=Mod);
+	GG = ideal_intersection(Iso,Emb,V,0|mod=Mod);
+	return gen_gb_comp(G,GG,Mod);
 }
 
 /* SYC primary decomositions */
@@ -83,12 +84,13 @@ def qd_check(B,V,QD)
 def syca_dec(B,V)
 {
 T00 = time();
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(Nolexdec=getopt(nolexdec)) == -1 ) Nolexdec = 0;
 	if ( type(SepIdeal=getopt(sepideal)) == -1 ) SepIdeal = 1;
 	if ( type(NoSimp=getopt(nosimp)) == -1 ) NoSimp = 0;
 	if ( type(Time=getopt(time)) == -1 ) Time = 0;
 	Ord = 0;
-	Gt = G0 = G = fast_gb(B,V,0,Ord); 
+	Gt = G0 = G = fast_gb(B,V,Mod,Ord); 
 	Q0 = Q = []; IntQ0 = IntQ = [1]; First = 1;
 	C = 0;
 
@@ -97,56 +99,63 @@ T00 = time();
 	while ( 1 ) {
 		if ( type(Gt[0])==1 ) break;
 		T0 = time();
-		Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec);
+		Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec,mod=Mod);
 		T1 = time(); Tass += T1[0]-T0[0]+T1[1]-T0[1]; Rass += T1[3]-T0[3];
 		T0 = time();
-		Qt = iso_comp(Gt,Pt,V,Ord);
+		Qt = iso_comp(Gt,Pt,V,Ord|mod=Mod,isgb=1);
 		T1 = time(); Tiso += T1[0]-T0[0]+T1[1]-T0[1]; Riso += T1[3]-T0[3];
-		IntQt = ideal_list_intersection(map(first_element,Qt),V,Ord);
-		IntPt = ideal_list_intersection(map(first_element,Pt),V,Ord);
+		IntQt = ideal_list_intersection(map(first_element,Qt),V,Ord|mod=Mod);
+		IntPt = ideal_list_intersection(map(first_element,Pt),V,Ord|mod=Mod);
 		if ( First ) {
 			IntQ0 = IntQ = IntQt; IntP = IntPt; Qi = Qt; First = 0;
 		} else {
-			IntQ1 = ideal_intersection(IntQ,IntQt,V,Ord);
-			if ( gb_comp(IntQ,IntQ1) ) {
+			IntQ1 = ideal_intersection(IntQ,IntQt,V,Ord|mod=Mod);
+			if ( gen_gb_comp(IntQ,IntQ1,Mod) ) {
 				G = Gt; IntP = IntPt; Q = []; IntQ = [1]; C = 0;
 				continue;
 			} else {
 				IntQ = IntQ1; 
-				IntQ1 = ideal_intersection(IntQ0,IntQt,V,Ord);
-				if ( !gb_comp(IntQ0,IntQ1) ) {
+				IntQ1 = ideal_intersection(IntQ0,IntQt,V,Ord|mod=Mod);
+				if ( !gen_gb_comp(IntQ0,IntQ1,Mod) ) {
+					Q = append(Qt,Q); 
+#if 1
+					for ( T = Qt; T != []; T = cdr(T) )
+						if ( !ideal_inclusion(IntQ0,car(T)[0],V,Ord|mod=Mod) )
+							Q0 = append(Q0,[car(T)]);
+#else
+					Q0 = append(Q0,Qt);
+#endif
 					IntQ0 = IntQ1;
-					Q = append(Qt,Q); Q0 = append(Qt,Q0);
 				}
 			}
 		}
-		if ( gb_comp(IntQt,Gt) || gb_comp(IntQ,G) || gb_comp(IntQ0,G0) ) break;
+		if ( gen_gb_comp(IntQt,Gt,Mod) || gen_gb_comp(IntQ,G,Mod) || gen_gb_comp(IntQ0,G0,Mod) ) break;
 		T0 = time();
-		C1 = ideal_colon(G,IntQ,V);
+		C1 = ideal_colon(G,IntQ,V|mod=Mod);
 		T1 = time(); Tcolon += T1[0]-T0[0]+T1[1]-T0[1]; Rcolon += T1[3]-T0[3];
-		if ( C && gb_comp(C,C1) ) {
+		if ( C && gen_gb_comp(C,C1,Mod) ) {
 			G = Gt; IntP = IntPt; Q = []; IntQ = [1]; C = 0;
 			continue;
 		} else C = C1;
 		T0 = time();
 		if ( SepIdeal == 0 )
-			Ok = find_separating_ideal0(C,G,IntQ,IntP,V,Ord);
+			Ok = find_separating_ideal0(C,G,IntQ,IntP,V,Ord|mod=Mod);
 		else if ( SepIdeal == 1 )
-			Ok = find_separating_ideal1(C,G,IntQ,IntP,V,Ord);
+			Ok = find_separating_ideal1(C,G,IntQ,IntP,V,Ord|mod=Mod);
 		else if ( SepIdeal == 2 )
-			Ok = find_separating_ideal2(C,G,IntQ,IntP,V,Ord);
+			Ok = find_separating_ideal2(C,G,IntQ,IntP,V,Ord|mod=Mod);
 		G1 = append(Ok,G);
-		Gt1 = fast_gb(G1,V,0,Ord);
+		Gt1 = fast_gb(G1,V,Mod,Ord);
 		T1 = time(); Tsep += T1[0]-T0[0]+T1[1]-T0[1]; Rsep += T1[3]-T0[3];
 #if 0
-		if ( ideal_inclusion(Gt1,Gt,V,Ord) ) {
+		if ( ideal_inclusion(Gt1,Gt,V,Ord|mod=Mod) ) {
 			G = Gt; IntP = IntPt; Q = []; IntQ = [1]; C = 0;
 		} else
 #endif
 			Gt = Gt1;
 	}		
 	T0 = time();
-	if ( !NoSimp ) Q1 = qd_remove_redundant_comp(G0,Qi,Q0,V,Ord);
+	if ( !NoSimp ) Q1 = qd_remove_redundant_comp(G0,Qi,Q0,V,Ord|mod=Mod);
 	else Q1 = Q0;
 	if ( Time ) {
 		T1 = time(); Tirred += T1[0]-T0[0]+T1[1]-T0[1]; Rirred += T1[3]-T0[3];
@@ -161,50 +170,51 @@ T00 = time();
 def syc_dec(B,V)
 {
 T00 = time();
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(Nolexdec=getopt(nolexdec)) == -1 ) Nolexdec = 0;
 	if ( type(SepIdeal=getopt(sepideal)) == -1 ) SepIdeal = 1;
 	if ( type(NoSimp=getopt(nosimp)) == -1 ) NoSimp = 0;
 	if ( type(Time=getopt(time)) == -1 ) Time = 0;
 	Ord = 0;
-	G = fast_gb(B,V,0,Ord);
+	G = fast_gb(B,V,Mod,Ord);
 	Q = []; IntQ = [1]; Gt = G; First = 1;
 	Tass = Tiso = Tcolon = Tsep = Tirred = 0;
 	Rass = Riso = Rcolon = Rsep = Rirred = 0;
 	while ( 1 ) {
 		if ( type(Gt[0])==1 ) break;
 		T0 = time();
-		Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec);
+		Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec,mod=Mod);
 		T1 = time(); Tass += T1[0]-T0[0]+T1[1]-T0[1]; Rass += T1[3]-T0[3];
 		T0 = time();
-		Qt = iso_comp(Gt,Pt,V,Ord);
+		Qt = iso_comp(Gt,Pt,V,Ord|mod=Mod,isgb=1);
 		T1 = time(); Tiso += T1[0]-T0[0]+T1[1]-T0[1]; Riso += T1[3]-T0[3];
-		IntQt = ideal_list_intersection(map(first_element,Qt),V,Ord);
-		IntPt = ideal_list_intersection(map(first_element,Pt),V,Ord);
+		IntQt = ideal_list_intersection(map(first_element,Qt),V,Ord|mod=Mod);
+		IntPt = ideal_list_intersection(map(first_element,Pt),V,Ord|mod=Mod);
 		if ( First ) {
 			IntQ = IntQt; Qi = Qt; First = 0;
 		} else {
-			IntQ1 = ideal_intersection(IntQ,IntQt,V,Ord);
-			if ( !gb_comp(IntQ1,IntQ) )
+			IntQ1 = ideal_intersection(IntQ,IntQt,V,Ord|mod=Mod);
+			if ( !gen_gb_comp(IntQ1,IntQ,Mod) )
 				Q = append(Qt,Q);
 		}
-		if ( gb_comp(IntQ,G) || gb_comp(IntQt,Gt) ) 
+		if ( gen_gb_comp(IntQ,G,Mod) || gen_gb_comp(IntQt,Gt,Mod) ) 
 			break;
 		T0 = time();
-		C = ideal_colon(Gt,IntQt,V);
+		C = ideal_colon(Gt,IntQt,V|mod=Mod);
 		T1 = time(); Tcolon += T1[0]-T0[0]+T1[1]-T0[1]; Rcolon += T1[3]-T0[3];
 		T0 = time();
 		if ( SepIdeal == 0 )
-			Ok = find_separating_ideal0(C,Gt,IntQt,IntPt,V,Ord);
+			Ok = find_separating_ideal0(C,Gt,IntQt,IntPt,V,Ord|mod=Mod);
 		else if ( SepIdeal == 1 )
-			Ok = find_separating_ideal1(C,Gt,IntQt,IntPt,V,Ord);
+			Ok = find_separating_ideal1(C,Gt,IntQt,IntPt,V,Ord|mod=Mod);
 		else if ( SepIdeal == 2 )
-			Ok = find_separating_ideal2(C,Gt,IntQt,IntPt,V,Ord);
+			Ok = find_separating_ideal2(C,Gt,IntQt,IntPt,V,Ord|mod=Mod);
 		G1 = append(Ok,Gt);
-		Gt = fast_gb(G1,V,0,Ord);
+		Gt = fast_gb(G1,V,Mod,Ord);
 		T1 = time(); Tsep += T1[0]-T0[0]+T1[1]-T0[1]; Rsep += T1[3]-T0[3];
 	}
 	T0 = time();
-	if ( !NoSimp ) Q1 = qd_remove_redundant_comp(G,Qi,Q,V,Ord);
+	if ( !NoSimp ) Q1 = qd_remove_redundant_comp(G,Qi,Q,V,Ord|mod=Mod);
 	else Q1 = Q;
 	T1 = time(); Tirred += T1[0]-T0[0]+T1[1]-T0[1]; Rirred += T1[3]-T0[3];
 	Tall = T1[0]-T00[0]+T1[1]-T00[1]; Rall += T1[3]-T00[3];
@@ -216,52 +226,55 @@ T00 = time();
 	return [Qi,Q1];
 }
 
+/* XXX */
 /* C=G:Q, Rad=rad(Q), return J s.t. Q cap (G+J) = G */
 
 def find_separating_ideal0(C,G,Q,Rad,V,Ord) {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	for ( CI = C, I = 1; ; I++ ) {
 		for ( T = CI, S = []; T != []; T = cdr(T) )
-			if ( nd_nf(car(T),Q,V,Ord,0) ) S = cons(car(T),S);
+			if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
 		if ( S == [] )
 			error("find_separating_ideal0 : cannot happen");
 		G1 = append(S,G);
-		Int = ideal_intersection(G1,Q,V,Ord);
+		Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
 		/* check whether (Q cap (G+S)) = G */
-		if ( gb_comp(Int,G) ) return reverse(S);
-		CI = ideal_product(CI,C,V);
+		if ( gen_gb_comp(Int,G,Mod) ) return reverse(S);
+		CI = ideal_product(CI,C,V|mod=Mod);
 	}
 }
 
 def find_separating_ideal1(C,G,Q,Rad,V,Ord) {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	for ( T = C, S = []; T != []; T = cdr(T) )
-		if ( nd_nf(car(T),Q,V,Ord,0) ) S = cons(car(T),S);
+		if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
 	if ( S == [] )
 		error("find_separating_ideal1 : cannot happen");
 	G1 = append(S,G);
-	Int = ideal_intersection(G1,Q,V,Ord);
+	Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
 	/* check whether (Q cap (G+S)) = G */
-	if ( gb_comp(Int,G) ) return reverse(S);
+	if ( gen_gb_comp(Int,G,Mod) ) return reverse(S);
 
 	/* or qsort(C,comp_tdeg) */
 	C = qsort(S,comp_tdeg);
 
 	Tmp = ttttt; TV = cons(Tmp,V); Ord1 = [[0,1],[Ord,length(V)]];
 	Int0 = incremental_gb(append(vtol(ltov(G)*Tmp),vtol(ltov(Q)*(1-Tmp))),
-		TV,Ord1|gbblock=[[0,length(G)]]);
+		TV,Ord1|gbblock=[[0,length(G)]],mod=Mod);
 	for ( T = C, S = []; T != []; T = cdr(T) ) {
-		if ( !nd_nf(car(T),Rad,V,Ord,0) ) continue;
+		if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
 		Ui = U = car(T);
 		for ( I = 1; ; I++ ) {
 			G1 = cons(Ui,G);
-			Int = ideal_intersection(G1,Q,V,Ord);
-			if ( gb_comp(Int,G) ) break;
+			Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
+			if ( gen_gb_comp(Int,G,Mod) ) break;
 			else
-				Ui = nd_nf(Ui*U,G,V,Ord,0);
+				Ui = gen_nf(Ui*U,G,V,Ord,Mod);
 		}
 		Int1 = incremental_gb(append(Int0,[Tmp*Ui]),TV,Ord1
-			|gbblock=[[0,length(Int0)]]);
+			|gbblock=[[0,length(Int0)]],mod=Mod);
 		Int = elimination(Int1,V);
-		if ( !gb_comp(Int,G) ) 
+		if ( !gen_gb_comp(Int,G,Mod) ) 
 			break;
 		else {
 			Int0 = Int1;
@@ -272,31 +285,34 @@ def find_separating_ideal1(C,G,Q,Rad,V,Ord) {
 }
 
 def find_separating_ideal2(C,G,Q,Rad,V,Ord) {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	for ( T = C, S = []; T != []; T = cdr(T) )
-		if ( nd_nf(car(T),Q,V,Ord,0) ) S = cons(car(T),S);
+		if ( gen_nf(car(T),Q,V,Ord,Mod) ) S = cons(car(T),S);
 	if ( S == [] )
 		error("find_separating_ideal2 : cannot happen");
 	G1 = append(S,G);
-	Int = ideal_intersection(G1,Q,V,Ord);
+	Int = ideal_intersection(G1,Q,V,Ord|mod=Mod);
 	/* check whether (Q cap (G+S)) = G */
-	if ( gb_comp(Int,G) ) return reverse(S);
+	if ( gen_gb_comp(Int,G,Mod) ) return reverse(S);
 
 	/* or qsort(S,comp_tdeg) */
 	C = qsort(C,comp_tdeg);
 	Dp = dp_gr_print(); dp_gr_print(0);
 	for ( T = C, S = []; T != []; T = cdr(T) ) {
-		print(length(T));
-		if ( !nd_nf(car(T),Rad,V,Ord,0) ) continue;
+		if ( !gen_nf(car(T),Rad,V,Ord,Mod) ) continue;
 		Ui = U = car(T);
 		for ( I = 1; ; I++ ) {
-			G1 = cons(Ui,G);
-			Int = ideal_intersection(G1,Q,V,Ord);
-			if ( gb_comp(Int,G) ) break;
+			G1 = append(G,[Ui]);
+			Int = ideal_intersection(G1,Q,V,Ord|mod=Mod,
+				gbblock=[[0,length(G)],[length(G1),length(Q)]]);
+			if ( gen_gb_comp(Int,G,Mod) ) break;
 			else
-				Ui = nd_nf(Ui*U,G,V,Ord,0);
+				Ui = gen_nf(Ui*U,G,V,Ord,Mod);
 		}
+		print([length(T),I],2);
 		S = cons(Ui,S);
 	}
+	print("");
 	S = qsort(S,comp_tdeg);
 	/* S = reverse(S); */
 	Len = length(S);
@@ -307,13 +323,13 @@ def find_separating_ideal2(C,G,Q,Rad,V,Ord) {
 		Cur = 2;
 		G1 = append(G,[S[0]]);
 		Int0 = incremental_gb(append(vtol(ltov(G1)*Tmp),vtol(ltov(Q)*(1-Tmp))),
-			TV,Ord1|gbblock=[[0,length(G)]]);
+			TV,Ord1|gbblock=[[0,length(G)]],mod=Mod);
 		while ( Prev < Cur ) {
 			for ( St = [], I = Prev; I < Cur; I++ ) St = cons(Tmp*S[I],St);
 			Int1 = incremental_gb(append(Int0,St),TV,Ord1
-				|gbblock=[[0,length(Int0)]]);
+				|gbblock=[[0,length(Int0)]],mod=Mod);
 			Int = elimination(Int1,V);
-			if ( gb_comp(Int,G) ) {
+			if ( gen_gb_comp(Int,G,Mod) ) {
 				print(Cur);
 				Prev = Cur;
 				Cur = Cur+idiv(Len-Cur+1,2);
@@ -335,43 +351,46 @@ def find_separating_ideal2(C,G,Q,Rad,V,Ord) {
 
 def sy_dec(B,V)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(Nolexdec=getopt(nolexdec)) == -1 ) Nolexdec = 0;
 	Ord = 0;
-	G = fast_gb(B,V,0,Ord);
+	G = fast_gb(B,V,Mod,Ord);
 	Q = [];
 	IntQ = [1];
 	Gt = G;
 	First = 1;
 	while ( 1 ) {
 		if ( type(Gt[0]) == 1 ) break;
-		Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec);
-		L = pseudo_dec(Gt,Pt,V,Ord);
+		Pt = prime_dec(Gt,V|indep=1,nolexdec=Nolexdec,mod=Mod);
+		L = pseudo_dec(Gt,Pt,V,Ord|mod=Mod);
 		Qt = L[0]; Rt = L[1]; St = L[2];
-		IntQt = ideal_list_intersection(Qt,V,Ord);
+		IntQt = ideal_list_intersection(map(first_element,Qt),V,Ord|mod=Mod);
 		if ( First ) {
 			IntQ = IntQt;
 			Qi = Qt;
 			First = 0;
 		} else {
-			IntQ = ideal_intersection(IntQ,IntQt,V,Ord);
+			IntQ = ideal_intersection(IntQ,IntQt,V,Ord|mod=Mod);
 			Q = append(Qt,Q);
 		}
-		if ( gb_comp(IntQ,G) ) break;
+		if ( gen_gb_comp(IntQ,G,Mod) ) break;
 		for ( T = Rt; T != []; T = cdr(T) ) {
 			if ( type(car(T)[0]) == 1 ) continue;
-			U = sy_dec(car(T),V|nolexdec=Nolexdec);
-			IntQ = ideal_list_intersection(cons(IntQ,U),V,Ord);
+			U = sy_dec(car(T),V|nolexdec=Nolexdec,mod=Mod);
+			IntQ = ideal_list_intersection(cons(IntQ,map(first_element,U)),
+				V,Ord|mod=Mod);
 			Q = append(U,Q);
-			if ( gb_comp(IntQ,G) ) break;
+			if ( gen_gb_comp(IntQ,G,Mod) ) break;
 		}
-		Gt = fast_gb(append(Gt,St),V,0,Ord);
+		Gt = fast_gb(append(Gt,St),V,Mod,Ord);
 	}
-	Q = qd_remove_redundant_comp(G,Qi,Q,V,Ord);
+	Q = qd_remove_redundant_comp(G,Qi,Q,V,Ord|mod=Mod);
 	return append(Qi,Q);
 }
 
 def pseudo_dec(G,L,V,Ord)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	N = length(L);
 	S = vector(N);
 	Q = vector(N);
@@ -379,54 +398,56 @@ def pseudo_dec(G,L,V,Ord)
 	L0 = map(first_element,L);
 	for ( I = 0; I < N; I++ ) {
 		LI = setminus(L0,[L0[I]]);
-		PI = ideal_list_intersection(LI,V,Ord);
+		PI = ideal_list_intersection(LI,V,Ord|mod=Mod);
 		PI = qsort(PI,comp_tdeg);
 		for ( T = PI; T != []; T = cdr(T) )
-			if ( p_nf(car(T),L0[I],V,Ord) ) break;
+			if ( gen_nf(car(T),L0[I],V,Ord,Mod) ) break;
 		if ( T == [] ) error("separator : cannot happen");
-		SI = sat_ind(G,car(T),V);	
+		SI = satind(G,car(T),V|mod=Mod);	
 		QI = SI[0];
 		S[I] = car(T)^SI[1];
 		PV = L[I][1];
 		V0 = setminus(V,PV);
 #if 0
-		GI = fast_gb(QI,append(V0,PV),0,
+		GI = fast_gb(QI,append(V0,PV),Mod,
 			[[Ord,length(V0)],[Ord,length(PV)]]);
 #else
-		GI = fast_gb(QI,append(V0,PV),0,
+		GI = fast_gb(QI,append(V0,PV),Mod,
 			[[2,length(V0)],[Ord,length(PV)]]);
 #endif
-		LCFI = lcfactor(GI,V0,Ord);
+		LCFI = lcfactor(GI,V0,Ord,Mod);
 		for ( F = 1, T = LCFI, Gt = QI; T != []; T = cdr(T) ) {
-			St = sat_ind(Gt,T[0],V);
+			St = satind(Gt,T[0],V|mod=Mod);
 			Gt = St[0]; F *= T[0]^St[1];
 		}
-		Q[I] = Gt;
-		R[I] = fast_gb(cons(F,QI),V,0,Ord);
+		Q[I] = [Gt,L0[I]];
+		R[I] = fast_gb(cons(F,QI),V,Mod,Ord);
 	}
 	return [vtol(Q),vtol(R),vtol(S)];
 }
 
 def iso_comp(G,L,V,Ord)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
 	N = length(L);
 	S = vector(N);
 	Ind = vector(N);
 	Q = vector(N);
 	L0 = map(first_element,L);
-	G = nd_gr(G,V,0,Ord);
+	if ( !IsGB ) G = nd_gr(G,V,Mod,Ord);
 	for ( I = 0; I < N; I++ ) {
 		LI = setminus(L0,[L0[I]]);
-		PI = ideal_list_intersection(LI,V,Ord);
+		PI = ideal_list_intersection(LI,V,Ord|mod=Mod);
 		for ( T = PI; T != []; T = cdr(T) )
-			if ( p_nf(car(T),L0[I],V,Ord) ) break;
+			if ( gen_nf(car(T),L0[I],V,Ord,Mod) ) break;
 		if ( T == [] ) error("separator : cannot happen");
 		S[I] = car(T);
-		QI = sat(G,S[I],V|isgb=1);
+		QI = sat(G,S[I],V|isgb=1,mod=Mod);
 		PV = L[I][1];
 		V0 = setminus(V,PV);
-		GI = elim_gb(QI,V0,PV,0,[[0,length(V0)],[0,length(PV)]]);
-		Q[I] = [contraction(GI,V0),L0[I]];
+		GI = elim_gb(QI,V0,PV,Mod,[[0,length(V0)],[0,length(PV)]]);
+		Q[I] = [contraction(GI,V0|mod=Mod),L0[I]];
 	}
 	return vtol(Q);
 }
@@ -435,71 +456,77 @@ def iso_comp(G,L,V,Ord)
 
 def prima_dec(B,V)
 {
-	G = nd_gr_trace(B,V,1,GBCheck,0);
-	G0 = G;
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	if ( type(Ord=getopt(ord)) == -1 ) Ord = 0;
+	G0 = fast_gb(B,V,Mod,0);
+	G = fast_gb(G0,V,Mod,Ord);
 	IntP = [1];
 	QD = [];
 	while ( 1 ) {
-		if ( ideal_inclusion(IntP,G0,V,0) ) 
-			return QD;
-		W = maxindep(G,V,0); NP = length(W);
+		if ( type(G[0])==1 || ideal_inclusion(IntP,G0,V,0|mod=Mod) ) 
+			break;
+		W = maxindep(G,V,Ord); NP = length(W);
 		V0 = setminus(V,W); N = length(V0);
 		V1 = append(V0,W);
-		G1 = fast_gb(G,V1,0,[[0,N],[0,NP]]);
-		LCF = lcfactor(G1,V0,0);
-		L = zprimacomp(G,V0);
+		G1 = fast_gb(G,V1,Mod,[[Ord,N],[Ord,NP]]);
+		LCF = lcfactor(G1,V0,Ord,Mod);
+		L = zprimacomp(G,V0|mod=Mod);
 		F = 1;
-		for ( T = LCF, G2 = G1; T != []; T = cdr(T) ) {
-			S = sat_ind(G2,T[0],V1);
+		for ( T = LCF, G2 = G; T != []; T = cdr(T) ) {
+			S = satind(G2,T[0],V1|mod=Mod);
 			G2 = S[0]; F *= T[0]^S[1];
 		}
 		for ( T = L, QL = []; T != []; T = cdr(T) )
 			QL = cons(car(T)[0],QL);
-		Int = ideal_list_intersection(QL,V,0);
-		IntP = ideal_intersection(IntP,Int,V,0);
+		Int = ideal_list_intersection(QL,V,0|mod=Mod);
+		IntP = ideal_intersection(IntP,Int,V,0|mod=Mod);
 		QD = append(QD,L);
-		F = p_nf(F,G,V,0);
-		G = cons(F,G);
+		F = gen_nf(F,G,V,0,Mod);
+		G = fast_gb(cons(F,G),V,Mod,Ord);
 	}
+	QD = qd_remove_redundant_comp(G0,[],QD,V,0);
+	return QD;
 }
 
 /* SL prime decomposition */
 
 def prime_dec(B,V)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
 	if ( type(NoLexDec=getopt(nolexdec)) == -1 ) NoLexDec = 0;
-	B = map(sq,B);
+	B = map(sq,B,Mod);
 	if ( !NoLexDec )
-		PD = lex_predec1(B,V);
+		PD = lex_predec1(B,V|mod=Mod);
 	else
 		PD = [B];
-	G = ideal_list_intersection(PD,V,0);
-	PD = pd_remove_redundant_comp(G,PD,V,0);
+	G = ideal_list_intersection(PD,V,0|mod=Mod);
+	PD = pd_remove_redundant_comp(G,PD,V,0|mod=Mod);
 	R = [];
 	for ( T = PD; T != []; T = cdr(T) )
-		R = append(prime_dec_main(car(T),V|indep=Indep),R);
+		R = append(prime_dec_main(car(T),V|indep=Indep,mod=Mod),R);
 	if ( Indep ) {
-		G = ideal_list_intersection(map(first_element,R),V,0);
-		if ( !NoLexDec ) R = pd_remove_redundant_comp(G,R,V,0|first=1);
+		G = ideal_list_intersection(map(first_element,R),V,0|mod=Mod);
+		if ( !NoLexDec ) R = pd_remove_redundant_comp(G,R,V,0|first=1,mod=Mod);
 	} else {
-		G = ideal_list_intersection(R,V,0);
-		if ( !NoLexDec ) R = pd_remove_redundant_comp(G,R,V,0);
+		G = ideal_list_intersection(R,V,0|mod=Mod);
+		if ( !NoLexDec ) R = pd_remove_redundant_comp(G,R,V,0|mod=Mod);
 	}
 	return R;
 }
 
 def prime_dec_main(B,V)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(Indep=getopt(indep)) == -1 ) Indep = 0;
-	G = nd_gr_trace(B,V,1,GBCheck,0);
+	G = fast_gb(B,V,Mod,0);
 	IntP = [1];
 	PD = [];
 	while ( 1 ) {
 		/* rad(G) subset IntP */		
 		/* check if IntP subset rad(G) */
 		for ( T = IntP; T != []; T = cdr(T) ) {
-			if ( (GNV = modular_radical_membership(car(T),G,V)) ) {
+			if ( (GNV = modular_radical_membership(car(T),G,V|mod=Mod)) ) {
 				F = car(T);
 				break;
 			}
@@ -507,20 +534,20 @@ def prime_dec_main(B,V)
 		if ( T == [] ) return PD;
 
 		/* GNV = [GB(<NV*F-1,G>),NV] */
-		G1 = nd_gr_trace(GNV[0],cons(GNV[1],V),1,GBCheck,[[0,1],[0,length(V)]]);
+		G1 = fast_gb(GNV[0],cons(GNV[1],V),Mod,[[0,1],[0,length(V)]]);
 		G0 = elimination(G1,V);
-		PD0 = zprimecomp(G0,V,Indep);
+		PD0 = zprimecomp(G0,V,Indep|mod=Mod);
 		if ( Indep ) {
-			Int = ideal_list_intersection(PD0[0],V,0);
+			Int = ideal_list_intersection(PD0[0],V,0|mod=Mod);
 			IndepSet = PD0[1];	
 			for ( PD1 = [], T = PD0[0]; T != []; T = cdr(T) )
 				PD1 = cons([car(T),IndepSet],PD1);
 			PD = append(PD,reverse(PD1));
 		} else {
-			Int = ideal_list_intersection(PD0,V,0);
+			Int = ideal_list_intersection(PD0,V,0|mod=Mod);
 			PD = append(PD,PD0);
 		}
-		IntP = ideal_intersection(IntP,Int,V,0);
+		IntP = ideal_intersection(IntP,Int,V,0|mod=Mod);
 	}
 }
 
@@ -528,14 +555,15 @@ def prime_dec_main(B,V)
 
 def lex_predec1(B,V)
 {
-	G = nd_gr_trace(B,V,1,GBCheck,2);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	G = fast_gb(B,V,Mod,2);
 	for ( T = G; T != []; T = cdr(T) ) {
-		F = fctr(car(T));
+		F = gen_fctr(car(T),Mod);
 		if ( length(F) > 2 || length(F) == 2 && F[1][1] > 1 ) {
 			for ( R = [], S = cdr(F); S != []; S = cdr(S) ) {
 				Ft = car(S)[0];
-				Gt = map(ptozp,map(p_nf,G,[Ft],V,0));
-				R1 = nd_gr_trace(cons(Ft,Gt),V,1,GBCheck,0);
+				Gt = map(ptozp,map(gen_nf,G,[Ft],V,0,Mod));
+				R1 = fast_gb(cons(Ft,Gt),V,Mod,0);
 				R = cons(R1,R);
 			}
 			return R;
@@ -749,7 +777,7 @@ def partial_decomp0(GD,V,PV,Ord,I,Mod)
 		Mt = car(car(T));
 		D1 = D*1;
 		D1[I] = Mt;
-		GIt = map(p_nf,GI,[Mt],V,Ord);
+		GIt = map(gen_nf,GI,[Mt],V,Ord,Mod);
 		G1 = cons(Mt,GIt);
 		Gelim = elim_gb(G1,V,PV,Mod,Ord);
 		D1[N] = LD = ldim(Gelim,V);
@@ -766,34 +794,36 @@ def partial_decomp0(GD,V,PV,Ord,I,Mod)
 /* prime/primary components over rational function field */
 
 def zprimacomp(G,V) {
-	L = zprimadec(G,V,0);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	L = zprimadec(G,V,0|mod=Mod);
 	R = [];
 	dp_ord(0);
 	for ( T = L; T != []; T = cdr(T) ) {
 		S = car(T);
-		UQ = contraction(S[0],V);
-		UP = contraction(S[1],V);
+		UQ = contraction(S[0],V|mod=Mod);
+		UP = contraction(S[1],V|mod=Mod);
 		R = cons([UQ,UP],R);
 	}
 	return R;
 } 
 
 def zprimecomp(G,V,Indep) {
-	W = maxindep(G,V,0);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	W = maxindep(G,V,0|mod=Mod);
 	V0 = setminus(V,W);
 	V1 = append(V0,W);
 #if 0
 	O1 = [[0,length(V0)],[0,length(W)]];
-	G1 = nd_gr_trace(G,V1,1,GBCheck,O1);
+	G1 = fast_gb(G,V1,Mod,O1);
 	dp_ord(0);
 #else
 	G1 = G;
 #endif
-	PD = zprimedec(G1,V0,0);
+	PD = zprimedec(G1,V0,Mod);
 	dp_ord(0);
 	R = [];
 	for ( T = PD; T != []; T = cdr(T) ) {
-		U = contraction(car(T),V0);
+		U = contraction(car(T),V0|mod=Mod);
 		R = cons(U,R);
 	}
 	if ( Indep ) return [R,W];
@@ -802,15 +832,16 @@ def zprimecomp(G,V,Indep) {
 
 def fast_gb(B,V,Mod,Ord)
 {
-	NoRA = (NoRA=getopt(nora))&&type(NoRA)!=-1 ? 1 : 0;
+	if ( type(Block=getopt(gbblock)) == -1 ) Block = 0;
+	if ( type(NoRA=getopt(nora)) == -1 ) NoRA = 0;
 	if ( Mod )
 		G = nd_f4(B,V,Mod,Ord|nora=NoRA);
-	else {
-		if ( F4 )
-			G = map(ptozp,f4_chrem(B,V,Ord));
-		else 
-			G = nd_gr_trace(B,V,1,GBCheck,Ord|nora=NoRA);
-	}
+	else if ( F4 )
+		G = map(ptozp,f4_chrem(B,V,Ord));
+	else if ( Block )
+		G = nd_gr_trace(B,V,1,GBCheck,Ord|nora=NoRA,gbblock=Block);
+	else
+		G = nd_gr_trace(B,V,1,GBCheck,Ord|nora=NoRA);
 	return G;
 }
 
@@ -818,9 +849,12 @@ def incremental_gb(A,V,Ord)
 {
 	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(Block=getopt(gbblock)) == -1 ) Block = 0;
-	if ( Mod )
-		G = nd_gr(A,V,Mod,Ord);
-	else if ( Procs ) {
+	if ( Mod ) {
+		if ( Block )
+			G = nd_gr(A,V,Mod,Ord|gbblock=Block);
+		else
+			G = nd_gr(A,V,Mod,Ord);
+	} else if ( Procs ) {
 		Arg0 = ["nd_gr",A,V,0,Ord];
 		Arg1 = ["nd_gr_trace",A,V,1,GBCheck,Ord];
 		G = competitive_exec(Procs,Arg0,Arg1);
@@ -837,9 +871,12 @@ def elim_gb(G,V,PV,Mod,Ord)
 	O1 = [[0,N],[0,PN]];
 	if ( Ord == O1 )
 		Ord = Ord[0][0];
-	if ( Mod ) /* XXX */
+	if ( Mod ) /* XXX */ {
+		for ( T = G, H = []; T != []; T = cdr(T) )
+			if ( car(T) ) H = cons(car(T),H);
+		G = reverse(H);
 		G = dp_gr_mod_main(G,V,0,Mod,Ord);
-	else if ( Procs ) {
+	} else if ( Procs ) {
 		Arg0 = ["nd_gr_trace",G,V,1,GBCheck,Ord];
 		Arg1 = ["nd_gr_trace_rat",G,V,PV,1,GBCheck,O1,Ord];
 		G = competitive_exec(Procs,Arg0,Arg1);
@@ -856,6 +893,8 @@ def ldim(G,V)
 	return D;
 }
 
+/* over Q only */
+
 def make_mod_subst(GD,V,PV,HC)
 {
 	N = length(V);
@@ -937,7 +976,8 @@ def gen_minipoly(G,V,PV,Ord,VI,Mod)
 	}
 #elif 1
 	if ( Mod ) {
-		G = nd_gr(G,V1,Mod,[[0,length(W)],[0,length(PV1)]]|nora=1);
+		V1 = append(W,PV1);
+		G = nd_gr(G,V1,Mod,[[0,length(W)],[0,length(PV1)]]);
 		G = elimination(G,PV1);
 	} else {
 		PV2 = setminus(PV1,[PV1[length(PV1)-1]]);
@@ -981,7 +1021,8 @@ def indepset(V,H)
 
 def maxindep(B,V,O)
 {
-	G = nd_gr_trace(B,V,1,GBCheck,O);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	G = fast_gb(B,V,Mod,O);
 	Old = dp_ord();
 	dp_ord(O);
 	H = map(dp_dtop,map(dp_ht,map(dp_ptod,G,V)),V);
@@ -1004,10 +1045,11 @@ def maxindep(B,V,O)
 /* ideal operations */
 def contraction(G,V)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	C = [];
 	for ( T = G; T != []; T = cdr(T) ) {
 		C1 = dp_hc(dp_ptod(car(T),V));
-		S = fctr(C1);
+		S = gen_fctr(C1,Mod);
 		for ( S = cdr(S); S != []; S = cdr(S) )
 			if ( !member(S[0][0],C) ) C = cons(S[0][0],C);
 	}
@@ -1017,7 +1059,7 @@ def contraction(G,V)
 	NV = ttttt;
 	for ( T = C, S = 1; T != []; T = cdr(T) )
 		S *= car(T);
-	G = saturation([G,NV],S,W);
+	G = saturation([G,NV],S,W|mod=Mod);
 	return G;
 }
 
@@ -1041,10 +1083,14 @@ def ideal_intersection(A,B,V,Ord)
 	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(Block=getopt(gbblock)) == -1 ) Block = 0;
 	T = ttttt;
-	if ( Mod )
-		G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
-			cons(T,V),Mod,[[0,1],[Ord,length(V)]]);
-	else
+	if ( Mod ) {
+		if ( Block )
+			G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
+				cons(T,V),Mod,[[0,1],[Ord,length(V)]]|gbblock=Block,nora=1);
+		else
+			G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
+				cons(T,V),Mod,[[0,1],[Ord,length(V)]]|nora=1);
+	} else
 	if ( Procs ) {
 		Arg0 = ["nd_gr",
 			append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
@@ -1056,37 +1102,35 @@ def ideal_intersection(A,B,V,Ord)
 	} else {
 		if ( Block )
 			G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
-				cons(T,V),0,[[0,1],[Ord,length(V)]]|gbblock=Block);
+				cons(T,V),0,[[0,1],[Ord,length(V)]]|gbblock=Block,nora=0);
 		else
 			G = nd_gr(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
-				cons(T,V),0,[[0,1],[Ord,length(V)]]);
+				cons(T,V),0,[[0,1],[Ord,length(V)]]|nora=0);
 	}
 	G0 = elimination(G,V);
+	if ( 0 && !Procs )
+		G0 = nd_gr_postproc(G0,V,Mod,Ord,0);
 	return G0;
 }
 
 /* returns GB if F notin rad(G) */
 
 def radical_membership(F,G,V) {
-	F = p_nf(F,G,V,0);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	F = gen_nf(F,G,V,0,Mod);
 	if ( !F ) return 0;
 	NV = ttttt;
-	T = nd_gr_trace(cons(NV*F-1,G),cons(NV,V),1,GBCheck,0);
+	T = fast_gb(cons(NV*F-1,G),cons(NV,V),Mod,0);
 	if ( type(car(T)) != 1 ) return [T,NV];
 	else return 0;
 }
 
-def quick_radical_membership(F,G,V) {
-	F = p_nf(F,G,V,0);
-	if ( !F ) return 1;
-	NV = ttttt;
-	T = nd_f4(cons(NV*F-1,G),cons(NV,V),lprime(0),0);
-	if ( type(car(T)) != 1 ) return 0;
-	else return 1;
-}
-
 def modular_radical_membership(F,G,V) {
-	F = p_nf(F,G,V,0);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	if ( Mod )
+		return radical_membership(F,G,V|mod=Mod);
+
+	F = gen_nf(F,G,V,0,0);
 	if ( !F ) return 0;
 	NV = ttttt;
 	for ( J = 0; ; J++ ) {
@@ -1109,7 +1153,7 @@ def radical_membership_rep(F,G,V,Max,Ord,Mod) {
 	Ft = F;
 	I = 1;
 	while ( Max < 0 || I <= Max ) {
-		Ft = nd_nf(Ft,G,V,Ord,Mod);
+		Ft = gen_nf(Ft,G,V,Ord,Mod);
 		if ( !Ft ) return I;
 		Ft *= F;
 		I++;
@@ -1119,6 +1163,7 @@ def radical_membership_rep(F,G,V,Max,Ord,Mod) {
 
 def ideal_product(A,B,V)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	dp_ord(0);
 	DA = map(dp_ptod,A,V);
 	DB = map(dp_ptod,B,V);
@@ -1139,20 +1184,23 @@ def ideal_product(A,B,V)
 	Len = length(A)>length(B)?length(A):length(B);
 	Len *= 2;
 	L = sep_list(T,Len); B0 = L[0]; B1 = L[1];
-	R = nd_gr_trace(B0,V,0,-1,0);
+	R = fast_gb(B0,V,Mod,0);
 	while ( B1 != [] ) {
 		print(length(B1));
 		L = sep_list(B1,Len);
 		B0 = L[0]; B1 = L[1];
-		R = nd_gr_trace(append(R,B0),V,0,-1,0|gbblock=[[0,length(R)]],nora=1);
+		R = fast_gb(append(R,B0),V,Mod,0|gbblock=[[0,length(R)]],nora=1);
 	}
 	return R;
 }
 
 def saturation(GNV,F,V) 
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	G = GNV[0]; NV = GNV[1];
-	if ( Procs ) {
+	if ( Mod )
+		G1 = nd_gr(cons(NV*F-1,G),cons(NV,V),Mod,[[0,1],[0,length(V)]]);
+	else if ( Procs ) {
 		Arg0 = ["nd_gr_trace", 
 		cons(NV*F-1,G),cons(NV,V),0,GBCheck,[[0,1],[0,length(V)]]];
 		Arg1 = ["nd_gr_trace", 
@@ -1165,9 +1213,12 @@ def saturation(GNV,F,V) 
 
 def sat(G,F,V) 
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
 	NV = ttttt;
-	if ( Procs ) {
+	if ( Mod )
+		G1 = nd_gr(cons(NV*F-1,G),cons(NV,V),Mod,[[0,1],[0,length(V)]]);
+	else if ( Procs ) {
 		Arg0 = ["nd_gr_trace", 
 		cons(NV*F-1,G),cons(NV,V),0,GBCheck,[[0,1],[0,length(V)]]];
 		Arg1 = ["nd_gr_trace", 
@@ -1188,12 +1239,25 @@ def sat(G,F,V) 
 
 def satind(G,F,V)
 {
+	if ( type(Block=getopt(gbblock)) == -1 ) Block = 0;
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	NV = ttttt;
 	N = length(V);
 	B = append(G,[NV*F-1]);
 	V1 = cons(NV,V); 
-	D = nd_gr_trace(B,V1,1,GBCheck,[[0,1],[0,N]]
-		|nora=1,gentrace=1,gbblock=[[0,length(G)]]);
+	Ord1 = [[0,1],[0,N]];
+	if ( Mod )
+		if ( Block )
+			D = nd_gr(B,V1,Mod,Ord1|nora=1,gentrace=1,gbblock=Block);
+		else 
+			D = nd_gr(B,V1,Mod,Ord1|nora=1,gentrace=1);
+	else
+		if ( Block )
+			D = nd_gr_trace(B,V1,SatHomo,GBCheck,Ord1
+				|nora=1,gentrace=1,gbblock=Block);
+		else
+			D = nd_gr_trace(B,V1,SatHomo,GBCheck,Ord1
+				|nora=1,gentrace=1);
 	G1 = D[0];
 	Len = length(G1);
 	Deg = compute_deg(B,V1,NV,D);
@@ -1212,11 +1276,13 @@ def satind(G,F,V)
 
 def sat_ind(G,F,V)
 {
+	if ( type(Ord=getopt(ord)) == -1 ) Ord = 0;
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	NV = ttttt;
-	F = p_nf(F,G,V,0);
+	F = gen_nf(F,G,V,Ord,Mod);
 	for ( I = 0, GI = G; ; I++ ) {
-		G1 = colon(GI,F,V);
-		if ( ideal_inclusion(G1,GI,V,0) )  {
+		G1 = colon(GI,F,V|mod=Mod,ord=Ord);
+		if ( ideal_inclusion(G1,GI,V,Ord|mod=Mod) )  {
 			return [GI,I];
 		}
 		else GI = G1;
@@ -1225,42 +1291,43 @@ def sat_ind(G,F,V)
 
 def colon(G,F,V)
 {
+	if ( type(Ord=getopt(ord)) == -1 ) Ord = 0;
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(IsGB=getopt(isgb)) == -1 ) IsGB = 0;
-	F = p_nf(F,G,V,0);
+	F = gen_nf(F,G,V,Ord,Mod);
 	if ( !F ) return [1];
 	if ( IsGB )
-		T = ideal_intersection(G,[F],V,0|gbblock=[[0,length(G)]]);
+		T = ideal_intersection(G,[F],V,Ord|gbblock=[[0,length(G)]],mod=Mod);
 	else
-		T = ideal_intersection(G,[F],V,0);
-	return map(ptozp,map(sdiv,T,F));
+		T = ideal_intersection(G,[F],V,Ord|mod=Mod);
+	return Mod?map(sdivm,T,F,Mod):map(ptozp,map(sdiv,T,F));
 }
 
 def ideal_colon(G,F,V)
 {
-	G = nd_gr(G,V,0,0);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	G = nd_gr(G,V,Mod,0);
 	for ( T = F, L = []; T != []; T = cdr(T) )
-		L = cons(colon(G,car(T),V|isgb=1),L);
+		L = cons(colon(G,car(T),V|isgb=1,mod=Mod),L);
 	L = reverse(L);
-	return ideal_list_intersection(L,V,0);
+	return ideal_list_intersection(L,V,0|mod=Mod);
 }
 
 def ideal_sat(G,F,V)
 {
-	G = nd_gr(G,V,0,0);
-	L = mapat(sat,1,G,F,V);
-	return ideal_list_intersection(L,V,0);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	G = nd_gr(G,V,Mod,0);
+	for ( T = F, L = []; T != []; T = cdr(T) )
+		L = cons(sat(G,car(T),V|mod=Mod),L);
+	L = reverse(L);
+	return ideal_list_intersection(L,V,0|mod=Mod);
 }
 
 def ideal_inclusion(F,G,V,O)
 {
 	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
-	if ( Mod ) {
-		for ( T = F; T != []; T = cdr(T) )
-			if ( p_nf_mod(car(T),G,V,O,Mod) ) return 0;
-	} else {
-		for ( T = F; T != []; T = cdr(T) )
-			if ( p_nf(car(T),G,V,O) ) return 0;
-	}
+	for ( T = F; T != []; T = cdr(T) )
+		if ( gen_nf(car(T),G,V,O,Mod) ) return 0;
 	return 1;
 }
 
@@ -1268,14 +1335,15 @@ def ideal_inclusion(F,G,V,O)
 
 def qd_simp_comp(QP,V)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	R = ltov(QP);
 	N = length(R);
 	for ( I = 0; I < N; I++ ) {
 		if ( R[I] ) {
 			QI = R[I][0]; PI = R[I][1];
 			for ( J = I+1; J < N; J++ )
-				if ( R[J] && gb_comp(PI,R[J][1]) ) {
-					QI = ideal_intersection(QI,R[J][0],V,0);
+				if ( R[J] && gen_gb_comp(PI,R[J][1],Mod) ) {
+					QI = ideal_intersection(QI,R[J][0],V,0|mod=Mod);
 					R[J] = 0;
 				}
 			R[I] = [QI,PI];
@@ -1288,18 +1356,20 @@ def qd_simp_comp(QP,V)
 
 def qd_remove_redundant_comp(G,Iso,Emb,V,Ord)
 {
-	IsoInt = ideal_list_intersection(map(first_element,Iso),V,Ord);
-	Emb = qd_simp_comp(Emb,V);
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
+	IsoInt = ideal_list_intersection(map(first_element,Iso),V,Ord|mod=Mod);
+	Emb = qd_simp_comp(Emb,V|mod=Mod);
 	Emb = reverse(qsort(Emb));
 	A = ltov(Emb); N = length(A);
 	Pre = IsoInt; Post = vector(N+1);
-	for ( Post[N] = [1], I = N-1; I >= 1; I-- )
-		Post[I] = ideal_intersection(Post[I+1],A[I][0],V,Ord);
+	for ( Post[N] = IsoInt, I = N-1; I >= 1; I-- )
+		Post[I] = ideal_intersection(Post[I+1],A[I][0],V,Ord|mod=Mod);
 	for ( I = 0; I < N; I++ ) {
-		Int = ideal_intersection(Pre,Post[I+1],V,Ord);
-		if ( gb_comp(Int,G) ) A[I] = 0;
+		print(".",2);
+		Int = ideal_intersection(Pre,Post[I+1],V,Ord|mod=Mod);
+		if ( gen_gb_comp(Int,G,Mod) ) A[I] = 0;
 		else
-			Pre = ideal_intersection(Pre,A[I][0],V,Ord);
+			Pre = ideal_intersection(Pre,A[I][0],V,Ord|mod=Mod);
 	}
 	for ( T = [], I = 0; I < N; I++ )
 		if ( A[I] ) T = cons(A[I],T);
@@ -1308,24 +1378,25 @@ def qd_remove_redundant_comp(G,Iso,Emb,V,Ord)
 
 def pd_remove_redundant_comp(G,P,V,Ord)
 {
+	if ( type(Mod=getopt(mod)) == -1 ) Mod = 0;
 	if ( type(First=getopt(first)) == -1 ) First = 0;
 	A = ltov(P); N = length(A);
 	for ( I = 0; I < N; I++ ) {
 		if ( !A[I] ) continue;
 		for ( J = I+1; J < N; J++ )
 			if ( A[J] && 
-				gb_comp(First?A[I][0]:A[I],First?A[J][0]:A[J]) ) A[J] = 0;
+				gen_gb_comp(First?A[I][0]:A[I],First?A[J][0]:A[J],Mod) ) A[J] = 0;
 	}
 	for ( I = 0, T = []; I < N; I++ ) if ( A[I] ) T = cons(A[I],T);
 	A = ltov(reverse(T)); N = length(A);
 	Pre = [1]; Post = vector(N+1);
 	for ( Post[N] = [1], I = N-1; I >= 1; I-- )
-		Post[I] = ideal_intersection(Post[I+1],First?A[I][0]:A[I],V,Ord);
+		Post[I] = ideal_intersection(Post[I+1],First?A[I][0]:A[I],V,Ord|mod=Mod);
 	for ( I = 0; I < N; I++ ) {
-		Int = ideal_intersection(Pre,Post[I+1],V,Ord);
-		if ( gb_comp(Int,G) ) A[I] = 0;
+		Int = ideal_intersection(Pre,Post[I+1],V,Ord|mod=Mod);
+		if ( gen_gb_comp(Int,G,Mod) ) A[I] = 0;
 		else
-			Pre = ideal_intersection(Pre,First?A[I][0]:A[I],V,Ord);
+			Pre = ideal_intersection(Pre,First?A[I][0]:A[I],V,Ord|mod=Mod);
 	}
 	for ( T = [], I = 0; I < N; I++ ) if ( A[I] ) T = cons(A[I],T);
 	return reverse(T);
@@ -1343,22 +1414,22 @@ def ppart(F,V,Mod)
 }
 
 
-def sq(F)
+def sq(F,Mod)
 {
 	if ( !F ) return 0;
-	A = cdr(fctr(F));
+	A = cdr(gen_fctr(F,Mod));
 	for ( R = 1; A != []; A = cdr(A) )
 		R *= car(car(A));
 	return R;
 }
 
-def lcfactor(G,V,O)
+def lcfactor(G,V,O,Mod)
 {
 	O0 = dp_ord(); dp_ord(O);
 	C = [];
 	for ( T = G; T != []; T = cdr(T) ) {
 		C1 = dp_hc(dp_ptod(car(T),V));
-		S = fctr(C1);
+		S = gen_fctr(C1,Mod);
 		for ( S = cdr(S); S != []; S = cdr(S) )
 			if ( !member(S[0][0],C) ) C = cons(S[0][0],C);
 	}
@@ -1366,6 +1437,44 @@ def lcfactor(G,V,O)
 	return C;
 }
 
+def gen_fctr(F,Mod)
+{
+	if ( Mod ) return modfctr(F,Mod);
+	else return fctr(F);
+}
+
+def gen_mptop(F)
+{
+	if ( !F ) return F;
+	else if ( type(F)==1 )
+		if ( ntype(F)==5 ) return mptop(F);
+		else return F;
+	else {
+		V = var(F);
+		D = deg(F,V);
+		for ( R = 0, I = 0; I <= D; I++ )
+			if ( C = coef(F,I,V) ) R += gen_mptop(C)*V^I;
+		return R;
+	}
+}
+
+def gen_nf(F,G,V,Ord,Mod)
+{
+	if ( !Mod ) return p_nf(F,G,V,Ord);
+
+	setmod(Mod);
+	dp_ord(Ord); DF = dp_mod(dp_ptod(F,V),Mod,[]);
+	N = length(G); DG = newvect(N);
+	for ( I = N-1, IL = []; I >= 0; I-- ) {
+		DG[I] = dp_mod(dp_ptod(G[I],V),Mod,[]);
+		IL = cons(I,IL);
+	}
+	T = dp_nf_mod(IL,DF,DG,1,Mod);
+	for ( R = 0; T; T = dp_rest(T) )
+		R += gen_mptop(dp_hc(T))*dp_dtop(dp_ht(T),V);
+	return R;
+}
+
 /* Ti = [D,I,M,C] */
 
 def compute_deg0(Ti,P,V,TV)
@@ -1504,5 +1613,16 @@ def comp_by_second(A,B)
 	else if ( A[1] < B[1] ) return -1;
 	else return 0;
 }
+
+def gen_gb_comp(A,B,Mod)
+{
+	if ( !Mod ) return gb_comp(A,B);
+	LA = length(A); LB = length(B);
+	if ( LA != LB ) return 0;
+	A = qsort(A); B = qsort(B);
+	if ( A != B ) return 0;
+	return 1;
+}
+
 endmodule$
 end$