\\ this is necessary to set the priorities of the variables for gp 
x;  \\ this is the variable for the polynomial to be factored 
y;  \\ this is the variable for polynomials which arise as resultants 
Z;  \\ extensions of Q_p are formed by adjoining a polynomial in Z 

INFINITY = 9999;  /* valuation of 0 */

DEBUG=0; /* higher DEBUG levels print out more information */

/* find extended gcd(u,v), return u1, u2 such that u1*u+u2*v = 1 */
{
xgcd(u,v)=
if(DEBUG>=10,print("in xgcd(",u,","v,")"));
    u1=1;    u2=0;    u3=u;
    v1=0;    v2=1;    v3=v;

    while(v3 != 0,
        l=divrem(u3,v3);
	q = l[1];
	t1=u1-v1*q;
	t2=u2-v2*q;
	t3=u3-v3*q;
	u1=v1;	u2=v2;	u3=v3;
	v1=t1;	v2=t2;	v3=t3;
    );
    u1 = u1/u3;  /* normalize to gcd = 1 */
    u2 = u2/u3;
    return([u1,u2]);
}


\\ find gcd of two p-adic polynomials, using Euclidean algorithm 
\\ this is vulnerable to precision problems, and should eventually be 
\\ replaced by the subresultant method 
{
Zp_gcd(u,v,w,p,prec, t,m,du,dv,i,uu)=
if(DEBUG>=10,print("in Zp_gcd(",u,",",v,")"));
    while(poldegree(v)>0,
      du = poldegree(u);
      dv = poldegree(v);
      if(du < dv,
        t = u;
	u = v;
	v = t,
/* else	*/
	m = polcoeff(u,du)/polcoeff(v,dv);
	u = u - m*v*x^(du-dv);
        uu = 0;
        for(i=0,du,
	    if(padicprec(polcoeff(u,i),p) != valuation(polcoeff(u,i),p),
	        uu += polcoeff(u,i)*x^i));
	u = uu;
    ));
    return(u/polcoeff(u,poldegree(u)));
}

\\ give all coefficients the same precision 
\\ coef is an element of an extension of Q_p, given as a poly mod w(Z)
{
my_truncate(coef,w,p,prec, qq, tt,ww,ii)=
        tt = lift(coef,Z);
	qq = 0;
        for(ii=0,poldegree(w),
	    ww = polcoeff(tt,ii);
	    if(type(ww) == "t_PADIC", ww = truncate(ww));
/*if(DEBUG>=0,print("ii = ",ii,", coef = ",polcoeff(tt,ii),", ww = ",ww));*/
            if(ww != 0,
	    qq = qq + (ww*(1+O(p^prec)))*Z^ii);
	);
	return(Mod(qq,w))
}

/* fix precisions, get rid of (0+O(p^prec)) terms */
{
clean(R,w,p,prec, vv,i,qq)=
        vv=0;
	for(i=0,poldegree(R),
	    qq = my_truncate(polcoeff(R,i),w,p,prec);
	    if(qq != 0, vv = vv+qq*y^i);
	);
	return(vv);
}

{
get_precision(R,w,p,prec, pp,minprec,c)=
	minprec = INFINITY;
	for(i=0,poldegree(R),
	    c = polcoeff(R,i);
	    pp = padicprec(c,p)-valuation(c,p);
	    if((pp < minprec) && (pp > 0), minprec = pp);
	);
	return(minprec);
}

{
myvaluation(c,p)=
        if(c==0,INFINITY,valuation(c,p))
}

{
val(R,lambda,p,prec,  vv,ww)=
if(DEBUG>=10,print("in val, R=",R));
        vv=0;
	for(i=0,poldegree(R),
	    vl = myvaluation(polcoeff(R,i),p);
	    ww = p^(-vl)*lambda^i;
	    if((ww > vv) && (vl < INFINITY), vv=ww);
	);
	vv
}

{
lemma6_10(B,C,u,v,e,lambda,mu,w,p,prec,  U1,V,e_val,l,i)=
        l = divrem(e*v,B);
        Q = l[1];
	V = l[2];
	U1 = (e*u) + (Q*C);
	e_val = val(e,lambda,p,prec);
	for(i=0,poldegree(U1),
/* the 0.99999 is to avoid roundoff errors when lambda is irrational */
	    if(val(polcoeff(U1,i)*y^i,lambda,p,prec) < (e_val*0.99999),
	         U1 = U1 - polcoeff(U1,i)*y^i;
	    );
	);
	U1=clean(U1,w,p,prec);
	return([U1,V]);
}

