===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/F.c,v
retrieving revision 1.1
retrieving revision 1.13
diff -u -p -r1.1 -r1.13
--- OpenXM_contrib2/asir2000/engine/F.c	1999/12/03 07:39:08	1.1
+++ OpenXM_contrib2/asir2000/engine/F.c	2018/03/29 01:32:51	1.13
@@ -1,1191 +1,1294 @@
-/* $OpenXM: OpenXM/src/asir99/engine/F.c,v 1.1.1.1 1999/11/10 08:12:26 noro Exp $ */
+/*
+ * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED 
+ * All rights reserved.
+ * 
+ * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
+ * non-exclusive and royalty-free license to use, copy, modify and
+ * redistribute, solely for non-commercial and non-profit purposes, the
+ * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
+ * conditions of this Agreement. For the avoidance of doubt, you acquire
+ * only a limited right to use the SOFTWARE hereunder, and FLL or any
+ * third party developer retains all rights, including but not limited to
+ * copyrights, in and to the SOFTWARE.
+ * 
+ * (1) FLL does not grant you a license in any way for commercial
+ * purposes. You may use the SOFTWARE only for non-commercial and
+ * non-profit purposes only, such as academic, research and internal
+ * business use.
+ * (2) The SOFTWARE is protected by the Copyright Law of Japan and
+ * international copyright treaties. If you make copies of the SOFTWARE,
+ * with or without modification, as permitted hereunder, you shall affix
+ * to all such copies of the SOFTWARE the above copyright notice.
+ * (3) An explicit reference to this SOFTWARE and its copyright owner
+ * shall be made on your publication or presentation in any form of the
+ * results obtained by use of the SOFTWARE.
+ * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
+ * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
+ * for such modification or the source code of the modified part of the
+ * SOFTWARE.
+ * 
+ * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
+ * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
+ * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
+ * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
+ * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
+ * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
+ * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
+ * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
+ * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
+ * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
+ * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
+ * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
+ * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
+ * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
+ * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
+ * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
+ * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
+ *
+ * $OpenXM: OpenXM_contrib2/asir2000/engine/F.c,v 1.12 2013/01/08 07:25:58 noro Exp $ 
+*/
 #include "ca.h"
 #include <math.h>
 
-void homfctr();
+int use_new_hensel;
 
-void fctrp(vl,f,dcp)
-VL vl;
-P f;
-DCP *dcp;
+void fctrp(VL vl,P f,DCP *dcp)
 {
-	VL nvl;
-	DCP dc;
+  VL nvl;
+  DCP dc;
 
-	if ( !f || NUM(f) ) {
-		NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
-		NEXT(dc) = 0; *dcp = dc;
-		return;
-	} else if ( !qpcheck((Obj)f) )
-		error("fctrp : invalid argument");	
-	else {
-		clctv(vl,f,&nvl);
-		if ( !NEXT(nvl) ) 
-			ufctr(f,1,dcp);
-		else
-			mfctr(nvl,f,dcp);
-	}
+  if ( !f || NUM(f) ) {
+    NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
+    NEXT(dc) = 0; *dcp = dc;
+    return;
+  } else if ( !qpcheck((Obj)f) )
+    error("fctrp : invalid argument");  
+  else {
+    clctv(vl,f,&nvl);
+    if ( !NEXT(nvl) ) 
+      ufctr(f,1,dcp);
+    else
+      mfctr(nvl,f,dcp);
+  }
 }
 
-void homfctr(vl,g,dcp)
-VL vl;
-P g;
-DCP *dcp;
+void fctr_wrt_v_p(VL vl,P f,V v,DCP *dcp)
 {
-	P s,t,u,f,h,x;
-	Q e;
-	int d,d0;
-	DCP dc,dct;
+  VL nvl;
+  DCP dc;
 
-	substp(vl,g,vl->v,(P)ONE,&t); fctrp(vl,t,&dc); MKV(vl->v,x);
-	for ( dct = dc; dct; dct = NEXT(dct) )
-		if ( !NUM(dc->c) ) {
-			for ( s = 0, f = dc->c, d = d0 = homdeg(f); f; d = homdeg(f) ) {
-				exthp(vl,f,d,&h); STOQ(d0-d,e); pwrp(vl,x,e,&t);
-				mulp(vl,t,h,&u); addp(vl,s,u,&t); s = t;
-				subp(vl,f,h,&u); f = u;
-			}
-			dc->c = s;
-		}
-	*dcp = dc;
+  if ( !f || NUM(f) ) {
+    NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
+    NEXT(dc) = 0; *dcp = dc;
+    return;
+  } else if ( !qpcheck((Obj)f) )
+    error("fctrp : invalid argument");  
+  else {
+    clctv(vl,f,&nvl);
+    if ( !NEXT(nvl) ) 
+      ufctr(f,1,dcp);
+    else
+      mfctr_wrt_v(nvl,f,v,dcp);
+  }
 }
 
-void mfctr(vl,f,dcp)
-VL vl;
-P f;
-DCP *dcp;
+void homfctr(VL vl,P g,DCP *dcp)
 {
-	DCP dc,dc0,dct,dcs,dcr;
-	P p,pmin,ppmin,cmin,t;
-	VL mvl;
-	Q c;
+  P s,t,u,f,h,x;
+  Q e;
+  int d,d0;
+  DCP dc,dct;
 
-	ptozp(f,1,&c,&p);
-	NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
-	msqfr(vl,p,&dct);
-	for ( ; dct; dct = NEXT(dct) ) {
-		mindegp(vl,COEF(dct),&mvl,&pmin); 
+  substp(vl,g,vl->v,(P)ONE,&t); fctrp(vl,t,&dc); MKV(vl->v,x);
+  for ( dct = dc; dct; dct = NEXT(dct) )
+    if ( !NUM(dc->c) ) {
+      for ( s = 0, f = dc->c, d = d0 = homdeg(f); f; d = homdeg(f) ) {
+        exthp(vl,f,d,&h); STOQ(d0-d,e); pwrp(vl,x,e,&t);
+        mulp(vl,t,h,&u); addp(vl,s,u,&t); s = t;
+        subp(vl,f,h,&u); f = u;
+      }
+      dc->c = s;
+    }
+  *dcp = dc;
+}
+
+void mfctr(VL vl,P f,DCP *dcp)
+{
+  DCP dc,dc0,dct,dcs,dcr;
+  P p,pmin,ppmin,cmin,t;
+  VL mvl;
+  Q c;
+
+  ptozp(f,1,&c,&p);
+  NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
+  msqfr(vl,p,&dct);
+  for ( ; dct; dct = NEXT(dct) ) {
+    mindegp(vl,COEF(dct),&mvl,&pmin); 
 #if 0
-		minlcdegp(vl,COEF(dct),&mvl,&pmin); 
-		min_common_vars_in_coefp(vl,COEF(dct),&mvl,&pmin); 
+    minlcdegp(vl,COEF(dct),&mvl,&pmin); 
+    min_common_vars_in_coefp(vl,COEF(dct),&mvl,&pmin); 
 #endif
-		pcp(mvl,pmin,&ppmin,&cmin);
-		if ( !NUM(cmin) ) {
-			mfctrmain(mvl,cmin,&dcs);
-			for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
-				DEG(dcr) = DEG(dct);
-				reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
-			}
-			for ( ; NEXT(dc); dc = NEXT(dc) );
-			NEXT(dc) = dcs;
-		}
-		mfctrmain(mvl,ppmin,&dcs);
-		for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
-			DEG(dcr) = DEG(dct);
-			reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
-		}
-		for ( ; NEXT(dc); dc = NEXT(dc) );
-		NEXT(dc) = dcs;
-	}
-	adjsgn(f,dc0); *dcp = dc0;
+    pcp(mvl,pmin,&ppmin,&cmin);
+    if ( !NUM(cmin) ) {
+      mfctr(mvl,cmin,&dcs);
+      for ( dcr = NEXT(dcs); dcr; dcr = NEXT(dcr) ) {
+        DEG(dcr) = DEG(dct);
+        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
+      }
+      for ( ; NEXT(dc); dc = NEXT(dc) );
+      NEXT(dc) = NEXT(dcs);
+    }
+    mfctrmain(mvl,ppmin,&dcs);
+    for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
+      DEG(dcr) = DEG(dct);
+      reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
+    }
+    for ( ; NEXT(dc); dc = NEXT(dc) );
+    NEXT(dc) = dcs;
+  }
+  adjsgn(f,dc0); *dcp = dc0;
 }
 
-void adjsgn(p,dc)
-P p;
-DCP dc;
+void mfctr_wrt_v(VL vl,P f,V v,DCP *dcp)
 {
-	int sgn;
-	DCP dct;
-	P c;
+  DCP dc,dc0,dct,dcs,dcr;
+  P p,pmin,ppmin,cmin,t;
+  VL nvl,mvl;
+  Q c;
 
-	for ( dct = dc, sgn = headsgn(p); dct; dct = NEXT(dct) )
-		if ( !EVENN(NM(dct->d)) )
-			sgn *= headsgn(COEF(dct));
-	if ( sgn < 0 ) {
-		chsgnp(COEF(dc),&c); COEF(dc) = c;
-	}
+  ptozp(f,1,&c,&p);
+  NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
+  msqfr(vl,p,&dct);
+  for ( ; dct; dct = NEXT(dct) ) {
+    clctv(vl,COEF(dct),&nvl);
+    reordvar(nvl,v,&mvl);
+    reorderp(mvl,vl,COEF(dct),&pmin);
+    pcp(mvl,pmin,&ppmin,&cmin);
+    if ( !NUM(cmin) ) {
+      mfctrmain(mvl,cmin,&dcs);
+      for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
+        DEG(dcr) = DEG(dct);
+        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
+      }
+      for ( ; NEXT(dc); dc = NEXT(dc) );
+      NEXT(dc) = dcs;
+    }
+    mfctrmain(mvl,ppmin,&dcs);
+    for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
+      DEG(dcr) = DEG(dct);
+      reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
+    }
+    for ( ; NEXT(dc); dc = NEXT(dc) );
+    NEXT(dc) = dcs;
+  }
+  adjsgn(f,dc0); *dcp = dc0;
 }
 
