Z:=IntegerRing(); 
Zxy<x,y>:=PolynomialRing(Z,2);

/* *********************************** */

IsIndecomp:=function(f);
/* Prend un polynome f de Zxy<x,y>:=PolynomialRing(IntegerRing(),2) ou de PolynomialRing(RationalField(),2) et dit si ce polynome est indecomp ou non */
/*  nous faisons le test sur le polytope de Newton
FALSE SIGNIFIE : JE NE SAIS PAS*/
t:=Cputime();
Qxy<x,y>:=PolynomialRing(RationalField(),2)!f;
Z:=IntegerRing();
ff:=f-Coefficient(Coefficient(f,1,0),2,0)+1;
N:=AllVertices(NewtonPolygon(ff));
S:=[];
for i:=1 to #N do
S:=Append(Append(S,Z!N[i][1]),Z!N[i][2]);
end for;
bool_gcd:=GCD(S);
if (bool_gcd eq 1) then return true; else return false; end if;
end function;

/*  ************************************************************ */



/* ********************************************************************* */
fabriq_poly:=function(n,N);
/* fabriq des polynomesde dgre total n avec des coefficients pris au hasard dans [-N,N] */
Zxy<x,y>:=PolynomialRing(IntegerRing(),2);
f:=Zxy!0;
for k:=0 to n do
  for l:=0 to n-k do
  a:=Random(-N,N);
   f:=f+a*x^k*y^l; 
end for;
end for;

return f;
end function;

/* ************************************************************************ */






 /* ******************************************************************** */
stat_direct_complet:=function(n,nb,N);
/* Effectue IsIndecomp sur nb poly de degré total n avec coeff dans [-N,N] */
compt:=0;
temps_moy:=RealField()!0;
temps_min:=RealField()!0;
temps_max:=RealField()!0;


f:=fabriq_poly(n,N);
t:=Cputime();
bool:=IsIndecomp(f);
t:=Cputime()-t;
temps_moy:=temps_moy+t;
temps_min:=t;
temps_max:=t;
if (bool eq true) then compt:=compt+1;end if;


for i:=2 to nb do
f:=fabriq_poly(n,N);
t:=Cputime();
bool:=IsIndecomp(f);
t:=Cputime()-t;
if t le temps_min then temps_min:=t; end if;
if t ge temps_max then temps_max:=t; end if;
temps_moy:=temps_moy+t;
if (bool eq true) then compt:=compt+1;end if;
end for;
temps_moy:=temps_moy/nb;
return "Nb de succes=",compt,"Temps moyen=",temps_moy,"temps_max=", temps_max,"temps_min=",temps_min;
end function;

/* ********************************************************************* */
fabriq_poly_creux:=function(n,N,prop);
/* fabriq des polynomesde dgre total n avec des coefficientsoù il ya prop fois plus de coefficients nuls que de coeff non nuls pris au hasard dans [-N,N] */
Zxy<x,y>:=PolynomialRing(IntegerRing(),2);
f:=Zxy!0;
for k:=0 to n do
  for l:=0 to n-k do
  creux:=Random(prop);
if (creux eq 0) then
  a:=Random(-N,N);
   f:=f+a*x^k*y^l; 
end if;
end for;
end for;

return f;
end function;

/* ************************************************************************* */
stat_direct_complet_creux:=function(n,nb,N,prop);
/* Effectue IsIndecomp sur nb poly CREUX de degré total n avec coeff dans [-N,N] */

compt:=0;
temps_moy:=RealField()!0;
temps_min:=RealField()!0;
temps_max:=RealField()!0;

f:=fabriq_poly_creux(n,N,prop);
t:=Cputime();
bool:=IsIndecomp(f);

t:=Cputime()-t;
temps_moy:=temps_moy+t;
temps_min:=t;
temps_max:=t;
if (bool eq true) then compt:=compt+1;end if;


for i:=2 to nb do
 f:=fabriq_poly_creux(n,N,prop);
t:=Cputime();
bool:=IsIndecomp(f);