{
hensel0(R,B0,w,p,prec, C0,l,u,v,r,s,n,a0,val_a0)=
if(DEBUG>=1,print("in hensel0, \nR=",R,"\nB=",B0));
    R = subst(make_padic(R,w,p,prec),x,y);
    n=poldegree(R);
    a0 = polcoeff(R,0);
    val_a0=myvaluation(a0,p);
    r=n/gcd(val_a0,n);
    s=val_a0/gcd(val_a0,n);
    R /= a0;
    B0 = B0/polcoeff(B,0);
    C0 = divrem(R,B0)[1];
    B0 = subst(clean(B0,w,p,prec),x,y);
    C0 = clean(C0,w,p,prec);
    l = xgcd(B0,C0);
    u = l[1];
    v = l[2];
    B0 = hensel(R,B0,C0,u,v,0,r,s,w,p,prec);
    return(hensel(R,B0,C0,u,v,0,r,s,w,p,prec));
}

{
hensel(R,B,C,u,v,h,r,s,w,p,prec, zzz,BB,CC,U2,V2,l,e)=
if(DEBUG>=10,
print("in hensel, R = ",R);
print("u = ",u);
print("v = ",v));
	mu = p^(-h);
	lambda = p^(-s/r);
        while(true,
	    e = R - B*C;
	    e = clean(e,w,p,prec);
	    l = lemma6_10(B,C,u,v,e,lambda,mu,w,p,prec);
	    U2 = l[1];
	    V2 = l[2];
	    B = B + V2;
	    C = C + U2;
	    B = clean(B,w,p,prec);
	    C = clean(C,w,p,prec);

	    if(clean(R-B*C,w,p,prec) == 0, break());
	);
	new_prec = get_precision(R,w,p,prec);
/*new_prec = prec; */
	BB = make_monic(B,w,p,new_prec);
	CC = make_monic(C,w,p,new_prec);
        return(BB); /* Hensel shouldn't do factoring, so quit here */
        BB = subst(BB,y,x);
        CC = subst(CC,y,x);
	zzz = find_factor(BB,w,p,prec);
	zzz = subst(zzz,x,y);
print("zzz = ",zzz);
	return(zzz);
}

{
cor6_34(R,b,g,w,p,prec,  ll,n,a0,r,s,l,m,B0,C0,u,v)=
if(DEBUG>=1,print("in Cor 6.34, R = ",R));
    ll= xgcd(b,g);
    mu=ll[1];
    nu=ll[2];  
    n=poldegree(R);
/*    R = R/polcoeff(R,0);  /* make constant term 1 */
    a0=myvaluation(polcoeff(R,0),p);
    r=n/gcd(a0,n);
    s=a0/gcd(a0,n);
    l=poldegree(b);
    m=poldegree(g);

    B0 = my_lift(b,s,r,w,p,prec);
    C0 = my_lift(g,s,r,w,p,prec);
    u = my_lift(mu,s,r,w,p,prec);
    v = my_lift(nu,s,r,w,p,prec);

    return(hensel(p^(-a0)*R,B0,C0,u,v,0,r,s,w,p,prec));
}

{
/* remove ith element of L */
remove_elt(L,i, j,L2)=
        L2 = [];
	for(j=1,i-1,L2 = concat(L2,[L[j]]));
	for(j=i+1,length(L),L2 = concat(L2,[L[j]]));
	L2
}

{
newtondiag(R,p,prec,  d,L,x1,y1,x2,y2,x3,y3,i,j,l)=
if(DEBUG==1,print("in newtondiag, R = ",R));
	L=[];
	for(i=0,poldegree(R),
		if(polcoeff(R,i) != 0,
			L = concat(L,[[i,myvaluation(polcoeff(R,i),p)]]);
		);
	);
	l = length(L);
	for(j=2,l-1,
		x1 = L[j-1][1];
		y1 = L[j-1][2];
		x2 = L[j][1];
		y2 = L[j][2];
		x3 = L[j+1][1];
		y3 = L[j+1][2];
		d = (x1*y2)+(y1*x3)+(x2*y3)-(x3*y2)-(x2*y1)-(x1*y3);
		if(d < 0,  \\ (x2,y2) not in convex hull, lose it
		   L = remove_elt(L,j);
		   l--;
		   if(j>2,j= j-2,j--); /* move back to test earlier point */
		);
		if(j <= 0 || j >=(l-1), break);
	);
	L
}