-int headsgn(p)
-P p;
+#if 0
+void adjsgn(P p,DCP dc)
 {
-	if ( !p )
-		return 0;
-	else {
-		while ( !NUM(p) )
-			p = COEF(DC(p));
-		return SGN((Q)p);
-	}
+  int sgn;
+  DCP dct;
+  P c;
+
+  for ( dct = dc, sgn = headsgn(p); dct; dct = NEXT(dct) )
+    if ( !EVENN(NM(dct->d)) )
+      sgn *= headsgn(COEF(dct));
+  if ( sgn < 0 ) {
+    chsgnp(COEF(dc),&c); COEF(dc) = c;
+  }
 }
+#else
+void adjsgn(P p,DCP dc)
+{     
+  DCP dct;
+  P c;
 
-void fctrwithmvp(vl,f,v,dcp)
-VL vl;
-P f;
-V v;
-DCP *dcp;
+  if ( headsgn(COEF(dc)) != headsgn(p) ) {
+    chsgnp(COEF(dc),&c); COEF(dc) = c;
+  }
+  for ( dct = NEXT(dc); dct; dct = NEXT(dct) )
+    if ( headsgn(COEF(dct)) < 0 ) {
+      chsgnp(COEF(dct),&c); COEF(dct) = c;
+    }
+}
+#endif
+
+int headsgn(P p)
 {
-	VL nvl;
-	DCP dc;
+  if ( !p )
+    return 0;
+  else {
+    while ( !NUM(p) )
+      p = COEF(DC(p));
+    return SGN((Q)p);
+  }
+}
 
-	if ( NUM(f) ) {
-		NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
-		NEXT(dc) = 0; *dcp = dc;
-		return;
-	}
+void fctrwithmvp(VL vl,P f,V v,DCP *dcp)
+{
+  VL nvl;
+  DCP dc;
 
-	clctv(vl,f,&nvl);
-	if ( !NEXT(nvl) ) 
-		ufctr(f,1,dcp);
-	else
-		mfctrwithmv(nvl,f,v,dcp);
+  if ( NUM(f) ) {
+    NEWDC(dc); COEF(dc) = f; DEG(dc) = ONE;
+    NEXT(dc) = 0; *dcp = dc;
+    return;
+  }
+
+  clctv(vl,f,&nvl);
+  if ( !NEXT(nvl) ) 
+    ufctr(f,1,dcp);
+  else
+    mfctrwithmv(nvl,f,v,dcp);
 }
 
-void mfctrwithmv(vl,f,v,dcp)
-VL vl;
-P f;
-V v;
-DCP *dcp;
+void mfctrwithmv(VL vl,P f,V v,DCP *dcp)
 {
-	DCP dc,dc0,dct,dcs,dcr;
-	P p,pmin,ppmin,cmin,t;
-	VL mvl;
-	Q c;
+  DCP dc,dc0,dct,dcs,dcr;
+  P p,pmin,ppmin,cmin,t;
+  VL mvl;
+  Q c;
 
-	ptozp(f,1,&c,&p);
-	NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
-	msqfr(vl,p,&dct);
-	for ( ; dct; dct = NEXT(dct) ) {
-		if ( !v )
-			mindegp(vl,COEF(dct),&mvl,&pmin);
-		else {
-			reordvar(vl,v,&mvl); reorderp(mvl,vl,COEF(dct),&pmin);
-		}
-		pcp(mvl,pmin,&ppmin,&cmin);
-		if ( !NUM(cmin) ) {
-			mfctrmain(mvl,cmin,&dcs);
-			for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
-				DEG(dcr) = DEG(dct);
-				reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
-			}
-			for ( ; NEXT(dc); dc = NEXT(dc) );
-			NEXT(dc) = dcs;
-		}
-		mfctrmain(mvl,ppmin,&dcs);
-		for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
-			DEG(dcr) = DEG(dct);
-			reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
-		}
-		for ( ; NEXT(dc); dc = NEXT(dc) );
-		NEXT(dc) = dcs;
-	}
-	*dcp = dc0;
+  ptozp(f,1,&c,&p);
+  NEWDC(dc0); dc = dc0; COEF(dc) = (P)c; DEG(dc) = ONE; NEXT(dc) = 0;
+  msqfr(vl,p,&dct);
+  for ( ; dct; dct = NEXT(dct) ) {
+    if ( !v )
+      mindegp(vl,COEF(dct),&mvl,&pmin);
+    else {
+      reordvar(vl,v,&mvl); reorderp(mvl,vl,COEF(dct),&pmin);
+    }
+    pcp(mvl,pmin,&ppmin,&cmin);
+    if ( !NUM(cmin) ) {
+      mfctrmain(mvl,cmin,&dcs);
+      for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
+        DEG(dcr) = DEG(dct);
+        reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
+      }
+      for ( ; NEXT(dc); dc = NEXT(dc) );
+      NEXT(dc) = dcs;
+    }
+    mfctrmain(mvl,ppmin,&dcs);
+    for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) {
+      DEG(dcr) = DEG(dct);
+      reorderp(vl,mvl,COEF(dcr),&t); COEF(dcr) = t;
+    }
+    for ( ; NEXT(dc); dc = NEXT(dc) );
+    NEXT(dc) = dcs;
+  }
+  *dcp = dc0;
 }
 
-void ufctr(f,hint,dcp)
-P f;
-int hint;
-DCP *dcp;
+void ufctr(P f,int hint,DCP *dcp)
 {
-	P p,c;
-	DCP dc,dct,dcs,dcr;
+  P p,c;
+  DCP dc,dct,dcs,dcr;
 
-	ptozp(f,SGN((Q)UCOEF(f)),(Q *)&c,&p);
-	NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = 0;
-	usqp(p,&dct);
-	for ( *dcp = dc; dct; dct = NEXT(dct) ) {
-		ufctrmain(COEF(dct),hint,&dcs);
-		for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) 
-			DEG(dcr) = DEG(dct);
-		for ( ; NEXT(dc); dc = NEXT(dc) );
-		NEXT(dc) = dcs;
-	}
+  ptozp(f,SGN((Q)UCOEF(f)),(Q *)&c,&p);
+  NEWDC(dc); COEF(dc) = c; DEG(dc) = ONE; NEXT(dc) = 0;
+  usqp(p,&dct);
+  for ( *dcp = dc; dct; dct = NEXT(dct) ) {
+    ufctrmain(COEF(dct),hint,&dcs);
+    for ( dcr = dcs; dcr; dcr = NEXT(dcr) ) 
+      DEG(dcr) = DEG(dct);
+    for ( ; NEXT(dc); dc = NEXT(dc) );
+    NEXT(dc) = dcs;
+  }
 }
 
-void mfctrmain(vl,p,dcp)
-VL vl;
-P p;
-DCP *dcp;
+void mfctrmain(VL vl,P p,DCP *dcp)
 {
-	int i,j,k,*win,np;
-	VL nvl,tvl;
-	VN vn,vnt,vn1;
-	P p0,f,g,f0,g0,s,t,lp,m;
-	P *fp0,*fpt,*l,*tl;
-	DCP dc,dc0,dcl;
-	int count,nv;
+  int i,j,k,*win,np,x;
+  VL nvl,tvl;
+  VN vn,vnt,vn1;
+  P p0,f,g,f0,g0,s,t,lp,m;
+  P *fp0,*fpt,*l,*tl;
+  DCP dc,dc0,dcl;
+  int count,nv;
+  int *nonzero;
 
-	if ( !cmpq(DEG(DC(p)),ONE) ) {
-		NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
-		*dcp = dc; return;
-	}
-	lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);
-	for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
-	W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);
-	W_CALLOC(nv,struct oVN,vn1);
-	for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ ) 
-		vn1[i].v = vn[i].v = tvl->v;
-	vn1[i].v = vn[i].v = 0;
-	count = 0;
-	while  ( 1 ) {
-		while ( 1 ) {
-			count++;
-			for ( i = 0, j = 0; vn[i].v; i++ ) 
-				if ( vn[i].n )
-					vnt[j++].v = (V)i;
-			vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);
-			for ( i = 0; vn1[i].v; i++ )
-				if ( vn1[i].n ) {
-					if ( vn1[i].n > 0 )
-						vn1[i].n = sprime[vn1[i].n];
-					else
-						vn1[i].n = -sprime[-vn1[i].n];
-				}
-			if ( valideval(nvl,dcl,vn1) ) {
-				substvp(nvl,p,vn1,&p0);
-				if ( sqfrchk(p0) ) {
-					ufctr(p0,1,&dc0); 
-					if ( NEXT(NEXT(dc0)) == 0 ) {
-						NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
-						*dcp = dc;
-						return;
-					} else
-						goto MAIN;
-				} 
-			}
-			if ( nextbin(vnt,j) ) 
-				break;
-		}
-		next(vn);
-	}
+  if ( !cmpq(DEG(DC(p)),ONE) ) {
+    NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+    *dcp = dc; return;
+  }
+  lp = LC(p); fctrp(vl,lp,&dcl); clctv(vl,p,&nvl);
+  for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
+  W_CALLOC(nv,struct oVN,vn); W_CALLOC(nv,struct oVN,vnt);
+  W_CALLOC(nv,struct oVN,vn1);
+  W_CALLOC(nv,int,nonzero);
+
+  for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ ) 
+    vn1[i].v = vn[i].v = tvl->v;
+  vn1[i].v = vn[i].v = 0;
+
+  /* determine a heuristic bound of deg(GCD(p,p')) */
+  while ( 1 ) {
+    for ( i = 0; vn1[i].v; i++ )
+      vn1[i].n = ((unsigned int)random())%256+1;
+    substvp(nvl,LC(p),vn1,&p0);
+    if ( p0 ) {
+      substvp(nvl,p,vn1,&p0);
+      if ( sqfrchk(p0) ) {
+        ufctr(p0,1,&dc0); 
+        if ( NEXT(NEXT(dc0)) == 0 ) {
+          NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+          *dcp = dc;
+          return;
+        } else {
+          for ( dc0 = NEXT(dc0), np = 0; dc0; dc0 = NEXT(dc0), np++ );
+          break;
+        }
+      }
+    }
+  }
+  
+  /* determine the position of variables which is not allowed to
+     be set to 0 */
+  for ( i = 0; vn1[i].v; i++ ) {
+    x = vn1[i].n; vn1[i].n = 0;
+    substvp(nvl,LC(p),vn1,&p0);
+    if ( !p0 )
+      vn1[i].n = x;
+    else {
+      substvp(nvl,p,vn1,&p0);
+      if ( !sqfrchk(p0) )
+        vn1[i].n = x;
+      else {
+        ufctr(p0,1,&dc0); 
+        for ( dc0 = NEXT(dc0), j = 0; dc0; dc0 = NEXT(dc0), j++ );
+        if ( j > np )
+          vn1[i].n = x;
+      }
+    }
+  }
+  for ( i = 0; vn1[i].v; i++ )
+    if (vn1[i].n )
+      nonzero[i] = 1;
+
+  count = 0;
+  while  ( 1 ) {
+    while ( 1 ) {
+      count++;
+      for ( i = 0, j = 0; vn[i].v; i++ ) 
+        if ( vn[i].n )
+          vnt[j++].v = (V)i;
+      vnt[j].n = 0; mulsgn(vn,vnt,j,vn1);
+      for ( i = 0; vn1[i].v; i++ )
+        if ( vn1[i].n ) {
+          if ( vn1[i].n > 0 )
+            vn1[i].n = sprime[vn1[i].n];
+          else
+            vn1[i].n = -sprime[-vn1[i].n];
+        }
+      if ( valideval(nvl,dcl,vn1) ) {
+        substvp(nvl,p,vn1,&p0);
+        if ( sqfrchk(p0) ) {
+          ufctr(p0,1,&dc0); 
+          if ( NEXT(NEXT(dc0)) == 0 ) {
+            NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+            *dcp = dc;
+            return;
+          } else {
+            for ( dc = NEXT(dc0), i = 0; dc; dc = NEXT(dc), i++ );
+            if ( i <= np )
+              goto MAIN;
+            if ( i < np )
+              np = i;
+          }
+        }
+      }
+      if ( nextbin(vnt,j) ) 
+        break;
+    }
+    while ( 1 ) {
+      next(vn);
+      for ( i = 0; vn[i].v; i++ )
+        if ( nonzero[i] && !vn[i].n )
+          break;
+      if ( !vn[i].v )
+        break;
+    }
+  }
 MAIN :
 #if 0