t:=Cputime()-t;
if t le temps_min then temps_min:=t; end if;
if t ge temps_max then temps_max:=t; end if;
temps_moy:=temps_moy+t;
if (bool eq true) then compt:=compt+1;end if;
end for;
temps_moy:=temps_moy/nb;
return "Nb de succes=",compt,"Temps moyen=",temps_moy,"temps_max=", temps_max,"temps_min=",temps_min;

end function;

/* ************************************************************************* */

/* *********************************************************************** */

/* ********************************************************************* */

/* ********************************************************************** */

/* ***************************************************************** */
IsIndecomp_mod:=function(f);
/* Prend un polynome f de Zxy<x,y>:=PolynomialRing(IntegerRing(),2) et dit si ce polynome est indecomp ou non en utilisant la méthode polytope Newton modulaire*/
temps:=Cputime();
bool_abs:=false;

Z:=IntegerRing();
Zxy<x,y>:=PolynomialRing(Z,2);
L:={};

ff:=f-Coefficient(Coefficient(f,1,0),2,0)+1;
N:=AllVertices(NewtonPolygon(ff));

for i:=1 to #N do
L1:=Factorization( Z!( Coefficient( UnivariatePolynomial(Coefficient(f,1,Z!N[i][1])),Z!N[i][2])));
for j:=1 to #L1 do
  Include(~L,L1[j][1]);
end for;
end for;

L:=IndexedSetToSequence(SetToIndexedSet(L));
Sort(~L);
long:=#L;
print"L=",L;



bool_gcd:=false;

i:=1;
p:=1;
F:=f;
while (bool_abs eq false) and (i le long) do
p:=L[i];
Fp:=FiniteField(p);
Fpxy<x,y>:=PolynomialRing(Fp,2);
F:=Fpxy!ff;
if (TotalDegree(F) eq TotalDegree(ff)) then