{
const_term(R,w,p,prec, aa,ii,bb,cc)=
if(DEBUG>=10,print("in const_term, R = ",R,"\nw=",w));
        aa=lift(polcoeff(R,0));
	bb = 0;
        for(ii=0,poldegree(w),
	    cc = polcoeff(aa,ii);
	    if(type(cc) == "t_PADIC",cc = truncate(cc));
	    bb = bb + cc*Z^ii;
        );
	bb = Mod(bb,w);
	bb
}

{
get_term(R,i,w,p,prec, aa,ii,bb,cc)=
if(DEBUG>=10,print("in get_term, R = ",R,"\ni=",i));
        aa=lift(polcoeff(R,i));
	bb = 0;
        for(ii=0,poldegree(w),
	    cc = polcoeff(aa,ii);
	    if(type(cc) == "t_PADIC",cc = truncate(cc));
	    bb = bb + cc*Z^ii;
        );
	bb = Mod(bb,w);
	bb
}

{
make_unpadic(R,w,p,prec, S,ii)=
if(DEBUG>=10,print("in make_unpadic, R = ",R));
    S = 0;
    for(ii=0,poldegree(R),
        S += get_term(R,ii,w,p,prec)*x^ii);
    return(S);
}

{
leading_term(R,w,p,prec, aa)=
        aa=lift(polcoeff(R,poldegree(R,y),y));
        aa
}

{
normalize(R,w,p,prec, a0,vv,i)=
        a0 = const_term(R,w,p,prec);
        vv=0;
	for(i=0,poldegree(R),
	    if(polcoeff(R,i) != 0, vv = vv + polcoeff(R,i)/a0*y^i);
	);
	vv
}

{
make_monic(R,w,p,prec, an,vv,i,pr)=
        an = leading_term(R,w,p,prec);
        an = polcoeff(R,poldegree(R,y),y);
        vv=0;
	for(i=0,poldegree(R,y),
/*	    pr = valuation(polcoeff(R,i,y),p)+prec;*/
	    if(polcoeff(R,i,y) != 0, vv = vv + polcoeff(R,i,y)/an*y^i);
	);
	vv
}

{
make_padic(F,w,p,prec, FF,i,cf)=
        FF=0;
	for(i=0,poldegree(F),
	cf = polcoeff(F,i);
	if(type(cf) == "t_POLMOD", cf = lift(cf,Z));
	    if(polcoeff(F,i) != 0, FF = FF + (Mod(cf,w)*x^i));
/*	    if(polcoeff(F,i) != 0, FF = FF + (Mod(cf+O(p^prec),w)*x^i));*/
	);
	FF
}

{
make_modw(F,w,p,prec, FF,i,cf)=
        FF=0;
	for(i=0,poldegree(F),
	cf = polcoeff(F,i);
	if(type(cf) == "t_POLMOD", cf = lift(cf,Z));
	    if(polcoeff(F,i) != 0, FF = FF + (Mod(cf,w)*x^i));
	);
	FF
}


{
newtonfac(R,w,p,prec,  L,j,x1,x2,x3,y1,y2,y3,B,C,u,v,h,l)=
if(DEBUG>=1,print("in newtonfac, R=",R));
        if(polcoeff(R,0) == 0, 
	    print("found factor y");  
	    return(y);
	);
        R=normalize(R,w,p,prec);
	L=newtondiag(R,p,prec);
	for(j=1,(length(L)-2),
		x1 = L[j][1];
		y1 = L[j][2];
		x2 = L[j+1][1];
		y2 = L[j+1][2];
		x3 = L[j+2][1];
		y3 = L[j+2][2];

		if(((x2-x1)*(y3-y1)) != ((x3-x1)*(y2-y1)),
\\ not colinear, factor R
                        r = L[j+1][1];
			s = -L[j+1][2];
			B=0;
			for(k=0,r,
			    B = B + polcoeff(R,k)*y^k;
			);

			B = clean(B,w,p,prec);
			C = 1;
			u = 1;
			v = 1 - B;
			h = 0;
			l = hensel(R,B,C,u,v,h,r,s,w,p,prec);
			l = subst(l,x,y);
			l = make_monic(l,w,p,prec);
			return(l);
		);
	);
	return(0);
}