-	for ( i = 0; vn1[i].v; i++ )
-		fprintf(stderr,"%d ",vn1[i].n);
-	fprintf(stderr,"\n");
+  for ( i = 0; vn1[i].v; i++ )
+    fprintf(stderr,"%d ",vn1[i].n);
+  fprintf(stderr,"\n");
 #endif
-	for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );
-	fp0 = (P *)ALLOCA((np + 1)*sizeof(P)); 
-	fpt = (P *)ALLOCA((np + 1)*sizeof(P));
-	l = tl = (P *)ALLOCA((np + 1)*sizeof(P));
-	for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) ) 
-		fp0[i-1] = COEF(dc);
+  for ( np = 0, dc = NEXT(dc0); dc; dc = NEXT(dc), np++ );
+  fp0 = (P *)ALLOCA((np + 1)*sizeof(P)); 
+  fpt = (P *)ALLOCA((np + 1)*sizeof(P));
+  l = tl = (P *)ALLOCA((np + 1)*sizeof(P));
+  for ( i = 1, dc = NEXT(dc0); dc; i++, dc = NEXT(dc) ) 
+    fp0[i-1] = COEF(dc);
 #if 0
-	sort_by_deg(np,fp0,fpt);
-	sort_by_deg_rev(np-1,fpt+1,fp0+1);
+  sort_by_deg(np,fp0,fpt);
+  sort_by_deg_rev(np-1,fpt+1,fp0+1);
 #endif
 #if 0
-	sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];
-	sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;
+  sort_by_deg_rev(np,fp0,fpt); fp0[0] = fpt[0];
+  sort_by_deg(np-1,fpt+1,fp0+1); fp0[np] = 0;
 #endif
-	fp0[np] = 0;
-	win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);
-	for ( k = 1, win[0] = 1, --np; ; ) {
-		P h0,lcg,lch;
-		Q c;
+  fp0[np] = 0;
+  win = W_ALLOC(np + 1); f = p; divsp(vl,p0,COEF(dc0),&f0);
+  for ( k = 1, win[0] = 1, --np; ; ) {
+    P h0,lcg,lch;
+    Q c;
 
-		g0 = fp0[win[0]];
-		for ( i = 1; i < k; i++ ) {
-			mulp(vl,g0,fp0[win[i]],&m); g0 = m;
-		}
-		mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lcg);
-		divsp(nvl,f0,g0,&h0);
-		mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lch);
-		mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);		
-		if ( g ) {
-			*tl++ = g; divsp(vl,f,g,&t); 
-			f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0); 
-			for ( i = 0; i < k - 1; i++ ) 
-				for ( j = win[i] + 1; j < win[i + 1]; j++ ) 
-					fp0[j - i - 1] = fp0[j];
-			for ( j = win[k - 1] + 1; j <= np; j++ ) 	
-					fp0[j - k] = fp0[j];
-			if ( ( np -= k ) < k ) break;
-			if ( np - win[0] + 1 < k ) 
-				if ( ++k <= np ) {
-					for ( i = 0; i < k; i++ ) win[i] = i + 1;
-					continue;
-				} else
-					break;
-			else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;
-		} else {
-			if ( ncombi(1,np,k,win) == 0 ) 
-				if ( k == np ) break;
-				else for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;
-		}
-	}
-	*tl++ = f; *tl = 0;
-	for ( dc0 = 0; *l; l++ ) {
-		NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;
-	}
-	NEXT(dc) = 0; *dcp = dc0;
+    g0 = fp0[win[0]];
+    for ( i = 1; i < k; i++ ) {
+      mulp(vl,g0,fp0[win[i]],&m); g0 = m;
+    }
+    mulq((Q)LC(g0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lcg);
+    divsp(nvl,f0,g0,&h0);
+    mulq((Q)LC(h0),(Q)COEF(dc0),&c); estimatelc(nvl,c,dcl,vn1,&lch);
+    mfctrhen2(nvl,vn1,f,f0,g0,h0,lcg,lch,&g);    
+    if ( g ) {
+      *tl++ = g; divsp(vl,f,g,&t); 
+      f = t; divsp(vl,f0,g0,&t); ptozp(t,1,(Q *)&s,&f0); 
+      for ( i = 0; i < k - 1; i++ ) 
+        for ( j = win[i] + 1; j < win[i + 1]; j++ ) 
+          fp0[j - i - 1] = fp0[j];
+      for ( j = win[k - 1] + 1; j <= np; j++ )   
+          fp0[j - k] = fp0[j];
+      if ( ( np -= k ) < k ) break;
+      if ( np - win[0] + 1 < k ) 
+        if ( ++k <= np ) {
+          for ( i = 0; i < k; i++ ) win[i] = i + 1;
+          continue;
+        } else
+          break;
+      else for ( i = 1; i < k; i++ ) win[i] = win[0] + i;
+    } else {
+      if ( ncombi(1,np,k,win) == 0 ) 
+        if ( k == np ) break;
+        else for ( i = 0, ++k; i < k; i++ ) win[i] = i + 1;
+    }
+  }
+  *tl++ = f; *tl = 0;
+  for ( dc0 = 0; *l; l++ ) {
+    NEXTDC(dc0,dc); DEG(dc) = ONE; COEF(dc) = *l;
+  }
+  NEXT(dc) = 0; *dcp = dc0;
 }
 
-void ufctrmain(p,hint,dcp)
-P p;
-int hint;
-DCP *dcp;
+void ufctrmain(P p,int hint,DCP *dcp)
 {
-	ML list;
-	DCP dc;
+  ML list;
+  DCP dc;
 
-	if ( NUM(p) || (UDEG(p) == 1) ) {
-		NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
-		*dcp = dc;
-	} else if ( iscycm(p) )
-		cycm(VR(p),UDEG(p),dcp);
-	else if ( iscycp(p) )
-		cycp(VR(p),UDEG(p),dcp);
-	else {
-		hensel(5,5,p,&list);
-		if ( list->n == 1 ) {
-			NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
-			*dcp = dc;
-		} else 
-			dtest(p,list,hint,dcp);
-	}
+  if ( NUM(p) || (UDEG(p) == 1) ) {
+    NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+    *dcp = dc;
+  } else if ( iscycm(p) )
+    cycm(VR(p),UDEG(p),dcp);
+  else if ( iscycp(p) )
+    cycp(VR(p),UDEG(p),dcp);
+  else {
+    if ( use_new_hensel )
+      hensel2(5,5,p,&list);
+    else
+      hensel(5,5,p,&list);
+    if ( list->n == 1 ) {
+      NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+      *dcp = dc;
+    } else 
+      dtest(p,list,hint,dcp);
+  }
 }
 
-struct oMF {
-	int m;
-	P f;
-};
-
-void calcphi();
-
-void cycm(v,n,dcp)
-V v;
-register int n;
-DCP *dcp;
+void cycm(V v,int n,DCP *dcp)
 {
-	register int i,j;
-	struct oMF *mfp;
-	DCP dc,dc0;
+  register int i,j;
+  struct oMF *mfp;
+  DCP dc,dc0;
 
-	if ( n == 1 ) {
-		NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;
-	} else {
-		mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
-		bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
-		for ( i = 1, j = 0; i <= n; i++ )
-			if ( !(n%i) ) mfp[j++].m = i;
-		mkssum(v,1,1,-1,&mfp[0].f);
-		for ( i = 1; i < j; i++ )
-			calcphi(v,i,mfp);
-		for ( dc0 = 0, i = 0; i < j; i++ ) {
-			NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
-		}
-	}
-	NEXT(dc) = 0; *dcp = dc0;
+  if ( n == 1 ) {
+    NEWDC(dc); mkssum(v,1,1,-1,&COEF(dc)); DEG(dc) = ONE;
+  } else {
+    mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
+    bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
+    for ( i = 1, j = 0; i <= n; i++ )
+      if ( !(n%i) ) mfp[j++].m = i;
+    mkssum(v,1,1,-1,&mfp[0].f);
+    for ( i = 1; i < j; i++ )
+      calcphi(v,i,mfp);
+    for ( dc0 = 0, i = 0; i < j; i++ ) {
+      NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
+    }
+  }
+  NEXT(dc) = 0; *dcp = dc0;
 }
 
