
`UP/Sqrfree` := proc() local P,R,U,env;

	P := args[1];
	if not hasCategory(P,UnivariatePolynomial) then
		ERROR(`argument must be a UnivariatePolynomial`) fi;
	R := P[CoefficientRing];
	env:=['C'=R, 'D'=P,'E'=U];

	# Sqrfree decomposition

	if hasCategory(R,FiniteField) then

	   # Musser's algo.  cf. [1]  algo 4.3 (non recursive form)

	   U[`p_expo`]:=subs(env, proc(a)     # highest power of p dividing
					      # all the exponents of a
	   local i,listcoef,ip,listexpo,g,p;

	   p:=D[Characteristic];
	   listcoef:=D[ListCoeffs](a);
	   i:=0;
	   listexpo:=[];

	   while nops(listcoef) > 0 do
	      if C[`<>`](op(1,listcoef),C[0]) then
		 listexpo:=[op(listexpo),i];
              fi;
	      listcoef:=[op(2..nops(listcoef),listcoef)];
	      i:=i+1
       	   od;

	   g:=igcd(op(listexpo));

	   ip:=0;
	   while irem(g,p) = 0 do
	      ip:=ip+1;
	      g:=iquo(q,p);
           od;

	   ip;

	   end);

           U[`p^n-th Root`] := subs(env,proc(pol,n)
	   local i,coefar,root,p,q;

	   coefar:=D[ArrayCoeffs](pol);
	   root:=D[0];
	   p:=D[Characteristic];
	   q:=p^n;

	   for i from 0 by q to D[degree](pol) do  # generic
	      root:=D[`+`](root,D[monomial](C[`p^n-th Root`](coefar[i],n),
					    iquo(i,q)));
	   od;      
           root;
	   end);

           U[Sqrfree]:=subs(env, proc(pol)
	   local lc,a,b,c,z,y,w,i,ip,ip1,Output,p;

	   p:=D[Characteristic];
	   lc:=D[Lcoeff](pol);
	   a:=D[`R/`](pol,lc);
	   Output:=[];
	   b:=D[Diff](a);
           ip:=1;

	   do
	   if D[`<>`](b,D[0]) then  # a is not a p-th power
	     i:=1;
	     c:=D[Gcd](a,b);
	     w:=D[Div](a,c);
	     while D[`<>`](w,D[1]) do
       	        y:=D[Gcd](w,c);
	        z:=D[Div](w,y);
		if D[Degree](z) <> 0 then
                   Output:=[op(Output),[z,i*ip]];
                fi;
	        w:=y;
	        c:=D[Div](c,y);
	        i:=i+1;
             od;
	     a:=c;
           fi;

           # a is a p-th power

           if D[`=`](a,D[1]) then break; fi;
	    
           # a is a non constant p-th power

	   ip1:=E[`p_expo`](a);
	   a:=E[`p^n-th Root`](a,ip1);
           b:=D[Diff](a);
	   ip:=ip*p*ip1;
           od;
	   
	   [lc,op(Output)];

           end);	


          elif R[Characteristic] = 0 and
	       hasCategory(R,UniqueFactorizationDomain) then

           # Yun's algorithm   cf.  [1] algo 8.2

	   U[Sqrfree]:=subs(env,proc(pol)
	   local a,b,c,i,g,w,y,z,cont,Output;

	   cont:=D[Content](pol,'a');
	   i:=1;
	   Output:=[];
	   b:=D[Diff](a);
	   c:=D[Gcd](a,b);

	   if D[`=`](c,D[1]) then w:=a
           else
              w:=D[Div](a,c);
              y:=D[Div](b,c);
	      z:=D[`-`](y,D[Diff](w));
	      while D[`<>`](z,D[0]) do
	         g:=D[Gcd](w,z);
		 if D[Degree](g) <> 0 then
	            Output:=[op(Output),[g,i]];
                 fi;
	         i:=i+1;
	         w:=D[Div](w,g);
	         y:=D[Div](z,g);
	         z:=D[`-`](y,D[Diff](w));
             od;
	  fi;

          if D[Degree](w) <> 0 then
	     [cont,op(Output),[w,i]];
          else
	     [cont,op(Output)]
          fi;

	  end);
	     
         else
       	    U[Sqrfree]:=notImplemented;
         fi;
	
	eval(U[Sqrfree]);
end:

save `Sqrfree.m`;
quit