{
is_straight(L)=
	for(j=1,(length(L)-2),
		x1 = L[j][1];
		y1 = L[j][2];
		x2 = L[j+1][1];
		y2 = L[j+1][2];
		x3 = L[j+2][1];
		y3 = L[j+2][2];
		if(((x2-x1)*(y3-y1)) != ((x3-x1)*(y2-y1)),  \\ not colinear
			return(0);
		);
	);
	return(1);
}

{
coef_reduce(coef,w,p)=
        tt = lift(coef,Z);
	qq = 0;
        for(ii=0,poldegree(w),
	    qq = qq + Mod(polcoeff(tt,ii),p)*Z^ii);
	return(Mod(qq,w))
}

\\ convert R to a polynomial in the residue class field 
{
reduce(R,w,p)=
      rbar = 0;
      for(j=0,poldegree(R),
           e = coef_reduce(polcoeff(R,j),w,p);
                rbar = rbar + e*y^j;
      );
      rbar
}

\\ convert R in the residue class field back to Q_p
{
my_lift(R,s,r,w,p,prec,  e,rbar)=
      rbar = 0;
      for(j=0,poldegree(R),
           e = Mod(lift(lift(polcoeff(R,j),dummy),Z),w);
                rbar = rbar + (e*p^((-j)*s)+O(p^prec))*y^(j*r);
      );
      rbar
}

/* go to splitting field of u.  Right now I'm cheating by using pari
   routines to do this for a number field.  Eventually I need to
   do this right for the finite field/p-adic field case, or prove that
   this will always work */
{
lift_ext(w,u)=
u=lift(lift(u));
      k=nfinit(w);
      if(poldegree(w)==1,return(subst(u,y,Z)));
      return(subst(rnfequation(k,u),y,Z));
}

{
/* reduce p-adic polynomial to residue class field */
get_star(R,w,p,prec, r)=
	rstar=0;
        L=newtondiag(R,p,prec);
        n=L[length(L)][1];
	m=L[1][2];
	if(length(L)>1,
	    slope = (L[length(L)][2] - L[1][2])/(L[length(L)][1] - L[1][1]),
	    slope = 1);
        r = denominator(slope);
	m = n/r;
	for(j=1,length(L),
	        xj = L[j][1];
		yj = L[j][2];
		rstar = rstar + polcoeff(R,xj)*y^(xj/r)/p^yj;
	);
	return(reduce(rstar,w,p));
}

{
my_resultant(F,AA,w,p,prec)=
    return(polresultant(F,y-AA,x));
}

{
get_norm(F,w,p,prec)=
    g = polresultant(lift(F),w,Z);
    g = make_padic(g,Z-1,p,prec);
    return(g);
}

{
check_rstar(F,AA,w,p,prec, R,L,n,m,g,rstar,f,WWW,f1,f2,i,w2)=
/* these were computed in Hensel_fac, but I'm redoing them here
   because I need F, and passing everything seemed silly */
/*        R = polresultant(F,y-AA,x);*/
        R = my_resultant(F,AA,w,p,prec);
if(DEBUG>=10,print("in check_rstar, R = ",R));
        L=newtondiag(R,p,prec);

        n=L[length(L)][1];
	m=L[1][2];

        if(gcd(n,m)==1, /* case 2d, irreducible */
	    if(poldegree(w)==1,
if(DEBUG>=2,print("in subcase 2d, ground field"));
	        return(F), /* base field, just return polynomial */
		/* extension field, need to take norm */
if(DEBUG>=2,print("in subcase 2d, extension field"));
		return(F);
		));

	rstar = get_star(R,w,p,prec);
	f=factorff(rstar,p,w);
        len=matsize(f);
        if((len[1] == 1), \\ only one factor
	   if((poldegree(f[1,1]) == 1),
\\ subcase 3, power of a linear factor
if(DEBUG>=2,print("in subcase 3"));
              return(0),
\\ subcase 2c, power of an irreducible of degree > 1
if(DEBUG>=2,print("in subcase 2c"));
              w2 = lift_ext(w,f[1,1]);
              f2 = find_factor(F,w2,p,prec);
              f2 = get_norm(f2,w2,p,prec);
              return(f2);
           ),
\\ subcase 2b, two relatively prime factors
if(DEBUG>=1,print("in subcase 2b, f = ",f));
\\ choose the smallest-degree factor
           dmin = poldegree(f[1,1])*f[1,2];
	   minterm = 1;
	   for(i=2,len[1],
	       d = poldegree(f[i,1])*f[i,2];
	       if(d < dmin,
	           dmin = d;
		   minterm = i;
	       );
	   );
           f1 = f[minterm,1]^f[minterm,2];
	   f2 = 1;
	   for(i=1,len[1],
	       if(i != minterm, f2 = f2*(f[i,1]^f[i,2])));
	   WWW = cor6_34(R,f1,f2,w,p,prec);
	   WWW = subst(WWW,y,AA);
	   return(WWW);
        );
print("in check_rstar, should not get here");
	return(rstar);
}