-void cycp(v,n,dcp)
-V v;
-register int n;
-DCP *dcp;
+void cycp(V v,int n,DCP *dcp)
 {
-	register int i,j;
-	int n0;
-	struct oMF *mfp;
-	DCP dc,dc0;
+  register int i,j;
+  int n0;
+  struct oMF *mfp;
+  DCP dc,dc0;
 
-	if ( n == 1 ) {
-		NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;
-	} else {
-		n0 = n; n *= 2;
-		mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
-		bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
-		for ( i = 1, j = 0; i <= n; i++ )
-			if ( !(n%i) ) mfp[j++].m = i;
-		mkssum(v,1,1,-1,&mfp[0].f);
-		for ( i = 1; i < j; i++ )
-			calcphi(v,i,mfp);
-		for ( dc0 = 0, i = 0; i < j; i++ ) 
-			if ( n0 % mfp[i].m ) {
-				NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
-			}
-	}
-	NEXT(dc) = 0; *dcp = dc0;
+  if ( n == 1 ) {
+    NEWDC(dc); mkssum(v,1,1,1,&COEF(dc)); DEG(dc) = ONE;
+  } else {
+    n0 = n; n *= 2;
+    mfp = (struct oMF *) ALLOCA((n+1)*sizeof(struct oMF));
+    bzero((char *)mfp,(int)(sizeof(struct oMF)*(n+1)));
+    for ( i = 1, j = 0; i <= n; i++ )
+      if ( !(n%i) ) mfp[j++].m = i;
+    mkssum(v,1,1,-1,&mfp[0].f);
+    for ( i = 1; i < j; i++ )
+      calcphi(v,i,mfp);
+    for ( dc0 = 0, i = 0; i < j; i++ ) 
+      if ( n0 % mfp[i].m ) {
+        NEXTDC(dc0,dc); COEF(dc) = mfp[i].f; DEG(dc) = ONE;
+      }
+  }
+  NEXT(dc) = 0; *dcp = dc0;
 }
 
-void calcphi(v,n,mfp)
-V v;
-int n;
-register struct oMF *mfp;
+void calcphi(V v,int n,struct oMF *mfp)
 {
-	register int i,m;
-	P t,s,tmp;
+  register int i,m;
+  P t,s,tmp;
 
-	for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )
-		if ( !(m%mfp[i].m) ) {
-			mulp(CO,t,mfp[i].f,&s); t = s;
-		}
-	mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);
-	if ( tmp ) 
-		error("calcphi: cannot happen");
+  for ( t = (P)ONE, i = 0, m = mfp[n].m; i < n; i++ )
+    if ( !(m%mfp[i].m) ) {
+      mulp(CO,t,mfp[i].f,&s); t = s;
+    }
+  mkssum(v,m,1,-1,&s); udivpz(s,t,&mfp[n].f,&tmp);
+  if ( tmp ) 
+    error("calcphi: cannot happen");
 }
 
-void mkssum(v,e,s,sgn,r)
-V v;
-int e,s,sgn;
-P *r;
+void mkssum(V v,int e,int s,int sgn,P *r)
 {
-	register int i,sgnt;
-	DCP dc,dc0;
-	Q q;
+  register int i,sgnt;
+  DCP dc,dc0;
+  Q q;
 
-	for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {
-		if ( !dc0 ) {
-			NEWDC(dc0); dc = dc0;
-		} else {
-			NEWDC(NEXT(dc)); dc = NEXT(dc);
-		}
-		STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;
-	}
-	NEXT(dc) = 0; MKP(v,dc0,*r);
+  for ( dc0 = 0, i = s, sgnt = 1; i >= 0; i--, sgnt *= sgn ) {
+    if ( !dc0 ) {
+      NEWDC(dc0); dc = dc0;
+    } else {
+      NEWDC(NEXT(dc)); dc = NEXT(dc);
+    }
+    STOQ(i*e,DEG(dc)); STOQ(sgnt,q),COEF(dc) = (P)q;
+  }
+  NEXT(dc) = 0; MKP(v,dc0,*r);
 }
 
-int iscycp(f)
-P f;
+int iscycp(P f)
 {
-	DCP dc;
-	dc = DC(f);
+  DCP dc;
+  dc = DC(f);
 
-	if ( !UNIQ((Q)COEF(dc)) ) 
-		return ( 0 );
-	dc = NEXT(dc);
-	if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) ) 
-		return ( 0 );
-	return ( 1 );	
+  if ( !UNIQ((Q)COEF(dc)) ) 
+    return ( 0 );
+  dc = NEXT(dc);
+  if ( NEXT(dc) || DEG(dc) || !UNIQ((Q)COEF(dc)) ) 
+    return ( 0 );
+  return ( 1 );  
 }
 
-int iscycm(f)
-P f;
+int iscycm(P f)
 {
-	DCP dc;
+  DCP dc;
 
-	dc = DC(f);
-	if ( !UNIQ((Q)COEF(dc)) ) 
-		return ( 0 );
-	dc = NEXT(dc);
-	if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) ) 
-		return ( 0 );
-	return ( 1 );	
+  dc = DC(f);
+  if ( !UNIQ((Q)COEF(dc)) ) 
+    return ( 0 );
+  dc = NEXT(dc);
+  if ( NEXT(dc) || DEG(dc) || !MUNIQ((Q)COEF(dc)) ) 
+    return ( 0 );
+  return ( 1 );  
 }
 
-void sortfs(dcp)
-DCP *dcp;
+void sortfs(DCP *dcp)
 {
-	int i,k,n,k0,d;
-	DCP dc,dct,t;
-	DCP *a;
+  int i,k,n,k0,d;
+  DCP dc,dct,t;
+  DCP *a;
 
-	dc = *dcp;
-	for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
-	a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
-	for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ ) 
-		a[i] = dct;
-	a[n] = 0;
+  dc = *dcp;
+  for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
+  a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
+  for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ ) 
+    a[i] = dct;
+  a[n] = 0;
 
-	for ( i = 0; i < n; i++ ) {
-		for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ ) 
-			if ( (int)UDEG(COEF(a[k])) < d ) {
-				k0 = k;
-				d = UDEG(COEF(a[k]));
-			}
-		if ( i != k0 ) {
-			t = a[i]; a[i] = a[k0]; a[k0] = t;
-		}
-	}
-	for ( *dcp = dct = a[0], i = 1; i < n; i++ )
-		dct = NEXT(dct) = a[i];
-	NEXT(dct) = 0;
+  for ( i = 0; i < n; i++ ) {
+    for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ ) 
+      if ( (int)UDEG(COEF(a[k])) < d ) {
+        k0 = k;
+        d = UDEG(COEF(a[k]));
+      }
+    if ( i != k0 ) {
+      t = a[i]; a[i] = a[k0]; a[k0] = t;
+    }
+  }
+  for ( *dcp = dct = a[0], i = 1; i < n; i++ )
+    dct = NEXT(dct) = a[i];
+  NEXT(dct) = 0;
 }
 
-void sortfsrev(dcp)
-DCP *dcp;
+void sortfsrev(DCP *dcp)
 {
-	int i,k,n,k0,d;
-	DCP dc,dct,t;
-	DCP *a;
+  int i,k,n,k0,d;
+  DCP dc,dct,t;
+  DCP *a;
 
-	dc = *dcp;
-	for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
-	a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
-	for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ ) 
-		a[i] = dct;
-	a[n] = 0;
+  dc = *dcp;
+  for ( n = 0, dct = dc; dct; dct = NEXT(dct), n++ );
+  a = (DCP *)ALLOCA((n+1)*(sizeof(DCP)));
+  for ( i = 0, dct = dc; dct; dct = NEXT(dct), i++ ) 
+    a[i] = dct;
+  a[n] = 0;
 
-	for ( i = 0; i < n; i++ ) {
-		for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ ) 
-			if ( (int)UDEG(COEF(a[k])) > d ) {
-				k0 = k;
-				d = UDEG(COEF(a[k]));
-			}
-		if ( i != k0 ) {
-			t = a[i]; a[i] = a[k0]; a[k0] = t;
-		}
-	}
-	for ( *dcp = dct = a[0], i = 1; i < n; i++ )
-		dct = NEXT(dct) = a[i];
-	NEXT(dct) = 0;
+  for ( i = 0; i < n; i++ ) {
+    for ( k0 = k = i, d = UDEG(COEF(a[i])); k < n; k++ ) 
+      if ( (int)UDEG(COEF(a[k])) > d ) {
+        k0 = k;
+        d = UDEG(COEF(a[k]));
+      }
+    if ( i != k0 ) {
+      t = a[i]; a[i] = a[k0]; a[k0] = t;
+    }
+  }
+  for ( *dcp = dct = a[0], i = 1; i < n; i++ )
+    dct = NEXT(dct) = a[i];
+  NEXT(dct) = 0;
 }
 
