===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/sp,v
retrieving revision 1.1
retrieving revision 1.5
diff -u -p -r1.1 -r1.5
--- OpenXM_contrib2/asir2000/lib/sp	1999/12/03 07:39:11	1.1
+++ OpenXM_contrib2/asir2000/lib/sp	2000/04/20 02:20:16	1.5
@@ -1,10 +1,11 @@
-/* $OpenXM: OpenXM/src/asir99/lib/sp,v 1.1.1.1 1999/11/10 08:12:30 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/lib/sp,v 1.1.1.1 1999/12/03 07:39:11 noro Exp $ */
 /* 
 	sp : functions related to algebraic number fields
 
 	Revision History:
 
-	99/08/24    noro    modified for 1999 release version
+	2000/04/20    noro    fixed bugs around gathering algebraic numbers
+	1999/08/24    noro    modified for 1999 release version
 */
 
 #include "defs.h"
@@ -84,7 +85,7 @@ def flatmf(L) {
 def af(P,AL)
 {
 	RESTIME=UFTIME=GCDTIME=SQTIME=0;
-	S = reverse(asq(P,AL));
+	S = reverse(asq(P));
 	for ( L = []; S != []; S = cdr(S) ) {
 		FM = car(S); F = FM[0]; M = FM[1];
 		G = af_sp(F,AL,1);
@@ -124,7 +125,7 @@ def af_spmain(P,AL,INIT,HINT,PP,SHIFT)
 	N = simpcoef(sp_norm(A0,V,subst(P,V,V-INIT*A0),AL));
 	RESTIME+=time()[0]-TTT;
 	TTT = time()[0];
-	DCSQ = sortfs(asq(N,AL));
+	DCSQ = sortfs(asq(N));
 	SQTIME+=time()[0]-TTT;
 	for ( G = P, A = V+INIT*A0, DCR = []; DCSQ != []; DCSQ = cdr(DCSQ) ) {
 		C = TT(DCSQ); D = TS(DCSQ);
@@ -142,7 +143,7 @@ def af_spmain(P,AL,INIT,HINT,PP,SHIFT)
 			else {
 				S = subst(car(DCT),V,A);
 				if ( pra(G,S,AL) )
-					U = cr_gcda(S,G,AL);
+					U = cr_gcda(S,G);
 				else
 					U = S;
 			}
@@ -270,26 +271,40 @@ def getalgp(P) {
 	else {
 		for ( V = var(P), I = deg(P,V), T = []; I >= 0; I-- )
 			if ( C = coef(P,I) )
-				T = union(T,getalgp(C));
+				T = union_sort(T,getalgp(C));
 		return T;
 	}
 }
 		
-def union(A,B)
-{
-	for ( T = B; T != []; T = cdr(T) )
-		A = union1(A,car(T));
-	return A;
+def getalgtreep(P) {
+	if ( type(P) <= 1 )
+		return getalgtree(P);
+	else {
+		for ( V = var(P), I = deg(P,V), T = []; I >= 0; I-- )
+			if ( C = coef(P,I) )
+				T = union_sort(T,getalgtreep(C));
+		return T;
+	}
 }
 
-def union1(A,E)
+/* C = union of A and B; A and B is sorted. C should also be sorted. */
+
+def union_sort(A,B)
 {
 	if ( A == [] )
-		return [E];
-	else if ( car(A) == E )
+		return B;
+	else if ( B == [] )
 		return A;
-	else
-		return cons(car(A),union1(cdr(A),E));
+	else {
+		A0 = car(A);
+		B0 = car(B);
+		if ( A0 == B0 )
+			return cons(A0,union_sort(cdr(A),cdr(B)));
+		else if ( A0 > B0 )
+			return cons(A0,union_sort(cdr(A),B));
+		else
+			return cons(B0,union_sort(A,cdr(B)));
+	}
 }
 
 def invalgp(A)
@@ -359,7 +374,7 @@ def sortfs(L)
 
 def pdiva(P1,P2)
 {
-	A = union(getalgp(P1),getalgp(P2));
+	A = union_sort(getalgp(P1),getalgp(P2));
 	P1 = algptorat(P1); P2 = algptorat(P2);
 	return simpalg(rattoalgp(sdiv(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2),A));
 }
@@ -371,7 +386,7 @@ def pqra(P1,P2)
 	else if ( (type(P1) != POLY) || (deg(P1,var(P1)) < deg(P2,var(P1))) )
 		return [0,P1];
 	else {
-		A = union(getalgp(P1),getalgp(P2));
+		A = union_sort(getalgp(P1),getalgp(P2));
 		P1 = algptorat(P1); P2 = algptorat(P2);
 		L = sqr(P1*LCOEF(P2)^(DEG(P1)-DEG(P2)+1),P2);
 		return [simpalg(rattoalgp(L[0],A)),simpalg(rattoalgp(L[1],A))];
@@ -410,7 +425,7 @@ def sort_alg(VL)
 	return L;
 }
 
-def asq(P,AL)
+def asq(P)
 {
 	P = simpalg(P);
 	if ( type(P) == NUM )
@@ -423,7 +438,7 @@ def asq(P,AL)
 			if ( type(F) == NUM ) 
 				break;
 			F1 = diff(F,V);
-			GCD = cr_gcda(F,F1,AL);
+			GCD = cr_gcda(F,F1);
 			FLAT = pdiva(F,GCD);
 			if ( type(GCD) == NUM ) {
 				A[I] = F; B[I] = 1; 
@@ -479,7 +494,7 @@ def ufctrhint_heuristic(P,HINT,PP,SHIFT) {
 		G = igcd(G,DT=deg(T,V));
 	if ( G == 1 ) {
 		if ( K*deg(PPP,V) != deg(P,V) )
-			PPP = cr_gcda(PPP,P,AL);
+			PPP = cr_gcda(PPP,P);
 		return ufctrhint2(P,HINT,PPP,AL);
 	} else {
 		for ( S = 0, T = P; T; T -= coef(T,DT)*V^DT ) {
@@ -489,7 +504,7 @@ def ufctrhint_heuristic(P,HINT,PP,SHIFT) {
 		L = fctr(S);
 		for ( DC = [car(L)], L = cdr(L); L != []; L = cdr(L) ) {
 			H = subst(car(car(L)),V,V^G);
-			HH = cr_gcda(PPP,H,AL);
+			HH = cr_gcda(PPP,H);
 			T = ufctrhint2(H,HINT,HH,AL);
 			DC = append(DC,T);
 		}
@@ -645,27 +660,13 @@ def norm_ch_lag(V,VM,P,P0) {
 	return S;
 }
 
-def cr_gcda(P1,P2,EXT)
+def cr_gcda(P1,P2)
 {
 	if ( !(V = var(P1)) || !var(P2) )
 		return 1;
-	AL = union(getalgp(P1),getalgp(P2));
-	if ( AL == [] )
+	EXT = union_sort(getalgtreep(P1),getalgtreep(P2));
+	if ( EXT == [] )
 		return gcd(P1,P2);
-	T = newvect(length(EXT));
-	for ( TAL = AL; TAL != []; TAL = cdr(TAL) ) {
-		A = getalg(car(TAL));
-		for ( TA = A; TA != []; TA = cdr(TA) ) {
-			B = car(TA);
-			for ( TEXT = EXT, I = 0; TEXT != []; TEXT = cdr(TEXT), I++ )
-				if ( car(TEXT) == B )
-					T[I] = B;
-		}
-	}
-	for ( I = length(EXT)-1, S = []; I >= 0; I-- )
-		if ( T[I] )
-			S = cons(T[I],S);
-	EXT = S;
 	NEXT = length(EXT);
 	if ( deg(P1,V) < deg(P2,V) ) {
 		T = P1; P1 = P2; P2 = T;
@@ -745,18 +746,6 @@ def member(E,L)
 		if ( E == car(L) )
 			return 1;
 	return 0;
-}
-
-def getallalg(A)
-{
-	T = cdr(vars(defpoly(A)));
-	if ( T == [] )
-		return [A];
-	else {
-		for ( S = [A]; T != []; T = cdr(T) )
-			S = union(S,getallalg(rattoalg(car(T))));
-		return S;
-	}
 }
 
 def discr(P) {