{
hensel_fac(F,AA,w,p,prec, WWW, jjj,R,L)=
/*      R = polresultant(F,y-AA,x);*/
      R = my_resultant(F,AA,w,p,prec);
if(DEBUG>=1,print("in hensel_fac, \nF=",F,"\nAA = ",AA,"\nR=",R,"\n"));
      if(const_term(R,w,p,prec)==0,
	  return(AA));
      L=newtondiag(R,p,prec);

      if(is_straight(L),
          WWW = check_rstar(F,AA,w,p,prec);
	  return(WWW),

if(DEBUG>=2,print("\nnewton diagram = ",L," is not straight\n"));
          WWW = newtonfac(R,w,p,prec);
	  WWW = subst(WWW,y,AA);
	  return(WWW);
      );
}


{
/* find root alpha of f in Q_p(w) with ord(alpha) = 0 */
get_alpha(f,w,p,prec)=
if(DEBUG>=10,print("in get_alpha, f=",f));
      LL = factorff(f,p,w);
      for(i=1,matsize(LL)[1],
          if(poldegree(LL[i,1])==1,
	      aa = -subst(LL[i,1],y,0);
	      if(myvaluation(aa,p)==0,return(make_padic(lift(aa,dummy),w,p,prec)))));
/* the lift brings us from the residue class field to K 
   the dummy variable makes pari lift the right way */
print("should not get here");
return(0);
}


{
/* reduce p-adic polynomial to residue class field */
/* this is get_star, but it also returns r,s,m,n */
get_all(R,w,p,prec)=
if(DEBUG>=10,print("in get_all, R = ",R));
	rstar=0;
        L=newtondiag(R,p,prec);
        n=L[length(L)][1];
	for(j=1,length(L),
	        xj = L[j][1];
		yj = L[j][2];
		rstar = rstar + polcoeff(R,xj)*y^xj/p^yj;
	);
	a0=myvaluation(polcoeff(R,0),p);
	m=gcd(a0,n);
	rr=n/m;
	ss=a0/m;
	return([get_alpha(reduce(rstar,w,p),w,p,prec),rr,ss,m,n]);
}

{
get_evec(h,r,s,t,u, L)=
if(DEBUG>=10,print("in get_evec"));
     if(h==1,
         if(Mod(s[1],r[1])==0,
	     return([0]),
	     return([lift(Mod(u,r[1])/Mod(s[1],r[1]))])));

     L = xgcd(s[h]*t[h+1]/r[h],t[h+1]/t[h]);
     e = L[1]*u;
     v = L[2]*u;
     k = floor(e*t[h]/t[h+1]);
     e = e - k*t[h+1]/t[h];
     v = v - k*s[h]*t[h+1]/r[h];
     return(concat(get_evec(h-1,r,s,t,v),[e]));
}