-void nthrootchk(f,dc,fp,dcp)
-P f;
-struct oDUM *dc;
-ML fp;
-DCP *dcp;
+void nthrootchk(P f,struct oDUM *dc,ML fp,DCP *dcp)
 {
-	register int i,k;
-	int e,n,dr,tmp,t;
-	int *tmpp,**tmppp;
-	int **pp,**wpp;
-	LUM fpa,tpa,f0l;
-	UM wt,wq,ws,dvr,f0,fe;
-	N lc;
-	int lcm;
-	int m,b;
+  register int i,k;
+  int e,n,dr,tmp,t;
+  int *tmpp,**tmppp;
+  int **pp,**wpp;
+  LUM fpa,tpa,f0l;
+  UM wt,wq,ws,dvr,f0,fe;
+  N lc;
+  int lcm;
+  int m,b;
 
-	m = fp->mod; b = fp->bound; fe = *((UM *)fp->c); 
-	e = dc->n; f0 = dc->f; nthrootn(NM((Q)COEF(DC(f))),e,&lc);
-	if ( !lc ) {
-		*dcp = 0;
-		return;
-	}
-	lcm = rem(lc,m); W_LUMALLOC(DEG(f0),b,f0l);
-	for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);
-		i >= 0; i-- ) 
-		*(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);
-	dtestroot(m,1,f,f0l,dc,dcp);
-	if ( *dcp ) 
-		return;
-	n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);
-	ptolum(m,b,f,fpa);
-	dvr = W_UMALLOC(n); wq = W_UMALLOC(n); 
-	wt = W_UMALLOC(n); ws = W_UMALLOC(n);
-	cpyum(fe,dvr); divum(m,dvr,f0,wq);
-	t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);
-	for ( k = 0; k <= DEG(wq); k++ ) 
-		COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);
-	for ( i = 1; i < b; i++ ) {
-		pwrlum(m,i+1,f0l,e,tpa);
-		for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa); 
-			k <= n; k++ )
-			COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
-		degum(wt,n); dr = divum(m,wt,dvr,ws);
-		if ( dr >= 0 ) {
-			*dcp = 0;
-			return;
-		}
-		for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )
-			pp[k][i] = COEF(ws)[k];
-		dtestroot(m,i+1,f,f0l,dc,dcp);
-		if ( *dcp )
-			return;
-	}
+  m = fp->mod; b = fp->bound; fe = *((UM *)fp->c); 
+  e = dc->n; f0 = dc->f; nthrootn(NM((Q)COEF(DC(f))),e,&lc);
+  if ( !lc ) {
+    *dcp = 0;
+    return;
+  }
+  lcm = rem(lc,m); W_LUMALLOC(DEG(f0),b,f0l);
+  for ( i = DEG(f0), tmppp = COEF(f0l), tmpp = COEF(f0);
+    i >= 0; i-- ) 
+    *(tmppp[i]) = dmb(m,tmpp[i],lcm,&tmp);
+  dtestroot(m,1,f,f0l,dc,dcp);
+  if ( *dcp ) 
+    return;
+  n = UDEG(f); W_LUMALLOC(n,b,fpa); W_LUMALLOC(n,b,tpa);
+  ptolum(m,b,f,fpa);
+  dvr = W_UMALLOC(n); wq = W_UMALLOC(n); 
+  wt = W_UMALLOC(n); ws = W_UMALLOC(n);
+  cpyum(fe,dvr); divum(m,dvr,f0,wq);
+  t = dmb(m,pwrm(m,lcm,e-1),e,&tmp); DEG(dvr) = DEG(wq);
+  for ( k = 0; k <= DEG(wq); k++ ) 
+    COEF(dvr)[k] = dmb(m,COEF(wq)[k],t,&tmp);
+  for ( i = 1; i < b; i++ ) {
+    pwrlum(m,i+1,f0l,e,tpa);
+    for ( k = 0, pp = COEF(fpa), wpp = COEF(tpa); 
+      k <= n; k++ )
+      COEF(wt)[k] = (pp[k][i]-wpp[k][i]+m)%m;
+    degum(wt,n); dr = divum(m,wt,dvr,ws);
+    if ( dr >= 0 ) {
+      *dcp = 0;
+      return;
+    }
+    for ( k = 0, pp = COEF(f0l); k <= DEG(ws); k++ )
+      pp[k][i] = COEF(ws)[k];
+    dtestroot(m,i+1,f,f0l,dc,dcp);
+    if ( *dcp )
+      return;
+  }
 }
 
-void sqfrp(vl,f,dcp)
-VL vl;
-P f;
-DCP *dcp;
+void sqfrp(VL vl,P f,DCP *dcp)
 {
-	P c,p;
-	DCP dc,dc0;
+  P c,p;
+  DCP dc,dc0;
 
-	if ( !f || NUM(f) ) {
-		NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
-		COEF(dc0) = f; NEXT(dc0) = 0;
-	} else if ( !qpcheck((Obj)f) )
-		error("sqfrp : invalid argument");	
-	else {
-		NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
-		ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);
-		COEF(dc0) = c; NEXT(dc0) = dc;
-		adjsgn(f,dc0);
-	}
+  if ( !f || NUM(f) ) {
+    NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
+    COEF(dc0) = f; NEXT(dc0) = 0;
+  } else if ( !qpcheck((Obj)f) )
+    error("sqfrp : invalid argument");  
+  else {
+    NEWDC(dc0); *dcp = dc0; DEG(dc0) = ONE;
+    ptozp(f,1,(Q *)&c,&p); msqfr(vl,p,&dc);
+    COEF(dc0) = c; NEXT(dc0) = dc;
+    adjsgn(f,dc0);
+  }
 }
 
 /* 
  * f : must be a poly with int coef, ignore int content
  */
-void msqfr(vl,f,dcp)
-VL vl;
-P f;
-DCP *dcp;
+void msqfr(VL vl,P f,DCP *dcp)
 {
-	DCP dc,dct,dcs;
-	P c,p,t,s,pc;
-	VL mvl;
+  DCP dc,dct,dcs;
+  P c,p,t,s,pc;
+  VL mvl;
 
-	ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);
-	if ( NUM(p) ) {
-		*dcp = dc;
-		return;
-	}
-	mindegp(vl,p,&mvl,&s);
+  ptozp(f,1,(Q *)&c,&t); monomialfctr(vl,t,&p,&dc);
+  if ( NUM(p) ) {
+    *dcp = dc;
+    return;
+  }
+  mindegp(vl,p,&mvl,&s);
 #if 0
-	minlcdegp(vl,p,&mvl,&s);
-	min_common_vars_in_coefp(vl,p,&mvl,&s);
+  minlcdegp(vl,p,&mvl,&s);
+  min_common_vars_in_coefp(vl,p,&mvl,&s);
 #endif
-	pcp(mvl,s,&p,&pc);
-	if ( !NUM(pc) ) {
-		msqfr(mvl,pc,&dcs);
-		if ( !dc )
-			dc = dcs;
-		else {
-			for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
-			NEXT(dct) = dcs;
-		}
-	}
-	msqfrmain(mvl,p,&dcs); 
-	for ( dct = dcs; dct; dct = NEXT(dct) ) {
-		reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;
-	}
-	if ( !dc )
-		*dcp = dcs;
-	else {
-		for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
-		NEXT(dct) = dcs; *dcp = dc;
-	}
+  pcp(mvl,s,&p,&pc);
+  if ( !NUM(pc) ) {
+    msqfr(mvl,pc,&dcs);
+    if ( !dc )
+      dc = dcs;
+    else {
+      for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
+      NEXT(dct) = dcs;
+    }
+  }
+  msqfrmain(mvl,p,&dcs); 
+  for ( dct = dcs; dct; dct = NEXT(dct) ) {
+    reorderp(vl,mvl,COEF(dct),&t); COEF(dct) = t;
+  }
+  if ( !dc )
+    *dcp = dcs;
+  else {
+    for ( dct = dc; NEXT(dct); dct = NEXT(dct) );
+    NEXT(dct) = dcs; *dcp = dc;
+  }
 }
 
-void usqp(f,dcp)
-P f;
-DCP *dcp;
+void usqp(P f,DCP *dcp)
 {
-	int index,nindex;
-	P g,c,h;
-	DCP dc;
+  int index,nindex;
+  P g,c,h;
+  DCP dc;
 
-	ptozp(f,1,(Q *)&c,&h); 
-	if ( SGN((Q)LC(h)) < 0 )
-		chsgnp(h,&g);
-	else
-		g = h;
-	for ( index = 0, dc = 0; !dc; index = nindex )
-		hsq(index,5,g,&nindex,&dc);
-	*dcp = dc;
+  ptozp(f,1,(Q *)&c,&h); 
+  if ( SGN((Q)LC(h)) < 0 )
+    chsgnp(h,&g);
+  else
+    g = h;
+  for ( index = 0, dc = 0; !dc; index = nindex )
+    hsq(index,5,g,&nindex,&dc);
+  *dcp = dc;
 }
 
-void msqfrmain(vl,p,dcp)
-VL vl;
-P p;
-DCP *dcp;
+void msqfrmain(VL vl,P p,DCP *dcp)
 {
-	int i,j;
-	VL nvl,tvl;
-	VN vn,vnt,vn1;
-	P gg,tmp,p0,g;
-	DCP dc,dc0,dcr,dcr0;
-	int nv,d,d1;
-	int found;
-	VN svn1;
-	P sp0;
-	DCP sdc0;
+  int i,j;
+  VL nvl,tvl;
+  VN vn,vnt,vn1;
+  P gg,tmp,p0,g;
+  DCP dc,dc0,dcr,dcr0;
+  int nv,d,d1;
+  int found;
+  VN svn1;
+  P sp0;
+  DCP sdc0;
 
-	if ( deg(VR(p),p) == 1 ) {
-		NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
-		*dcp = dc;
-		return;
-	}
-	clctv(vl,p,&nvl);
-	for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
-	if ( nv == 1 ) {
-		usqp(p,dcp);
-		return;
-	}
+  if ( deg(VR(p),p) == 1 ) {
+    NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+    *dcp = dc;
+    return;
+  }
+  clctv(vl,p,&nvl);
+  for ( nv = 0, tvl = nvl; tvl; tvl = NEXT(tvl), nv++);
+  if ( nv == 1 ) {
+    usqp(p,dcp);
+    return;
+  }
 #if 0
-	if ( heusqfrchk(nvl,p) ) {
-		NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
-		*dcp = dc;
-		return;
-	}
+  if ( heusqfrchk(nvl,p) ) {
+    NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+    *dcp = dc;
+    return;
+  }
 #endif
-	W_CALLOC(nv,struct oVN,vn);
-	W_CALLOC(nv,struct oVN,vnt);
-	W_CALLOC(nv,struct oVN,vn1);
-	W_CALLOC(nv,struct oVN,svn1);
-	for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ ) 
-		vn1[i].v = vn[i].v = tvl->v;
-	vn1[i].v = vn[i].v = 0;
+  W_CALLOC(nv,struct oVN,vn);
+  W_CALLOC(nv,struct oVN,vnt);
+  W_CALLOC(nv,struct oVN,vn1);
+  W_CALLOC(nv,struct oVN,svn1);
+  for ( i = 0, tvl = NEXT(nvl); tvl; tvl = NEXT(tvl), i++ ) 
+    vn1[i].v = vn[i].v = tvl->v;
+  vn1[i].v = vn[i].v = 0;
 