N:=AllVertices(NewtonPolygon(F));
print"p=",p;
print"N=",N;
S:=[];
for j:=1 to #N do
n1:=Z!N[j][1];
n2:=Z!N[j][2];
if (n1 ne 0) then S:=Append(S,n1); end if;
if (n2 ne 0) then S:=Append(S,n2); end if;
end for;
/*print"S=",S;*/
if (#S ne 0) then
  if (GCD(S) eq 1) then bool_abs:=true; end if;
    
end if;

end if;

i:=i+1;
end while;
temps:=Cputime()-temps;
print"Ce test a necessite", temps, "secondes";

return bool_abs;
end function;
/* ***************************************************************** */
IsIndecomp_mod_2:=function(f);
/* Prend un polynome f de Zxy<x,y>:=PolynomialRing(IntegerRing(),2) et dit si ce polynome est indecomp ou non en utilisant la méthode polytope Newton modulaire*/


/*  COMME PRECEDEMMENT MAIS SANS LES PRINT   */

temps:=Cputime();
bool_abs:=false;

Z:=IntegerRing();
Zxy<x,y>:=PolynomialRing(Z,2);
L:={};

ff:=f-Coefficient(Coefficient(f,1,0),2,0)+1;
N:=AllVertices(NewtonPolygon(ff));

for i:=1 to #N do
L1:=Factorization( Z!( Coefficient( UnivariatePolynomial(Coefficient(f,1,Z!N[i][1])),Z!N[i][2])));
for j:=1 to #L1 do
  Include(~L,L1[j][1]);
end for;
end for;

L:=IndexedSetToSequence(SetToIndexedSet(L));
Sort(~L);
long:=#L;




bool_gcd:=false;

i:=1;
p:=1;
F:=f;
while (bool_abs eq false) and (i le long) do
p:=L[i];
Fp:=FiniteField(p);
Fpxy<x,y>:=PolynomialRing(Fp,2);
F:=Fpxy!ff;
if (TotalDegree(F) eq TotalDegree(ff)) then

N:=AllVertices(NewtonPolygon(F));

S:=[];
for j:=1 to #N do
n1:=Z!N[j][1];
n2:=Z!N[j][2];
if (n1 ne 0) then S:=Append(S,n1); end if;
if (n2 ne 0) then S:=Append(S,n2); end if;
end for;

if (#S ne 0) then
  if (GCD(S) eq 1) then bool_abs:=true; end if;
    
end if;

end if;

i:=i+1;
end while;
temps:=Cputime()-temps;


return bool_abs;




end function;
/* ******************************************************************* */

/* *************************************************************************/
stat_mod_complet2:=function(n,nb,N);
/* test sur des exemples la fonction IsIndecomp_mod*/
compt:=0;
comptR:=0;
temps_moy:=RealField()!0;
temps_moyR:=RealField()!0;
Zxy<x,y>:=PolynomialRing(IntegerRing(),2);


f:=Zxy!0;
for k:=0 to n do
  for l:=0 to n-k do
  a:=Random(-N,N);
   f:=f+a*x^k*y^l; 
end for;
end for;
/*print"f=",f;*/
t:=Cputime();
bool_irred_abs:=IsIndecomp_mod_2(f);
t:=Cputime()-t;
temps_moy:=temps_moy+ t;
if (bool_irred_abs eq true) then
compt:=compt+1; end if;

t_min:=t;
t_max:=t;
temps_moy:=t;


for i:=2 to nb do
f:=Zxy!0;
for k:=0 to n do
  for l:=0 to n-k do
  a:=Random(-N,N);
   f:=f+a*x^k*y^l; 
end for;
end for;
/*print"f=",f;*/
t:=Cputime();
bool_irred_abs:=IsIndecomp_mod_2(f);
t:=Cputime()-t;
temps_moy:=temps_moy+ t;
if (t ge t_max) then t_max:=t; end if;
if (t le t_min) then t_min:=t; end if;
if (bool_irred_abs eq true) then
compt:=compt+1; end if;


end for;
temps_moy:=temps_moy/nb;

print "NB de succes_IsIndecomp_mod=",compt,"temps_moy=",temps_moy,"temps max=",t_max,"temps min=",t_min;

return compt;

end function;
/* *********************************************************************** */
remplir_tableau_test_mod:=function(nb);

print"Test avec degre 10";
print"";
stat_mod_complet2(10,nb,10^12);
print"";
print"";
print"Test avec degre 30";
print"";
stat_mod_complet2(10,nb,10^12);
print"";
print"";
print"Test avec degre 50";
print"";
stat_mod_complet2(50,nb,10^12);
print"";
print"";
print"Test avec degre 100";
print"";
stat_mod_complet2(100,nb,10^12);
print"";
print"";
print"Test avec degre 200";
print"";
stat_mod_complet2(200,nb,10^12);
return nb;
end function;
/* ************************************************* */
remplir_tableau_test_direct_creux:=function(nb);

print"Test avec degre 10, prop=0";
print"";
stat_direct_complet_creux(10,nb,10^12,0);
print"";
print"";
print"Test avec degre 10, prop=1";
print"";
stat_direct_complet_creux(10,nb,10^12,1);
print"";
print"";
print"Test avec degre 10, prop=2";
print"";
stat_direct_complet_creux(10,nb,10^12,2);
print"";
print"";
print"Test avec degre 10, prop=10";
print"";
stat_direct_complet_creux(10,nb,10^12,10);
print"";
print"";
print"Test avec degre 100, prop=2";
print"";
stat_direct_complet_creux(100,nb,10^12,2);
print"";
print"";
print"Test avec degre 200, prop=2";
print"";
stat_direct_complet_creux(200,nb,10^12,2);
return nb;
end function;
/* ************************************************* */
faire_test:=function(nb);

print"Je remplis le tableau avec des poly creux";
print"";
remplir_tableau_test_direct_creux(nb);

print"";
print"";
print"";

print"Je remplis le tableau avec le test modulaire";
print"";
remplir_tableau_test_mod(nb);

return nb;
end function; 
/* *************************************** */
