===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/al.c,v
retrieving revision 1.1.1.1
retrieving revision 1.5
diff -u -p -r1.1.1.1 -r1.5
--- OpenXM_contrib2/asir2000/builtin/al.c	1999/12/03 07:39:07	1.1.1.1
+++ OpenXM_contrib2/asir2000/builtin/al.c	2001/10/09 01:36:04	1.5
@@ -1,10 +1,10 @@
-/* $OpenXM: OpenXM/src/asir99/builtin/al.c,v 1.2 1999/11/18 05:42:01 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/builtin/al.c,v 1.4 2001/03/09 01:14:13 noro Exp $ */
 /* ----------------------------------------------------------------------
-   $Id: al.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $
+   $Id: al.c,v 1.5 2001/10/09 01:36:04 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
    ----------------------------------------------------------------------
@@ -14,6 +14,7 @@
 #include <parse.h>
 #include <al.h>
 
+void Preverse();
 void Phugo();
 void Pex();
 void Pall();
@@ -23,6 +24,7 @@ void Pfargs();
 void Pfopargs();
 void Pcompf();
 void Patnum();
+int gauss_abc();
 int compf();
 void Patl();
 void Pqevar();
@@ -38,6 +40,8 @@ void simpl_th2atl();
 int simpl_gand_udnargls();
 int simpl_gand_thupd();
 int simpl_gand_thprsism();
+int simpl_gand_smtbdlhs();
+int simpl_gand_smtbelhs();
 void lbc();
 void replaceq();
 void deleteq();
@@ -61,6 +65,7 @@ void qeblock_verbose0();
 int getmodulus();
 int qevar();
 int gausselim();
+int delv();
 int translate();
 int translate_a();
 void translate_a1();
@@ -70,11 +75,12 @@ void mkqgp();
 void getqcoeffs();
 void mkdiscr();
 int al_reorder();
-int indices();
+void indices();
 void mkeset();
 int selectside();
 int cmp2n();
 void add2eset();
+void seproots();
 void sp_add2eset();
 void subgpf();
 void subref();
@@ -116,12 +122,15 @@ void rep();
 void gpp();
 void esetp();
 void nodep();
+void gauss_mkeset1();
+void gauss_mkeset2();
 
 extern Verbose;
 
 struct oRE {
 	P p;
 	P discr;
+        int rootno;
 	int itype;
 };
 
@@ -161,9 +170,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
@@ -444,7 +454,7 @@ F f,*pnf;
 NODE th;
 int n;
 {
-	F h,hh;
+	F h;
 	oFOP op=FOP(f);
 
 	if (AL_ATOMIC(op)) {
@@ -566,7 +576,6 @@ NODE th,*patl,*patlc;
 int n;
 {
   NODE atl=NULL,atlc=NULL;
-  LBF h;
   F at,negat;
 
   switch (gand) {
@@ -1250,7 +1259,8 @@ int *pleft,*pmodulus;
 	int i=0;
 	
 	if (!Verbose)
-		return;
+	/* added by noro */
+		return 0;
 	if (*pleft == 0) {
 		for (; cvl; cvl=NEXT(cvl))
 			i++;
@@ -1352,7 +1362,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;
@@ -1393,7 +1403,7 @@ Q deg;
 	return (NZNUMBER(*pa) || NZNUMBER(*pb) || NZNUMBER(*pc));
 }
 
-gauss_mkeset1(rlhs,b,peset)
+void gauss_mkeset1(rlhs,b,peset)
 P rlhs,b;
 NODE *peset;
 {
@@ -1403,15 +1413,12 @@ NODE *peset;
 	MKNODE(*peset,hgp,NULL);
 }
 
-gauss_mkeset2(rlhs,a,b,c,peset)
+void gauss_mkeset2(rlhs,a,b,c,peset)
 P rlhs,a,b,c;
 NODE *peset;
 {
-	RE hre;
-	F hf;
 	GP hgp;
-	P discr;
-	NODE esetc;
+	NODE esetc=NULL;
 
 	*peset = NULL;
 	if (!NUM(a)) {
@@ -1420,25 +1427,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;
 }
 
@@ -1448,10 +1463,7 @@ V x;
 NODE trans[];
 {
 	NODE sc,transc[8];
-	RE hre;
-	GP hgp;
 	int bt,w=0;
-	P h;
 
 	for (bt=BTMIN; bt<=BTMAX; bt++)
 		trans[bt] = NULL;
@@ -1482,6 +1494,8 @@ NODE trans[],transc[];
 		return 1;
 	};
 	error("degree violation in translate_a");
+	/* XXX : NOTREACHED */
+	return -1;
 }
 
 void translate_a1(op,mp,trans,transc)
