===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/al.c,v
retrieving revision 1.1.1.1
retrieving revision 1.4
diff -u -p -r1.1.1.1 -r1.4
--- OpenXM_contrib2/asir2000/builtin/al.c	1999/12/03 07:39:07	1.1.1.1
+++ OpenXM_contrib2/asir2000/builtin/al.c	2001/03/09 01:14:13	1.4
@@ -1,10 +1,10 @@
-/* $OpenXM: OpenXM/src/asir99/builtin/al.c,v 1.2 1999/11/18 05:42:01 noro Exp $ */
+/* $OpenXM$ */
 /* ----------------------------------------------------------------------
-   $Id: al.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $
+   $Id: al.c,v 1.4 2001/03/09 01:14:13 noro Exp $
    ----------------------------------------------------------------------
    File al.c: Real quantifier elimination code for RISA/ASIR
 
-   Copyright (c) 1996, 1997, 1998 by
+   Copyright (c) 1996-2001 by
    Andreas Dolzmann and Thomas Sturm, University of Passau, Germany
    dolzmann@uni-passau.de, sturm@uni-passau.de
    ----------------------------------------------------------------------
@@ -61,6 +61,7 @@ void qeblock_verbose0();
 int getmodulus();
 int qevar();
 int gausselim();
+int delv();
 int translate();
 int translate_a();
 void translate_a1();
@@ -75,6 +76,7 @@ void mkeset();
 int selectside();
 int cmp2n();
 void add2eset();
+void seproots();
 void sp_add2eset();
 void subgpf();
 void subref();
@@ -122,6 +124,7 @@ extern Verbose;
 struct oRE {
 	P p;
 	P discr;
+        int rootno;
 	int itype;
 };
 
@@ -161,9 +164,10 @@ typedef struct oQVL *QVL;
 #define MKGP(x,g,p) (NEWGP(x),GUARD(x)=g,POINT(x)=p)
 
 #define NEWRE(x) ((x)=(RE)MALLOC(sizeof(struct oRE)))
-#define MKRE(x,pp,d,i) \
-(NEWRE(x),(x)->p=pp,(x)->discr=d,(x)->itype=i)
+#define MKRE(x,pp,d,rn,i) \
+(NEWRE(x),(x)->p=pp,(x)->discr=d,(x)->rootno=rn,(x)->itype=i)
 #define DISC(re) ((re)->discr)
+#define ROOTNO(re) ((re)->rootno)
 #define ITYPE(re) ((re)->itype)
 
 #define STD 0
@@ -1352,7 +1356,7 @@ NODE *peset;
 				if (!w)
 					continue;
 				*px = v;
-				delvip(v,pvl);
+				delv(v,*pvl,pvl);
 				if (a) {
 					gauss_mkeset2(rlhs,a,b,c,peset);
 					return 2;
@@ -1420,25 +1424,33 @@ NODE *peset;
 		BDY(esetc) = (pointer)hgp;
 	}
 	NEXTNODE(*peset,esetc);
-	mkqgp(rlhs,a,b,c,STD,&hgp);
+	mkqgp(rlhs,a,b,c,1,STD,&hgp);
 	BDY(esetc) = (pointer)hgp;
+	NEXTNODE(*peset,esetc);
+	mkqgp(rlhs,a,b,c,2,STD,&hgp);
+	BDY(esetc) = (pointer)hgp;
 }
 
-int delvip(v,pvl)
+int delv(v,vl,pnvl)
 V v;
-VL *pvl;
+VL vl,*pnvl;
 {
-	VL prev;
+        VL nvl=NULL,nvlc;
 
-	if (v == VR(*pvl)) {
-		*pvl = NEXT(*pvl);
+	if (v == VR(vl)) {
+	        *pnvl = NEXT(vl);
 		return 1;
 	}
-	for (prev=*pvl; NEXT(prev); prev=NEXT(prev))
-		if (VR(NEXT(prev)) == v) {
-			NEXT(prev) = NEXT(NEXT(prev));
-			return 1;
-		}
+	for (; vl && (VR(vl) != v); vl=NEXT(vl)) {
+                NEXTVL(nvl,nvlc);
+		VR(nvlc) = VR(vl);
+	}
+	if (vl) {
+                NEXT(nvlc) = NEXT(vl);
+		*pnvl = nvl;
+		return 1;
+	}
+	*pnvl = nvl;
 	return 0;
 }
 
@@ -1508,7 +1520,7 @@ GP *pgp;
 	RE hre;
 	F hf;
 	
-	MKRE(hre,mp,(P)ONE,itype);
+	MKRE(hre,mp,(P)ONE,1,itype);
 	MKAF(hf,AL_NEQ,b);
 	MKGP(*pgp,hf,hre);
 }
@@ -1529,12 +1541,13 @@ NODE trans[],transc[];
 	}
 	indices(op,0,&itype,&btype);
 	NEXTNODE(trans[btype],transc[btype]);
-	mkqgp(mp,a,b,c,itype,&hgp);
+	mkqgp(mp,a,b,c,-1,itype,&hgp);
 	BDY(transc[btype]) = (pointer)hgp;
 }
 
-void mkqgp(mp,a,b,c,itype,pgp)
+void mkqgp(mp,a,b,c,rootno,itype,pgp)
 P mp,a,b,c;
+int rootno;
 int itype;
 GP *pgp;
 {
@@ -1544,7 +1557,7 @@ GP *pgp;
 	NODE n=NULL,nc;
 	
 	mkdiscr(a,b,c,&discr);
-	MKRE(hre,mp,discr,itype);
+	MKRE(hre,mp,discr,rootno,itype);
 	NEXTNODE(n,nc);
 	MKAF(hf,AL_NEQ,a);
 	BDY(nc) = (pointer)hf;
@@ -1670,13 +1683,13 @@ V x;
 		add2eset(trans[cs],peset,&esetc);
 		sp_add2eset(trans[BTSO],deps,peset,&esetc);
 		NEXTNODE(*peset,esetc);
-		MKRE(hre,0,0,dinf);
+		MKRE(hre,0,0,0,dinf);
 		MKGP(hgp,F_TRUE,hre);
 		BDY(esetc) = (pointer)hgp;
 	} else {
 		NEXTNODE(*peset,esetc);
 		MKV(x,h);
-		MKRE(hre,h,(P)ONE,STD);
+		MKRE(hre,h,(P)ONE,1,STD);
 		MKGP(hgp,F_TRUE,hre);
 		BDY(esetc) = (pointer)hgp;
 	}
@@ -1733,16 +1746,48 @@ NODE n1a,n1b,n2a,n2b;
 void add2eset(trfield,peset,pesetc)
 NODE trfield,*peset,*pesetc;
 {
+        NODE ntrfield,ntrfieldc;
+
 	if (trfield == NULL)
 		return;
-	if (*peset == NULL)
-		*peset = *pesetc = trfield;
-	else
-		NEXT(*pesetc) = trfield;
-	while (NEXT(*pesetc))
-		*pesetc = NEXT(*pesetc);
+	seproots(trfield,&ntrfield,&ntrfieldc);
+	if (*peset == NULL) {
+		*peset = ntrfield;
+		*pesetc = ntrfieldc;
+	} else {
+		NEXT(*pesetc) = ntrfield;
+		*pesetc = ntrfieldc;
+	}
 }
 
+void seproots(trfield,pntrfield,pntrfieldc)
+NODE trfield,*pntrfield,*pntrfieldc;
+{
+      	NODE sc;
+	NODE ntrf=NULL,ntrfc;
+	RE hre,hre2;
+	GP hgp,hgp2;
+
+	for (sc=trfield; sc; sc=NEXT(sc)) {
+                hgp = (GP)BDY(sc);
+		hre = POINT(hgp);
+                if (ROOTNO(hre) == -1) {
+		        NEXTNODE(ntrf,ntrfc);
+			MKRE(hre2,PL(hre),DISC(hre),1,ITYPE(hre));
+                        MKGP(hgp2,GUARD(hgp),hre2);
+			BDY(ntrfc) = (pointer)hgp2;
+		        NEXTNODE(ntrf,ntrfc);
+			ROOTNO(hre) = 2;
+			BDY(ntrfc) = (pointer)hgp;
+		} else {
+		        NEXTNODE(ntrf,ntrfc);
+			BDY(ntrfc) = (pointer)hgp;
+		}
+	}
+	*pntrfield = ntrf;
+	*pntrfieldc = ntrfc;
+}
+
 void sp_add2eset(trfield,itype,peset,pesetc)
 NODE trfield,*peset,*pesetc;
 int itype;
@@ -1844,18 +1889,13 @@ RE re;
 F *pf;
 {
 	P a,b,c,ld;
-	F hf;
-	NODE d=NULL,dc;
-	
+
 	getrecoeffs(prem,fdeg,re,&a,&b,&c,&ld);
-	NEXTNODE(d,dc);
-	substd_a21(op,a,b,c,ld,&hf);
-	BDY(dc) = (pointer)hf;
-	NEXTNODE(d,dc);
-	chsgnp(b,&b);
-	substd_a21(op,a,b,c,ld,&hf);
-	BDY(dc) = (pointer)hf;
-	MKJF(*pf,AL_OR,d);
+	if (ROOTNO(re) == 1)
+	  chsgnp(b,&b);
+	else if (ROOTNO(re) != 2)
+	  error("unspecified quadratic root in substd_a2");
+	substd_a21(op,a,b,c,ld,pf);
 }
 
 void substd_a21(op,a,b,c,d,pf)
@@ -2113,7 +2153,7 @@ RE r;
 	RE stdre;
 	
 	subpme_a_o1(FOP(af),FPL(af),v,ITYPE(r)==MEPS,&h);
-	MKRE(stdre,PL(r),DISC(r),STD);
+	MKRE(stdre,PL(r),DISC(r),ROOTNO(r),STD);
 	subref(h,v,stdre,pnf);
 }