-	for  ( dcr0 = 0, g = p, d = deg(VR(g),g), found = 0; ; ) {
-		while ( 1 ) {
-			for ( i = 0, j = 0; vn[i].v; i++ ) 
-				if ( vn[i].n ) vnt[j++].v = (V)i;
-			vnt[j].n = 0;
+  /* determine a heuristic bound of deg(GCD(p,p')) */
+  while ( 1 ) {
+    for ( i = 0; vn1[i].v; i++ )
+      vn1[i].n = ((unsigned int)random())%256+1;
+    substvp(nvl,LC(p),vn1,&tmp);
+    if ( tmp ) {
+      substvp(nvl,p,vn1,&p0);
+      usqp(p0,&dc0);
+      for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
+        if ( DEG(dc) )
+          d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
+      if ( d1 == 0 ) {
+        /* p is squarefree */
+        NEWDC(dc); DEG(dc) = ONE; COEF(dc) = p; NEXT(dc) = 0;
+        *dcp = dc;
+        return;
+      } else {
+        d = d1+1; /* XXX : try searching better evaluation */
+        found = 0;
+        break;
+      }
+    }
+  }
 
-			mulsgn(vn,vnt,j,vn1);
-			substvp(nvl,LC(g),vn1,&tmp);
-			if ( tmp ) {
-				substvp(nvl,g,vn1,&p0);
-				usqp(p0,&dc0);
-				for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
-					if ( DEG(dc) )
-						d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
+  for  ( dcr0 = 0, g = p; ; ) {
+    while ( 1 ) {
+      for ( i = 0, j = 0; vn[i].v; i++ ) 
+        if ( vn[i].n ) vnt[j++].v = (V)i;
+      vnt[j].n = 0;
 
-				if ( d1 == 0 ) {
-					NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;
-					if ( !dcr0 ) 
-						dcr0 = dc;
-					else
-						NEXT(dcr) = dc;
-					*dcp = dcr0;
-					return;
-				}
+      mulsgn(vn,vnt,j,vn1);
+      substvp(nvl,LC(g),vn1,&tmp);
+      if ( tmp ) {
+        substvp(nvl,g,vn1,&p0);
+        usqp(p0,&dc0);
+        for ( d1 = 0, dc = dc0; dc; dc = NEXT(dc) )
+          if ( DEG(dc) )
+            d1 += (QTOS(DEG(dc))-1)*UDEG(COEF(dc));
 
-				if ( d < d1 )
-					goto END;
-				if ( d > d1 ) {
-					d = d1;
-					found = 1;
-					bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));
-					sp0 = p0; sdc0 = dc0;
-					goto END;
-				} 
-				/* d1 == d */
-				if ( found ) {
-					found = 0;
+        if ( d1 == 0 ) {
+          NEWDC(dc); DEG(dc) = ONE; COEF(dc) = g; NEXT(dc) = 0;
+          if ( !dcr0 ) 
+            dcr0 = dc;
+          else
+            NEXT(dcr) = dc;
+          *dcp = dcr0;
+          return;
+        }
+
+        if ( d < d1 )
+          goto END;
+        if ( d > d1 ) {
+          d = d1;
+          found = 1;
+          bcopy((char *)vn1,(char *)svn1,(int)(sizeof(struct oVN)*nv));
+          sp0 = p0; sdc0 = dc0;
+          goto END;
+        } 
+        /* d1 == d */
+        if ( found ) {
+          found = 0;
 #if 0
-				if ( d > d1 ) {
-					d = d1;
-							/*} */
+        if ( d > d1 ) {
+          d = d1;
+              /*} */
 #endif
-					msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;
-					if ( dc ) {
-						if ( !dcr0 )
-							dcr0 = dc;
-						else
-							NEXT(dcr) = dc;
-						for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );
-						if ( NUM(g) ) {
-							NEXT(dcr) = 0; *dcp = dcr0;
-							return;
-						} 
-						d = deg(VR(g),g);
-					}
-				}
-			}
+          msqfrmainmain(nvl,g,svn1,sp0,sdc0,&dc,&gg); g = gg;
+          if ( dc ) {
+            if ( !dcr0 )
+              dcr0 = dc;
+            else
+              NEXT(dcr) = dc;
+            for ( dcr = dc; NEXT(dcr); dcr = NEXT(dcr) );
+            if ( NUM(g) ) {
+              NEXT(dcr) = 0; *dcp = dcr0;
+              return;
+            } 
+            d = deg(VR(g),g);
+          }
+        }
+      }
 END:
-			if ( nextbin(vnt,j) ) 
-				break;
-		}
-		next(vn);
-	}
+      if ( nextbin(vnt,j) ) 
+        break;
+    }
+    next(vn);
+  }
 }
 
-void msqfrmainmain(vl,p,vn,p0,dc0,dcp,pp)
-VL vl;
-P p;
-VN vn;
-P p0;
-DCP dc0;
-DCP *dcp;
-P *pp;
+void msqfrmainmain(VL vl,P p,VN vn,P p0,DCP dc0,DCP *dcp,P *pp)
 {
-	int i,j,k,np;
-	DCP *a;
-	DCP dc,dcr,dcr0,dct;
-	P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;
-	Q q;
-	V v;
+  int i,j,k,np;
+  DCP *a;
+  DCP dc,dcr,dcr0,dct;
+  P g,t,s,u,t0,f,f0,d,d0,g0,h0,x,xx;
+  Q q;
+  V v;
 
-	for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );
-	a = (DCP *)ALLOCA((np + 1)*sizeof(DCP)); 
-	for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) ) 
-		a[np-i-1] = dc;
+  for ( np = 0, dc = dc0; dc; dc = NEXT(dc), np++ );
+  a = (DCP *)ALLOCA((np + 1)*sizeof(DCP)); 
+  for ( i = 0, dc = dc0; dc; i++, dc = NEXT(dc) ) 
+    a[np-i-1] = dc;
 
-	for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p); 
-		i < np; i++ ) { 
-		if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {
-			NEXTDC(dcr0,dcr);
-			DEG(dcr) = DEG(a[i]);
-			COEF(dcr) = f;
-			f = (P)ONE;
-		} else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {
-			diffp(vl,f,v,&s); pcp(vl,s,&t,&u);
-			if ( divtpz(vl,f,t,&s) ) {
-				NEXTDC(dcr0,dcr);
-				DEG(dcr) = DEG(a[i]);
-				COEF(dcr) = s;
-				f = (P)ONE;
-			} else 
-				break;
-		} else { 
-			for ( t = f, t0 = f0,
-				j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {
-				diffp(vl,t,v,&s); t = s;
-				diffp(vl,t0,v,&s); t0 = s;
-			}
-			factorial(k,&q); 
-			divsp(vl,t,(P)q,&s); 
-			for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );
-			if ( DEG(dct) ) {
-				MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx); 
-				divsp(vl,s,xx,&d); 
-			} else {
-				xx = (P)ONE; d = s;
-			}
-			divsp(vl,t0,xx,&t); 
-			divsp(vl,t,(P)q,&s); 
-			ptozp(s,1,(Q *)&t,&d0); 
+  for ( i = 0, dcr0 = 0, f = p, f0 = p0, v = VR(p); 
+    i < np; i++ ) { 
+    if ( (i == (np-1))&&UNIQ(DEG(a[i])) ) {
+      NEXTDC(dcr0,dcr);
+      DEG(dcr) = DEG(a[i]);
+      COEF(dcr) = f;
+      f = (P)ONE;
+    } else if ( (i == (np-1))&&UNIQ(DEG(DC(COEF(a[i])))) ) {
+      diffp(vl,f,v,&s); pcp(vl,s,&t,&u);
+      if ( divtpz(vl,f,t,&s) ) {
+        NEXTDC(dcr0,dcr);
+        DEG(dcr) = DEG(a[i]);
+        COEF(dcr) = s;
+        f = (P)ONE;
+      } else 
+        break;
+    } else { 
+      for ( t = f, t0 = f0,
+        j = 0, k = QTOS(DEG(a[i]))-1; j < k; j++ ) {
+        diffp(vl,t,v,&s); t = s;
+        diffp(vl,t0,v,&s); t0 = s;
+      }
+      factorial(k,&q); 
+      divsp(vl,t,(P)q,&s); 
+      for ( dct = DC(s); NEXT(dct); dct = NEXT(dct) );
+      if ( DEG(dct) ) {
+        MKV(VR(s),x); pwrp(vl,x,DEG(dct),&xx); 
+        divsp(vl,s,xx,&d); 
+      } else {
+        xx = (P)ONE; d = s;
+      }
+      divsp(vl,t0,xx,&t); 
+      divsp(vl,t,(P)q,&s); 
+      ptozp(s,1,(Q *)&t,&d0); 
 
-			for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );
-			if ( DEG(dct) )
-				divsp(vl,COEF(a[i]),xx,&g0);
-			else {
-				xx = (P)ONE; g0 = COEF(a[i]);
-			}
+      for ( dct = DC(COEF(a[i])); NEXT(dct); dct = NEXT(dct) );
+      if ( DEG(dct) )
+        divsp(vl,COEF(a[i]),xx,&g0);
+      else {
+        xx = (P)ONE; g0 = COEF(a[i]);
+      }
 
-			pcp(vl,d,&t,&u); d = t;
-			ptozp(g0,1,(Q *)&u,&t); g0 = t;
+      pcp(vl,d,&t,&u); d = t;
+      ptozp(g0,1,(Q *)&u,&t); g0 = t;
 
-		{
-			DCP dca,dcb;
+    {
+      DCP dca,dcb;
 
-			fctrp(vl,LC(d),&dca);
-			for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {
-				if ( NUM(COEF(dcb)) ) {
-					mulp(vl,u,COEF(dcb),&t); u = t;
-				} else {
-					Q qq;
-					P tt;
+      fctrp(vl,LC(d),&dca);
+      for ( dcb = dca, u = (P)ONE; dcb; dcb = NEXT(dcb) ) {
+        if ( NUM(COEF(dcb)) ) {
+          mulp(vl,u,COEF(dcb),&t); u = t;
+        } else {
+          Q qq;
+          P tt;
 
-					pwrp(vl,COEF(dcb),DEG(a[i]),&s);
-					for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );
-					STOQ(j,qq);
-					if ( cmpq(qq,DEG(dcb)) > 0 )
-						qq = DEG(dcb);
-					pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;
-				}
-			}
-			divsp(vl,d0,g0,&h0);
-		}
-			mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s); 
-			if ( s ) {
-				mulp(vl,s,xx,&g); 
-				pwrp(vl,g,DEG(a[i]),&t);
-				if ( divtpz(vl,f,t,&s) ) {
-					NEXTDC(dcr0,dcr);
-					DEG(dcr) = DEG(a[i]); COEF(dcr) = g;
-					f = s; substvp(vl,f,vn,&f0);
-				} else
-					break;
-			} else 
-				break;
-		}
-	}
-	*pp = f;	
-	if ( dcr0 ) 
-		NEXT(dcr) = 0;
-	*dcp = dcr0;
+          pwrp(vl,COEF(dcb),DEG(a[i]),&s);
+          for ( t = LC(f), j = 0; divtpz(vl,t,s,&tt); j++, t = tt );
+          STOQ(j,qq);
+          if ( cmpq(qq,DEG(dcb)) > 0 )
+            qq = DEG(dcb);
+          pwrp(vl,COEF(dcb),qq,&t); mulp(vl,u,t,&s); u = s;
+        }
+      }
+      divsp(vl,d0,g0,&h0);
+    }
+      mfctrhen2(vl,vn,d,d0,g0,h0,u,LC(d),&s); 
+      if ( s ) {
+        mulp(vl,s,xx,&g); 
+        pwrp(vl,g,DEG(a[i]),&t);
+        if ( divtpz(vl,f,t,&s) ) {
+          NEXTDC(dcr0,dcr);
+          DEG(dcr) = DEG(a[i]); COEF(dcr) = g;
+          f = s; substvp(vl,f,vn,&f0);
+        } else
+          break;
+      } else 
+        break;
+    }
+  }
+  *pp = f;  
+  if ( dcr0 ) 
+    NEXT(dcr) = 0;
+  *dcp = dcr0;
 }
 