@@ -1508,7 +1522,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,22 +1543,23 @@ 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;
 {
 	P discr;
 	RE hre;
 	F hf;
-	NODE n=NULL,nc;
+	NODE n=NULL,nc=NULL;
 	
 	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;
@@ -1600,7 +1615,7 @@ V v;
 		return 0;
 }
 
-int indices(op,s,pit,pbt)
+void indices(op,s,pit,pbt)
 oFOP op;
 int s,*pit,*pbt;
 {
@@ -1656,7 +1671,7 @@ void mkeset(trans,x,peset)
 NODE trans[],*peset;
 V x;
 {
-	NODE esetc;
+	NODE esetc=NULL;
 	P h;
 	RE hre;
 	GP hgp;
@@ -1670,13 +1685,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 +1748,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;
@@ -1762,7 +1809,7 @@ F f,*pnf;
 V v;
 GP gp;
 {
-	NODE argl=NULL,arglc;
+	NODE argl=NULL,arglc=NULL;
 
 	NEXTNODE(argl,arglc);
 	BDY(arglc) = (pointer)GUARD(gp);
@@ -1812,7 +1859,7 @@ V v;
 RE re;
 {
 	VL no;
-	P rlhs,prem,bdn,nlhs;
+	P rlhs,prem,nlhs;
 	Q dd,dndeg;
 	
 	reordvar(CO,v,&no);
@@ -1844,18 +1891,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)
@@ -1895,7 +1937,7 @@ P a,b,c,d;
 F *pf;
 {
 	F hf;
-	NODE cj=NULL,cjc;
+	NODE cj=NULL,cjc=NULL;
 	P hp1,hp2;
 
 	NEXTNODE(cj,cjc);
@@ -1917,7 +1959,7 @@ P a,b,c,d;
 F *pf;
 {
 	F hf;
-	NODE cj=NULL,cjc,dj=NULL,djc;
+	NODE cj=NULL,cjc=NULL,dj=NULL,djc=NULL;
 	P hp1,hp2;
 
 	NEXTNODE(dj,djc);
@@ -1953,7 +1995,7 @@ P a,b,c,d;
 F *pf;
 {
 	F hf,hf0;
-	NODE cj=NULL,cjc,d1=NULL,d1c,d2=NULL,d2c;
+	NODE cj=NULL,cjc=NULL,d1=NULL,d1c=NULL,d2=NULL,d2c=NULL;
 	P hp1,hp2;
 
 	NEXTNODE(d1,d1c);
@@ -1994,7 +2036,7 @@ P prem,*pa,*pb,*pc,*pld;
 Q fdeg;
 RE re;
 {
-	P a,b,c,alpha,beta,h1,h2,h3;
+	P a,b,c,alpha,beta,h1,h2;
 	Q two;
 
 	alpha = COEF(DC(prem));
@@ -2042,7 +2084,7 @@ F *pnf;
 {
 	P an;
 	F h;
-	NODE c=NULL,cc,d=NULL,dc;
+	NODE c=NULL,cc=NULL,d=NULL,dc=NULL;
 
 	if (lhsdcp == 0) {
 		MKAF(*pnf,op,0);
@@ -2113,7 +2155,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);
 }
 
@@ -2126,7 +2168,7 @@ F *pnf;
 {
 	Q deg;
 	F h;
-	NODE c=NULL,cc,d=NULL,dc;
+	NODE c=NULL,cc=NULL,d=NULL,dc=NULL;
 	P df;
 	
 	degp(v,lhs,&deg);
@@ -2337,7 +2379,7 @@ F f,*pf;
 int neg,disj;
 {
 	F h;
-	NODE sc,nargl=NULL,narglc;
+	NODE sc,nargl=NULL,narglc=NULL;
 	oFOP op=FOP(f);
 
 	if (AL_ATOMIC(op) || AL_TVAL(op)) {