multifac_cont_improve:=proc(f,ord,e)
  local f0,x,cont1,prim_f,FF1,FF,VARS,cont2;
  f0:=evalf(expand(f/norm(f,infinity)));
  x:=ord[2];
  # The following is to used to compute the content of the polynomial
  cont1:=cont_mul(f0,ord,x,e);
  prim_f:=mcofactor(f0,cont1,ord,e);
  prim_f:=cutoff(prim_f,ord,e);
  VARS:=[op(indets(prim_f))];
  x:=ord[1];
  cont2:=cont_mul(prim_f,ord,x,e);
  prim_f:=mcofactor(prim_f,cont2,ord,e);
  prim_f:=cutoff(prim_f,ord,e);
  FF1:=multifac(prim_f,VARS);
  FF1:= [seq(seq(FF1[i][1],j=1..FF1[i][2]), i=1..nops(FF1))]:
  FF1:=improve_factorization(prim_f,FF1,VARS):
  FF:=[cont1,cont2,op(FF1)];
  RETURN(FF);
end proc;