{
find_factor(F,w,p,prec,  r,s,m,n,L,A,R,t,h,i,S,u,evec,e0,j,E,C,T,g,Alist)=
/*     F = make_padic(F,w,p,prec);*/
if(DEBUG>=10,print("in find_factor, F = ",F));
     if(polcoeff(F,0)==0,return(y));
     F = make_modw(F,w,p,prec);
     g=hensel_fac(F,x,w,p,prec);
     if(g!=0, return(g)); /* succeeded */
     /* now the hard part */
     L = get_all(F,w,p,prec);
     alpha = L[1];
     r = [L[2]];
     s = [L[3]];
     m = L[4];
     n = L[5];

     /* outer loop */
     A = [x];
     R = [F];
     t = concat([1],r);

     for(h=1,9999,
          /* inner loop */
	  B = A[h]^(t[h+1]/t[h]);
          S = my_resultant(F,B,w,p,prec);
	  u = [s[h]*t[h+1]^2/(r[h]*t[h])];

	  for(i=0,9999,
	      evec = get_evec(h,r,s,t,u[i+1]);
	      e0 = u[i+1]/t[h+1];
	      for(j=1,h,e0 = e0 - evec[j]*s[j]/r[j]);
	      E=p^e0;
	      for(j=1,h,E = E*A[j]^evec[j]);
              /* Step 8 */
	      C = lift(Mod(B,F)/Mod(E,F));
	      T = my_resultant(F,C,w,p,prec);
	      g = hensel_fac(F,C,w,p,prec);
	      if(g!=0, return(subst(g,x,C)));

	      Alist = get_all(T,w,p,prec);
	      alpha = Alist[1];
	      Bnew = B - alpha*E;
	      Snew = my_resultant(F,Bnew,w,p,prec);
	      g = hensel_fac(F,Bnew,w,p,prec);
	      if(g!=0, 
	          return(subst(g,y,Bnew)));

	      Alist = get_all(Snew,w,p,prec);
	      ordB = Alist[3]/Alist[2];
	      if((t[h+1] % denominator(ordB))==0,
	           B = Bnew;
		   S = Snew;
		   u = concat(u,[ordB*t[h+1]]),
		   /* else */
		   break;
	      );
	  );
	      /* end inner loop */
	  /* Step 11 */
	  snew = numerator(ordB);
	  rnew = denominator(ordB);
	  A = concat(A,[Bnew]);
	  R = concat(R,[S]);
	  r = concat(r,[rnew]);
	  s = concat(s,[snew]);
	  t = concat(t,[lcm(t[h+1],r[h+1])]);
	  if(t[h+2] >= n, /* Step 11(b) */
	      evec = get_evec(h+1,r,s,t,1);
	      e0 = 1/t[h+2];
	      for(j=1,h+1,e0 = e0 - evec[j]*s[j]/r[j]);
	      E=p^e0;
	      for(j=1,h+1,E = E*(A[j]^evec[j]));
	      g = hensel_fac(F,E,w,p,prec);
	      if(g!=0, return(g));
	      print("should never get here");
	  );
      );
}

/* check if polynomials are same, up to available precision */
{
compare_poly(f1,f2,p,prec, i,diff,pr)=
    if(poldegree(f1) != poldegree(f2),return(0));
    for(i=0,poldegree(f1),
	diff = polcoeff(f1,i)-polcoeff(f2,i);
        pr = min(prec,padicprec(diff,p));
        if(valuation(diff,p) < pr,
	    return(0)));
    return(1);
}

{insert_poly(f,list,p,prec,  i,l)=
if(DEBUG>=10,print("in insert_poly, f = ",f));
    l = matsize(list)[1];
    for(i=1,l,
        if(compare_poly(f,list[i,1],p,prec)==1,
	    list[i,2]++;
	    return(list);
        );
    );
    list = concat(list~,[f,1]~)~;
    return(list);
}

{
my_factorpadic(F,p,in_prec,out_prec, g, fac_list, df,repeats)=
    df = deriv(F);
    repeats = gcd(F,df);
    if(repeats != 1,  /* nontrivial gcd, factor separately */
       fac_list = my_factorpadic(repeats,p,in_prec,out_prec);
       if(fac_list == 0, return(0));
       F = F/gcd(F,df),
       fac_list = [;]; /* no repeated roots */
    );
    w=Z-1;  
    while(true,
        g = subst(lift(find_factor(F,w,p,in_prec),Z),y,x);
        g = Zp_gcd(F,g,w,p,prec);
	g = g/polcoeff(g,poldegree(g));  /* make it monic */
	/* precision was too low, start over with higher precision */
	if(get_precision(g,w,p) < out_prec,
	    return(0);
	);
	F = F/g;
	g = g*(1+O(p^out_prec));  
	fac_list = insert_poly(g,fac_list,p,out_prec);
	if(type(F) != "t_POL",
	    print("error: factor ",g," of ",F);
	    break);
	if(poldegree(F)==0, break);
    );
    return(fac_list);
}

{
mf(F,p,out_prec, in_prec,l)=
    in_prec = out_prec;
    while(1,
        l = my_factorpadic(F,p,in_prec,out_prec);
	if(l==0,
	    print("raised precision to ",in_prec*2);
	    in_prec = in_prec*2, /* not enough precision, try again */
	    return(l);
	);
    );
}