===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/PU.c,v
retrieving revision 1.11
retrieving revision 1.15
diff -u -p -r1.11 -r1.15
--- OpenXM_contrib2/asir2000/engine/PU.c	2004/09/14 07:23:34	1.11
+++ OpenXM_contrib2/asir2000/engine/PU.c	2018/02/28 03:55:38	1.15
@@ -45,7 +45,7 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/engine/PU.c,v 1.10 2004/06/25 03:07:51 noro Exp $ 
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/PU.c,v 1.14 2010/01/31 03:25:54 noro Exp $ 
 */
 #include "ca.h"
 
@@ -114,6 +114,57 @@ void substp(VL vl,P p,V v0,P p0,P *pr)
 	}
 }
 
+void substpp(VL vl,P p,V *vvect,P *svect,int nv,P *pr);
+
+void substpp(VL vl,P p,V *vvect,P *svect,int nv,P *pr)
+{
+	P x,t,m,c,s,a,p0,c1;
+	DCP dc;
+	Q d;
+	V v;
+	int i;
+
+	if ( !p ) 
+		*pr = 0;
+	else if ( NUM(p) ) 
+		*pr = p;
+	else {
+		v = VR(p);
+		for ( i = 0; i < nv; i++ ) if ( vvect[i] == v ) break;
+		if ( svect[i] && OID(svect[i]) < 0 ) {
+			MKV(VR(p),x);
+			for ( c = 0, dc = DC(p); dc; dc = NEXT(dc) ) {
+				substpp(vl,COEF(dc),vvect,svect,nv,&t);
+				if ( DEG(dc) ) {
+					pwrp(vl,x,DEG(dc),&s); mulp(vl,s,t,&m); 
+					addp(vl,m,c,&a); 
+					c = a;
+				} else {
+					addp(vl,t,c,&a); 
+					c = a;
+				}
+			}
+			*pr = c;
+		} else {
+			p0 = svect[i];
+			dc = DC(p);
+			substpp(vl,COEF(dc),vvect,svect,nv,&c);
+			for ( d = DEG(dc), dc = NEXT(dc); 
+				dc; d = DEG(dc), dc = NEXT(dc) ) {
+					subq(d,DEG(dc),(Q *)&t); pwrp(vl,p0,(Q)t,&s); 
+					mulp(vl,s,c,&m); 
+					substpp(vl,COEF(dc),vvect,svect,nv,&c1);
+					addp(vl,m,c1,&c); 
+			}
+			if ( d ) {
+				pwrp(vl,p0,d,&t); mulp(vl,t,c,&m); 
+				c = m;
+			}
+			*pr = c;
+		}
+	}
+}
+
 void detp(VL vl,P **rmat,int n,P *dp)
 {
 	int i,j,k,l,sgn,nmin,kmin,lmin,ntmp;
@@ -131,7 +182,7 @@ void detp(VL vl,P **rmat,int n,P *dp)
 			*dp = 0; return;
 		}
 		nmin = nmonop(mat[i][j]);
-		kmin=lmin=j;
+		kmin=i; lmin=j;
 		for ( k = j; k < n; k++ )
 			for ( l = j; l < n; l++ )
 				if ( mat[k][l] && ((ntmp=nmonop(mat[k][l])) < nmin) ) {
@@ -678,6 +729,10 @@ void premp(VL vl,P p1,P p2,P *pr)
 		*pr = 0;
 	else if ( ( v1 = VR(p1) ) == ( v2 = VR(p2) ) ) {
 		n1 = deg(v1,p1); n2 = deg(v1,p2);
+    if ( n1 < n2 ) {
+      *pr = p1;
+      return;
+    }
 		pw = (P *)ALLOCA((n1+1)*sizeof(P));
 		bzero((char *)pw,(int)((n1+1)*sizeof(P)));