read("multifac_1.3.mpl"):
infolevel[appfac]:=0;

randp:=proc(vars::list               # the variables 
           ,list_deg_factor::list    # the list of the degree of the factors   
           ,N::integer               # absolute maximum of coefficients of the factors
           ,e::integer               # the perturbation size: 10^(-e)
           ,dn)                      # Optional: number between 0 and 1 giving the how dense the perturbation should be compared to product default = 1.
 
       local i,f0,fac,pert_pol,res, deg_pert, dd;

if nops([args]) = 4 then dd:=1 else dd:=dn end if;

       deg_pert := convert(list_deg_factor, `+`);
       
       f0:=1;
       for i from 1 to nops(list_deg_factor) do                                                                         #fac:=10^(-7)*randpoly(vars,degree=list_deg_factor[i],coeffs=rand(-N*10^7..N*10^7),dense);
fac:=randpoly(vars,degree=list_deg_factor[i],coeffs=rand(-N..N),dense);          

f0:=f0*(fac);
       end do;

dd:=round(nops(convert(expand(f0),list))*dd);

print('noise_terms' = dd);
     pert_pol:=randpoly(vars,degree=deg_pert,coeffs=rand(-10^(e)..10^(e)), terms=dd); 

       return(f0,'`*`'(10^(-e)*norm(expand(f0),2)/norm(pert_pol, 2),pert_pol));
       #return(res);

end proc;


N:=200;

agr:=[9,7],5,4;

times6 := vector(N,0);
#times6:=convert(times6,list);

for i from 1 to N do
  cleanF6, noiseF6:=randp([x,y],agr,.25):
  F6 := expand(evalf(cleanF6 + noiseF6)):
try
  times6[i]:=time():
  facF6:=appfac(F6,[x,y]):
  times6[i]:=time()-times6[i];
  numfacs6[i]:=nops(facF6);
catch:
  times6[i]:=0;
end try;

tmp:=select(x->x<>0,convert(times6,list));
print('i'=i, 'Min'=min(op(tmp)), 'Med'=stats[describe, median](tmp), 'Max'=max(op(tmp)))

 end do:

times6:=convert(times6,list);
times6:=select(x->x<>0,times6);

stats[describe,mean](times6);
stats[describe,median](times6);