-void mfctrhen2(vl,vn,f,f0,g0,h0,lcg,lch,gp)
-VL vl;
-VN vn;
-P f;
-P f0,g0,h0,lcg,lch;
-P *gp;
+void mfctrhen2(VL vl,VN vn,P f,P f0,P g0,P h0,P lcg,P lch,P *gp)
 {
-	V v;
-	P f1,lc,lc0,lcg0,lch0;
-	P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;
-	Q bb,qk,s;
-	Q cbd; 
-	int dbd;
-	int d;
+  V v;
+  P f1,lc,lc0,lcg0,lch0;
+  P m,ff0,gg0,hh0,gk,hk,ggg,gggr,hhh,ak,bk,tmp;
+  Q bb,qk,s;
+  Q cbd; 
+  int dbd;
+  int d;
 
-	if ( NUM(g0) ) {
-		*gp = (P)ONE;
-		return;
-	} 
+  if ( NUM(g0) ) {
+    *gp = (P)ONE;
+    return;
+  } 
 
-	v = VR(f); d = deg(v,f);
-	if ( d == deg(v,g0) ) {
-		pcp(vl,f,gp,&tmp); 
-		return;
-	}		
+  v = VR(f); d = deg(v,f);
+  if ( d == deg(v,g0) ) {
+    pcp(vl,f,gp,&tmp); 
+    return;
+  }    
 
-	mulp(vl,lcg,lch,&lc);
-	if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {
-		*gp = 0;
-		return;
-	}
-	mulp(vl,(P)s,f,&f1);
-	dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);
+  mulp(vl,lcg,lch,&lc);
+  if ( !divtpz(vl,lc,LC(f),(P *)&s) ) {
+    *gp = 0;
+    return;
+  }
+  mulp(vl,(P)s,f,&f1);
+  dbd = dbound(VR(f1),f1) + 1; cbound(vl,f1,&cbd);
 
-	substvp(vl,lc,vn,&lc0);
-	divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0); 
-	substvp(vl,lcg,vn,&lcg0);
-	divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0); 
-	substvp(vl,lch,vn,&lch0);
-	divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0); 
-	addq(cbd,cbd,&bb);
-	henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;
-	henmv(vl,vn,f1,gk,hk,ak,bk,
-		lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh); 
+  substvp(vl,lc,vn,&lc0);
+  divq((Q)lc0,(Q)LC(f0),(Q *)&m); mulp(vl,f0,m,&ff0); 
+  substvp(vl,lcg,vn,&lcg0);
+  divq((Q)lcg0,(Q)LC(g0),(Q *)&m); mulp(vl,g0,m,&gg0); 
+  substvp(vl,lch,vn,&lch0);
+  divq((Q)lch0,(Q)LC(h0),(Q *)&m); mulp(vl,h0,m,&hh0); 
+  addq(cbd,cbd,&bb);
+  henzq1(gg0,hh0,bb,&bk,&ak,&qk); gk = gg0; hk = hh0;
+  henmv(vl,vn,f1,gk,hk,ak,bk,
+    lcg,lch,lcg0,lch0,qk,dbd,&ggg,&hhh); 
 
-	if ( divtpz(vl,f1,ggg,&gggr) )
-		pcp(vl,ggg,gp,&tmp);
-	else
-		*gp = 0;
+  if ( divtpz(vl,f1,ggg,&gggr) )
+    pcp(vl,ggg,gp,&tmp);
+  else
+    *gp = 0;
 }
 
-int sqfrchk(p)
-P p;
+int sqfrchk(P p)
 {
-	Q c;
-	P f;
-	DCP dc;
+  Q c;
+  P f;
+  DCP dc;
 
-	ptozp(p,SGN((Q)UCOEF(p)),&c,&f); usqp(f,&dc); 
-	if ( NEXT(dc) || !UNIQ(DEG(dc)) ) 
-		return ( 0 );
-	else 
-		return ( 1 );
+  ptozp(p,SGN((Q)UCOEF(p)),&c,&f); usqp(f,&dc); 
+  if ( NEXT(dc) || !UNIQ(DEG(dc)) ) 
+    return ( 0 );
+  else 
+    return ( 1 );
 }
 
-int cycchk(p)
-P p;
+int cycchk(P p)
 {
-	Q c;
-	P f;
+  Q c;
+  P f;
 
-	ptozp(p,SGN((Q)UCOEF(p)),&c,&f);
-	if ( iscycp(f) || iscycm(f) )
-		return 0;
-	else 
-		return 1;
+  ptozp(p,SGN((Q)UCOEF(p)),&c,&f);
+  if ( iscycp(f) || iscycm(f) )
+    return 0;
+  else 
+    return 1;
 }
 
-int zerovpchk(vl,p,vn)
-VL vl;
-P p;
-VN vn;
+int zerovpchk(VL vl,P p,VN vn)
 {
-	P t;
+  P t;
 
-	substvp(vl,p,vn,&t);
-	if ( t )
-		return ( 0 );
-	else
-		return ( 1 );
+  substvp(vl,p,vn,&t);
+  if ( t )
+    return ( 0 );
+  else
+    return ( 1 );
 }
 
-int valideval(vl,dc,vn)
-VL vl;
-DCP dc;
-VN vn;
+int valideval(VL vl,DCP dc,VN vn)
 {
-	DCP dct;
-	Q *a;
-	int i,j,n;
-	N q,r;
+  DCP dct;
+  Q *a;
+  int i,j,n;
+  N q,r;
 
-	for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );
-	W_CALLOC(n,Q,a);
-	for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {
-		substvp(vl,COEF(dct),vn,(P *)&a[i]);
-		if ( !a[i] )
-			return ( 0 );
+  for ( dct = NEXT(dc), n = 0; dct; dct = NEXT(dct), n++ );
+  W_CALLOC(n,Q,a);
+  for ( dct = NEXT(dc), i = 0; dct; dct = NEXT(dct), i++ ) {
+    substvp(vl,COEF(dct),vn,(P *)&a[i]);
+    if ( !a[i] )
+      return ( 0 );
 
-		for ( j = 0; j < i; j++ ) {
-			divn(NM(a[j]),NM(a[i]),&q,&r);
-			if ( !r )
-				return ( 0 );
-			divn(NM(a[i]),NM(a[j]),&q,&r);
-			if ( !r )
-				return ( 0 );
-		}
-	}
-	return ( 1 );
+    for ( j = 0; j < i; j++ ) {
+      divn(NM(a[j]),NM(a[i]),&q,&r);
+      if ( !r )
+        return ( 0 );
+      divn(NM(a[i]),NM(a[j]),&q,&r);
+      if ( !r )
+        return ( 0 );
+    }
+  }
+  return ( 1 );
 }
 
-void estimatelc(vl,c,dc,vn,lcp)
-VL vl;
-Q c;
-DCP dc;
-VN vn;
-P *lcp;
+void estimatelc(VL vl,Q c,DCP dc,VN vn,P *lcp)
 {
-	int i;
-	DCP dct;
-	P r,s,t;
-	Q c0,c1,c2;
+  int i;
+  DCP dct;
+  P r,s,t;
+  Q c0,c1,c2;
 
-	for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
-		if ( NUM(COEF(dct)) ) {
-			mulp(vl,r,COEF(dct),&s); r = s;
-		} else {
-			substvp(vl,COEF(dct),vn,(P *)&c0);
-			for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {
-				divq(c1,c0,&c2);
-				if ( !INT(c2) )
-					break;
-				else
-					c1 = c2;
-			}
-			if ( i ) {
-				STOQ(i,c1);
-				pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;
-			}
-		}
-	}
-	*lcp = r;
+  for ( dct = dc, r = (P)ONE; dct; dct = NEXT(dct) ) {
+    if ( NUM(COEF(dct)) ) {
+      mulp(vl,r,COEF(dct),&s); r = s;
+    } else {
+      substvp(vl,COEF(dct),vn,(P *)&c0);
+      for ( i = 0, c1 = c; i < (int)QTOS(DEG(dct)); i++ ) {
+        divq(c1,c0,&c2);
+        if ( !INT(c2) )
+          break;
+        else
+          c1 = c2;
+      }
+      if ( i ) {
+        STOQ(i,c1);
+        pwrp(vl,COEF(dct),c1,&s); mulp(vl,r,s,&t); r = t;
+      }
+    }
+  }
+  *lcp = r;
 }
 
-void monomialfctr(vl,p,pr,dcp)
-VL vl;
-P p;
-P *pr;
-DCP *dcp;
+void monomialfctr(VL vl,P p,P *pr,DCP *dcp)
 {
-	VL nvl,avl;
-	Q d;
-	P f,t,s;
-	DCP dc0,dc;
+  VL nvl,avl;
+  Q d;
+  P f,t,s;
+  DCP dc0,dc;
 
-	clctv(vl,p,&nvl);
-	for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {
-		getmindeg(avl->v,f,&d);
-		if ( d ) {
-			MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;
-			pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;
-		}
-	}
-	if ( dc0 )
-		NEXT(dc) = 0;
-	*pr = f; *dcp = dc0;
+  clctv(vl,p,&nvl);
+  for ( dc0 = 0, avl = nvl, f = p; avl; avl = NEXT(avl) ) {
+    getmindeg(avl->v,f,&d);
+    if ( d ) {
+      MKV(avl->v,t); NEXTDC(dc0,dc); DEG(dc) = d; COEF(dc) = t;
+      pwrp(vl,t,d,&s); divsp(vl,f,s,&t); f = t;
+    }
+  }
+  if ( dc0 )
+    NEXT(dc) = 0;
+  *pr = f; *dcp = dc0;
 }
 
-void afctr(vl,p0,p,dcp)
-VL vl;
-P p,p0;
-DCP *dcp;
+void afctr(VL vl,P p0,P p,DCP *dcp)
 {
-	DCP dc,dc0,dcr,dct,dcs;
-	P t;
-	VL nvl;
+  DCP dc,dc0,dcr,dct,dcs;
+  P t;
+  VL nvl;
 
-	if ( VR(p) == VR(p0) ) {
-		NEWDC(dc);
-		DEG(dc) = ONE;
-		COEF(dc) = p;
-		NEXT(dc) = 0;
-		*dcp = dc;
-		return;
-	}
+  if ( VR(p) == VR(p0) ) {
+    NEWDC(dc);
+    DEG(dc) = ONE;
+    COEF(dc) = p;
+    NEXT(dc) = 0;
+    *dcp = dc;
+    return;
+  }
 
-	clctv(vl,p,&nvl);
-	if ( !NEXT(nvl) ) 
-		ufctr(p,1,&dc);
-	else {
-		sqa(vl,p0,p,&dc);
-		for ( dct = dc; dct; dct = NEXT(dct) ) {
-			pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;
-		}
-	}
-	if ( NUM(COEF(dc)) ) 
-		dcr = NEXT(dc);
-	else
-		dcr = dc;
-	for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {
-		afctrmain(vl,p0,COEF(dcr),1,&dcs);
+  clctv(vl,p,&nvl);
+  if ( !NEXT(nvl) ) 
+    ufctr(p,1,&dc);
+  else {
+    sqa(vl,p0,p,&dc);
+    for ( dct = dc; dct; dct = NEXT(dct) ) {
+      pmonic(vl,p0,COEF(dct),&t); COEF(dct) = t;
+    }
+  }
+  if ( NUM(COEF(dc)) ) 
+    dcr = NEXT(dc);
+  else
+    dcr = dc;
+  for ( dc0 = 0; dcr; dcr = NEXT(dcr) ) {
+    afctrmain(vl,p0,COEF(dcr),1,&dcs);
 
-		for ( dct = dcs; dct; dct = NEXT(dct) ) 
-			DEG(dct) = DEG(dcr);
-		if ( !dc0 )
-			dc0 = dcs;
-		else {
-			for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );
-			NEXT(dct) = dcs;
-		}
-	}
-	*dcp = dc0;
+    for ( dct = dcs; dct; dct = NEXT(dct) ) 
+      DEG(dct) = DEG(dcr);
+    if ( !dc0 )
+      dc0 = dcs;
+    else {
+      for ( dct = dc0; NEXT(dct); dct = NEXT(dct) );
+      NEXT(dct) = dcs;
+    }
+  }
+  *dcp = dc0;
 }
 
-void afctrmain(vl,p0,p,init,dcp)
-VL vl;
-P p,p0;
-int init;
-DCP *dcp;
+void afctrmain(VL vl,P p0,P p,int init,DCP *dcp)
 {
-	P x,y,s,m,a,t,u,pt,pt1,res,g;
-	Q q;
-	DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;
-	V v,v0;
+  P x,y,s,m,a,t,u,pt,pt1,res,g;
+  Q q;
+  DCP dc,dc0,dcsq0,dcr0,dcr,dct,dcs;
+  V v,v0;
 
-	if ( !cmpq(DEG(DC(p)),ONE) ) {
-		NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;
-		pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;
-		return;
-	}
+  if ( !cmpq(DEG(DC(p)),ONE) ) {
+    NEWDC(dc); DEG(dc) = ONE; NEXT(dc) = 0;
+    pmonic(vl,p0,p,&COEF(dc)); *dcp = dc;
+    return;
+  }
 
-	v = VR(p); MKV(v,x);
-	v0 = VR(p0); MKV(v0,y);
-	STOQ(init,q),s = (P)q;
-	mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a); 
-	substp(vl,p,v,t,&pt); 
-	remsdcp(vl,pt,p0,&pt1); 
+  v = VR(p); MKV(v,x);
+  v0 = VR(p0); MKV(v0,y);
+  STOQ(init,q),s = (P)q;
+  mulp(vl,s,y,&m); subp(vl,x,m,&t); addp(vl,x,m,&a); 
+  substp(vl,p,v,t,&pt); 
+  remsdcp(vl,pt,p0,&pt1); 
 
 /*
-	if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) ) 
-		resultp(vl,v0,p0,pt1,&res);
-	else 
-		srcrnorm(vl,v0,pt1,p0,&res);
+  if ( ( deg(v0,p0) <= 3 ) || ( TCPQ(p0) <= 2 ) ) 
+    resultp(vl,v0,p0,pt1,&res);
+  else 
+    srcrnorm(vl,v0,pt1,p0,&res);
 */
 #if 0
-	srcr(vl,v0,pt1,p0,&res);
+  srcr(vl,v0,pt1,p0,&res);
 #endif
-	resultp(vl,v0,p0,pt1,&res); 
-	usqp(res,&dcsq0); 
-	for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {
-		if ( UNIQ(DEG(dct)) ) 
-			ufctr(COEF(dct),deg(v0,p0),&dcs);
-		else
-			ufctr(COEF(dct),1,&dcs);
-		for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
-			DEG(dcr) = DEG(dct);
-		if ( !dc0 ) {
-			dc0 = NEXT(dcs);
-			dc = dc0;	
-		} else {
-			for ( ; NEXT(dc); dc = NEXT(dc) );
-			NEXT(dc) = NEXT(dcs);
-		}
-	}
-	sortfs(&dc0);
+  resultp(vl,v0,p0,pt1,&res); 
+  usqp(res,&dcsq0); 
+  for ( dc0 = 0, dct = dcsq0; dct; dct = NEXT(dct) ) {
+    if ( UNIQ(DEG(dct)) ) 
+      ufctr(COEF(dct),deg(v0,p0),&dcs);
+    else
+      ufctr(COEF(dct),1,&dcs);
+    for ( dcr = dcs; dcr; dcr = NEXT(dcr) )
+      DEG(dcr) = DEG(dct);
+    if ( !dc0 ) {
+      dc0 = NEXT(dcs);
+      dc = dc0;  
+    } else {
+      for ( ; NEXT(dc); dc = NEXT(dc) );
+      NEXT(dc) = NEXT(dcs);
+    }
+  }
+  sortfs(&dc0);
 
-	for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {
-		if ( !UNIQ(DEG(dc)) ) {
-			substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s); 
-			gcda(vl,p0,s,g,&u); 
-			if ( !NUM(u) && (VR(u) == v)) {
-				afctrmain(vl,p0,u,init+1,&dct);
-				for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {
-					mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);
-				}
-				pdiva(vl,p0,g,t,&s); g = s;
-				if ( !dcr0 )
-					dcr0 = dct;
-				else 
-					NEXT(dcr) = dct;
-				for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );
-			}
-		} else {
-			substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s); 
-			gcda(vl,p0,s,g,&u); 
-			if ( !NUM(u) && (VR(u) == v)) {
-				NEXTDC(dcr0,dcr);
-				DEG(dcr) = ONE;
-				COEF(dcr) = u;
-				pdiva(vl,p0,g,u,&t); g = t;
-			} 
-		}
-	}
-	if ( dcr0 )
-		NEXT(dcr) = 0;
-	*dcp = dcr0;
+  for ( g = p, dcr = dcr0 = 0, dc = dc0; dc; dc = NEXT(dc) ) {
+    if ( !UNIQ(DEG(dc)) ) {
+      substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s); 
+      gcda(vl,p0,s,g,&u); 
+      if ( !NUM(u) && (VR(u) == v)) {
+        afctrmain(vl,p0,u,init+1,&dct);
+        for ( dcs = dct, t = (P)ONE; dcs; dcs = NEXT(dcs) ) {
+          mulp(vl,t,COEF(dcs),&s); remsdcp(vl,s,p0,&t);
+        }
+        pdiva(vl,p0,g,t,&s); g = s;
+        if ( !dcr0 )
+          dcr0 = dct;
+        else 
+          NEXT(dcr) = dct;
+        for ( dcr = dct; NEXT(dcr); dcr = NEXT(dcr) );
+      }
+    } else {
+      substp(vl,COEF(dc),v,a,&pt); remsdcp(vl,pt,p0,&s); 
+      gcda(vl,p0,s,g,&u); 
+      if ( !NUM(u) && (VR(u) == v)) {
+        NEXTDC(dcr0,dcr);
+        DEG(dcr) = ONE;
+        COEF(dcr) = u;
+        pdiva(vl,p0,g,u,&t); g = t;
+      } 
+    }
+  }
+  if ( dcr0 )
+    NEXT(dcr) = 0;
+  *dcp = dcr0;
 }