/*
  Absolute factorization of bivariate polynomials over Q.

  pre-version 0.1 alpha 1

  (c) G. Chèze, April 2005.

  Reference:

  @inproceedings{Cheze2004,
    author = {G. Ch{\`e}ze},
    title = {Absolute polynomial factorization in two variables
             and the knapsack problem},
    booktitle = {Proceedings of ISSAC 2004},
    year = {2004},
    isbn = {1-58113-827-X},
    pages = {87--94},
    location = {Santander, Spain},
    doi = {http://doi.acm.org/10.1145/1005285.1005300},
    publisher = {ACM Press},
  }

  WARNING: this package is experimental !

  Please report to Guillaume.Cheze[at]math.unice.fr

  Export:

intrinsic Ab_Fac_Knap(P,prec)

   P must be an irreducible polynomial in Z[x,y] monic in y.
   P must have no linear factor.
   prec must be 1
 Abs_Fac_Knap return
 K the number field generated by all the coefficients of one absolute factor P_i
P1 one absolute factor
and P.

*/


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

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

chgvar:=function(P,x0,lambda,x,y)
Ptemp:=Evaluate(P,1,x-x0+lambda*y);
return Ptemp;
end function;

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

choixfibre:=function(P)
Ptemp:=Evaluate(P,1,0);
Ptemp:=UnivariatePolynomial(Ptemp);
return Ptemp;
end function;
 
/* ******************************************************* */

estimprec:=function(n)
l:=Ceiling( n^(1/2) );
temp:=n+l;
C:= Ceiling(2^(1/2)/(temp^(1/2))*( (temp/2+(3/4*temp)^2 )^(temp/4)) );
eta:=Ceiling( Log(C)/Log(10));
C:=10^(eta);
eta:=eta+20;
return C,eta;
end function;

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

estimracine:=function(f,eta)
  L:=Roots(f:Digits:=eta);
return L;
end function;

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

estim_ai_bi:=function(p,CCestimy,n,L)
alpha:=Derivative(p,1);
beta:=Derivative(p,2);
gamma:=Derivative(alpha,1);
delta:=Derivative(beta,2);
epsilon:=Derivative(alpha,2);

alpha0:=Evaluate(alpha,1,0);
alpha0:=UnivariatePolynomial(alpha0);
alpha0:=CCestimy!alpha0;
beta0:=Evaluate(beta,1,0);
beta0:=UnivariatePolynomial(beta0);
beta0:=CCestimy!beta0;
gamma0:=Evaluate(gamma,1,0);
gamma0:=UnivariatePolynomial(gamma0);
gamma0:=CCestimy!gamma0;
delta0:=Evaluate(delta,1,0);
delta0:=UnivariatePolynomial(delta0);
delta0:=CCestimy!delta0;
epsilon0:=Evaluate(epsilon,1,0);
epsilon0:=UnivariatePolynomial(epsilon0);
epsilon0:=CCestimy!epsilon0;

A:=[];
B:=[];

for i:=1 to n do
 alpha00:=Evaluate(alpha0,L[i][1]);
 beta00:=Evaluate(beta0,L[i][1]);
 gamma00:=Evaluate(gamma0,L[i][1]);
 delta00:=Evaluate(delta0,L[i][1]);
 epsilon00:=Evaluate(epsilon0,L[i][1]);
 a:=-alpha00/beta00;
A:=A cat [a];
 b:=-(1/(2*beta00))*(gamma00+2*epsilon00*a+delta00*a^2);
 B:= B cat [b];
end for;
return A,B,alpha0,beta0,gamma0,delta0,epsilon0;
end function; 

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

 
estim_l_et_pretri:=function(L,A,B,n,eta) /*estime le nb de bi réel et rend la 
                                          liste des bi réels ainsi que la liste
                                          des yi et ai correspondant*/
temp:=0;
AR:=[];
BR:=[];
LR:=[];
for i:=1 to n do 
 if Abs(Im(B[i])) lt 10^(-eta+20) then temp:=temp+1;
                                  BR:=BR cat [B[i]];
                                  LR:=LR cat [L[i][1]];
                                  AR:=AR cat [A[i]];
 end if;
end for;
return temp,LR,AR,BR;
end function;

/* ********************************************************************* */
facteur_lineaire:=function(B,n,eta)

/*test si tous les bi sont nuls dans ce cas cela signifierait que nous avond des facteurs lineaires
rend true si il est possible d'avoir des facteurs lineaires*/
  compt:=0;
for i:=1 to n do 
  if (Abs(B[i])) lt 10^(-eta+20) then compt:=compt+1; end if;
end for;
if (compt eq n) then return true;
else return false;
end if;
end function;

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

tribi:=function(L,LR,A,AR,B,BR,n,eta) /*tri bi avec bi réel d'abord puis 
                                      complex imag>0 et ensuite les conj,
                                   rend la liste Ly et La qui est cohérente 
                                   avec cette numérotation*/
Lb:=[];
Lbconj:=[];
Ly:=[];
Lyconj:=[];
La:=[];
Laconj:=[];
Lb:=Lb cat [BR[i]: i in [1..#BR]];
Ly:=Ly cat [LR[i]: i in [1..#LR]];
La:=La cat [AR[i]: i in [1..#AR]];
for i:=1 to n do
 if Im(B[i]) gt 10^(-eta) then Lb:=Lb cat [B[i]];
                               Lbconj:=Lbconj cat [ComplexConjugate(B[i])];
                               Ly:=Ly cat [L[i][1]];
                               Lyconj:=Lyconj cat [ComplexConjugate(L[i][1])];
                               La:=La cat [A[i]];
                               Laconj:=Laconj cat [ComplexConjugate(A[i])];
 end if;
end for;
Ly:=Ly cat Lyconj;
La:=La cat Laconj; 
return Ly,La,Lb,Lbconj;
end function;
 
/* ******************************************************************* */

fabriq_mat_pour_LLL:=function(Lb,n,l,C)
temp:=#Lb;
L:=[Round(C*Lb[i]): i in [1..l]];
L:=L cat [Round(2*C*Re(Lb[i])): i in [l+1..#Lb]];
X:=RMatrixSpace(IntegerRing(),temp,temp+1)!0;
for i:=1 to temp do X[i,i]:=1; end for;
for i:=1 to temp do X[i,temp+1]:=L[i]; end for;
return X;
end function;

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

GramSchmidt:=function(M,n,c)    /*Matrice coeff entiers et dim de M=n*c */
G:=RMatrixSpace(RationalField(),n,c)!M;
for k:=2 to n do


    for j:=1 to k-1 do
       num:=0;
       denom:=0;	
         for i:=1 to c do
         num:=num+G[j,i]*M[k,i];
         denom:=denom+G[j,i]^2;
         end for;
       G[k]:=G[k]-(num/denom)*G[j];
    end for;
end for;
return G; 
end function;


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

nouv_GramSchmidt:=function(M,n,c)
/* M est une matrice a coeff entiers avec n lines et c columns*/
  G:=RMatrixSpace(RationalField(),n,c)!M;
for k:=1 to n do
d:=0;
  for i:=1 to c do
d:=G[k,i]*G[k,i]+d;
  end for;

for j:=k+1 to n do
num:=0;
  for i:=1 to c do
num:=num+G[j,i]*G[k,i];
end for;
for i:=1 to c do
  G[j,i]:=G[j,i]-(num/d)*G[k,i];
end for;
end for;
end for;
return G;
end function;

/* ******************************************************************** */
nouv_nouv_GramSchmidt:=function(LX,G)
/* LX  est la matrice avec en ligne les vecteurs LLL reduit,
G est la matrice de Gram,
On rend la matice avec en ligne les coord des vecteurs orthog de gramschmidt associe a LX*/
  M:=RMatrixSpace(RationalField(),NumberOfRows(LX),NumberOfColumns(LX))!LX;
for i:=2 to NumberOfRows(M) do
  /*M[i]:=LX[i];*/
  for j:=1 to i-1 do
  M[i]:=M[i]-G[i,j]/G[j,j]*M[j];
    end for;
end for;
return M;
end function;


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

Norm_vec_Gram_rend_k0:=function(G,l,c,M) /* prend une matrice où les lignes
                                            sont coord vect GramSchmidt,
                                           l lignes et c colonnes,
                                           M est une constante,
					rend  nb k0 tq pour tout k>k0 la norme
                                        au carreé de G[k] (kieme ligne de G) 
					soit  supérieur à M*/ 
compt:=0;
for i:=1 to l do
    n:=0;
    for j:=1 to c do
       n:=n+G[i,j]^2;
    end for;
   if (n gt M) then compt:=compt+1; 
               else compt:=0;
   end if;
end for;
k0:=l-compt;
return k0;
end function;

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

Form_Echelon_que_des_0_1:=function(M,m,n) /*on donne une matrice  et son nombre                                            de ligne m et son nombre col n*/
                                        /*on teste si la matrice est du type
                                         un 1 par colonne et que des 0*/

if (m eq 1) then 
   nb1:=0;
   for i:=1 to n do
     if (M[1,i] eq 1) then nb1:=nb1+1; end if;
   end for;
if (nb1 eq n) then b:=true;
              else b:=false;
end if;

else


b:=true;
i:=1;
while b and (i le n) do
nb0:=0;
nb1:=0;
   for j:=1 to m do
   if M[j,i] eq 0 then nb0:=nb0+1; end if;
   if M[j,i] eq 1 then nb1:=nb1+1; end if;
   end for;
if (nb0 ne m-1) or (nb1 ne 1) then b:=false; end if;
i:=i+1;
end while;
end if;
return b;
end function;

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

test_situation_generique_pour_les_relations:=function(RREF)
/*prend la matice RREF et regarde si une ligne n'est pas 0_1 vecteur et si la relation obtenu pouraait etre nulle*/
m:=NumberOfRows(RREF);
n:=NumberOfColumns(RREF);
bool:=true;
for i:=1 to m do
    compt:=0;
    for j:=1 to n do 
       if (Abs(RREF[i,j]) gt 1) then compt:=1; end if;
    end for;
    if (compt eq 1) then 
                        N:=0;
                        for j:=1 to n-1 do
                        N:=N+RREF[i,j];
                        end for;
                        if (N*(1.5) gt RREF[i,n]) then bool:=false; end if;
    end if;
end for;
return bool;
end function;

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

liste_0_1_vecteur_avec_bi_reel:=function(nbligne,l,M)
/*prend la 0-1 matrice M qui a nbligne ligne et nous sommes dans une situation avec l bi reel
rend L1 la liste des 0-1 vecteurs de M "contenant" des bi reel et rend L liste des 0-1 vecteurs "ne contenant pas" de bi reels*/

L1:=[];L:=[];
if (l ne 0) then 
	for i:=1 to nbligne do
	temp:=0;
		for  j:=1 to l do
		if (M[i,j] eq 1) then temp:=1; end if;
		end for;
	if (temp ne 0) then L1:=L1 cat [M[i]];
			else L:=L cat [M[i]];
	end if;
	end for;
else 
	for i:= 1 to nbligne do L:=L cat [M[i]]; end for;
end if;
return L1,L;
end function;  

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

trouve_deg_m_qd_l_est_diff_de_zero:=function(v,l,nbcol)
/* lorsque l est different de 0 alors grace a v=L1[1] un 0-1 vecteur "contenant" des bi reel, nbcol est le nb de colonne dans v;
 on rend m le degre des facteurs absolus*/

m:=0;
for i:=1 to l do 
	if (v[i] eq 1) then m:=m+1; end if;
end for;
for i:=l+1 to nbcol do 
	if (v[i] eq 1) then m:=m+2; end if;
end for;
return m;
end function;

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

verif_coherence_RREF_qd_l_est_diff_de_zero_et_rend_L2_L3:=function(L1,L,l,m,nbcol)
/*verif la coherence du deg m obtenu,
 rend un booleen qui dit si coherent, L2 et L3*/

L2:=[];L3:=[];

for j:=2 to #L1 do
v:=L1[j];
n:=0;
	for i:=1 to l do 
		if (v[i] eq 1) then n:=n+1; end if;
	end for;
	for i:=l+1 to nbcol do 
		if (v[i] eq 1) then n:=n+2; end if;
	end for;
if (n ne m) then return false,L2,L3; end if;
end for;

for j:=1 to #L do
v:=L[j];
n:=0;
	for i:=l+1 to nbcol do
		if (v[i] eq 1) then n:=n+2; end if;
	end for;
if (n eq m) then L2:=L2 cat [v]; end if;
if (n eq (2*m)) then L3:=L3 cat [v]; end if;
if(n ne m) and (n ne (2*m)) then return false,L2,L3; end if;
end for;
return true,L2,L3;
end function;

/* ****************************************************************** */
 decoup_en_2_somme_bi_complex:=function(Lbcomplex,C)
/* prend un ensemble de bi complexe de sommme nulle sans bi reel que l'on doit redecouper en deux sous sommes nulles
rend la matrice sur laquelle on a effectuee LLL afin d'obtenir deux 0_1 vecteurs*/

c:=#Lbcomplex;
X:=RMatrixSpace(IntegerRing(),c,c+2)!0;
for i:=1 to c do X[i,i]:=1;
                 X[i,c+1]:=Round(C*Re(Lbcomplex[i]));
                 X[i,c+2]:=Round(C*Im(Lbcomplex[i]));
end for;
F:=X*Transpose(X);
G1,T:=LLLGram(F:Sort:=false,Method:="Integral");
LX:=T*X;
G2:=nouv_nouv_GramSchmidt(LX,G1);
return LX,G2;
end function; 

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

coherence_RREF_rend_liste_0_1_vecteur_selon_type:=function(M,l,n,C,Lb)
/*Prend M RREF avec que des 0 et des 1
l nb bi reel,
n est le degre du polynome,
X est la matrice de depart pour le premier LLL dans la derniere colonne se trouve les bi bien ordonnés
dimX est le nombre de ligne de X
test la coherence avec prop de divisibilite donc rend un bool,
m le degre attendue pour les facteurs,
3 listes:
L1 liste des 0-1 vesteurs "avec" bi reel,
L2                         "sans" bi reel pas a redecomposer,
L3                                            a redecomposer.*/


nbligne:=NumberOfRows(M);
nbcol:=NumberOfColumns(M);
bool:=true;

if (l ne 0) then  L1,L:=liste_0_1_vecteur_avec_bi_reel(nbligne,l,M);
       		  m:=trouve_deg_m_qd_l_est_diff_de_zero(L1[1],l,nbcol);
		  bool,L2,L3:=verif_coherence_RREF_qd_l_est_diff_de_zero_et_rend_L2_L3(L1,L,l,m,nbcol);
end if;

if (l eq 0) then 
L1:=[];L2:=[];L3:=[];
L:=[];
test1:=0;
for i:=1 to nbcol do if (M[1,i] eq 1) then test1:=test1+1; end if; end for;

for i:=1 to nbligne do L:=L cat [M[i]]; end for;
list_norm:=[];
different:=false; /*different est un booleen qui dit s'il y a des 0-1 vect de 				norm differentes*/

for i:=1 to nbligne do
ntemp:=0;
	for j:=1 to nbcol do
		if (M[i,j]) eq 1 then ntemp:=ntemp+1; end if;
	end for;
list_norm:=list_norm cat [ntemp];
	if (i gt 1) and (list_norm[i] ne list_norm[i-1]) then different:=true;
							 test2:=list_norm[i];
							 test1:=list_norm[i-1];								end if;
end for;



if (different eq true) then 
	m:=Max([test1,test2]);
	for i:=1 to nbligne do
		if (list_norm[i] eq m) then L3:=L3 cat [M[i]]; end if;
		if (list_norm[i] eq (m div 2)) then L2:=L2 cat [M[i]]; end if;
		if (list_norm[i] ne m) and (list_norm[i] ne (m div 2)) then bool:=false; end if;

	end for;
else 


 
LB:=[];
for i:=1 to nbcol do
		  if (M[1,i] eq 1) then LB:=LB cat [Lb[i],ComplexConjugate(Lb[i])]; end if;
end for;
LXX,GG:=decoup_en_2_somme_bi_complex(LB,C);
    dimXX:=NumberOfRows(LXX);
    MM:=dimXX+((1.5)*dimXX)^2;
    kk:=Norm_vec_Gram_rend_k0(GG,dimXX,dimXX+2,MM);
print"kk=",kk;
if (kk eq 1) then  L2:=L; m:=2*test1;
		   L3:=[]; 
end if;
if (kk eq 2) then L3:=L;m:=test1;
		  L2:=[];
end if;
if (kk gt 2) then bool:=false; m:=0;end if;


end if;

end if;

if ( ((#L1)*m + (#L2)*m + (#L3)*m*2) ne n) then return false,0,L1,L2,L3;
else return bool,m,L1,L2,L3; end if;
end function;

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

fabriq_fact_pour_paquet_reel:=function(v,l,dim,Ly,La,Lb,CCestimxy,x,y)
/* v une matrice ligne extraite  de RREF et elle correspond a une somme bi nulle "comportant" des bi reels donc la sommme ne doit pas etre redecomposee on rend donc le facteur correspondant,
dim est le nombre de colonne de v,
l est le nombre de bi reel*/
F:=CCestimxy!1;
for j:=1 to l do
       if (v[j] eq 1) then F:=F*(y-Ly[j]-La[j]*(x)-Lb[j]*(x)^2); end if;
end for;
for j:=l+1 to dim do
       if (v[j] eq 1) then F:=F*(y-Ly[j]-La[j]*(x)-Lb[j]*(x)^2)*(y-ComplexConjugate(Ly[j])-ComplexConjugate(La[j])*(x)-ComplexConjugate(Lb[j])*(x)^2); end if;
end for;
return F;
end function;

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

list_fact_reel:=function(L1,l,n,Ly,La,Lb,CCestimxy,x,y)
/*n est le degre de P,
prend liste des 0-1 vecteurs avec bi reel rend la liste des facteurs correspondants*/

fact:=[];
if (l ne 0) then
dim:=NumberOfColumns(L1[1]);
for i:=1 to #L1 do
fact:=fact cat [fabriq_fact_pour_paquet_reel(L1[i],l,dim,Ly,La,Lb,CCestimxy,x,y)];
end for;
end if;
return fact;
end function;

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

list_fact_complex_bon_deg:=function(L2,l,n,Ly,La,Lb,CCestimxy,x,y)
/*prend liste des 0-1 vecteurs sans bi reel mais tel qu 'il n'est pas necessaire de  redecomp,
rend la liste des facteurs corespondants*/

fact:=[];
stop:=(n+l) div 2;
for i:=1 to #L2 do
 	F:=CCestimxy!1;
	for j:=l+1 to stop do 
		if (L2[i][j] eq 1) then F:=F*(y-Ly[j]-La[j]*(x)-Lb[j]*(x)^2)*(y-ComplexConjugate(Ly[j])-ComplexConjugate(La[j])*(x)-ComplexConjugate(Lb[j])*(x)^2); end if;
	end for;
fact:=fact cat [F];
end for;
return fact;
end function;

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

recommence_augmente_prec:=function(P,Px0,n,C,eta,LX,k0)
/*prend C,eta qui correspondent a la nouvelle precision souhaitee,
LX la matrice --------------------------------base du reseau ou la derniere coord est une combinaison des bi avec l'ancienne precision,
k0 le nombre de vecteurs a garder;
rend etad le nb de digits utlise pour les calculs,CCestim,CCestimy,CCestimxy,
l le nb de bi reel,
Ly,La,Lb,Lbconj la liste des yi,ai,bi avec nouvelle prec,
nouvLX le nouveau reseau ou la derniere coord est une combinaison lineaier des bi avec la "nouvelle" precision,
nouvk0 le nb de vect a garder dans le nouveau reseau*/

  etad:=eta;
CCestim<i>:=ComplexField();
CCestimy<y>:=PolynomialRing(CCestim);
CCestimxy<x,y>:=PolynomialRing(CCestim,2);
nouvPx0:=CCestimy!Px0;
print"Je calcule les racines";
L:=estimracine(nouvPx0,etad);         /* a ameliorer car n'utilise pas Newton*/
A,B:=estim_ai_bi(P,CCestimy,n,L);
l,LR,AR,BR:=estim_l_et_pretri(L,A,B,n,eta);
Ly,La,Lb,Lbconj:=tribi(L,LR,A,AR,B,BR,n,eta);
M:=RowSubmatrix(LX,1,k0);
c:=NumberOfColumns(M);
for i:=1 to k0 do
	temp:=0;
	for j:=1 to l do 
		temp:=M[i,j]*Round(C*Lb[j])+temp;
	end for;
	for j:=l+1 to c-1 do
		temp:=M[i,j]*Round(2*C*Re(Lb[j]))+temp;
	end for;
M[i,c]:=temp;
end for;
print"J'effectue LLL";
F:=M*Transpose(M);
G1,T:=LLLGram(F:Sort:=false,Method:="Integral");
nouvLX:=T*M;
print"Je calcule GramSchmidt";
G:=nouv_nouv_GramSchmidt(nouvLX,G1);
N:=(c-1)+(1.5*(c-1))^2;
nouvk0:=Norm_vec_Gram_rend_k0(G,k0,c,N);
return etad,CCestim,CCestimy,CCestimxy,l,Ly,La,Lb,Lbconj,nouvLX,nouvk0;
end function;

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

Newton_1_iteration:=function(Px0,yapp)
/*prend le polynome Px0 et une racine approchee yapp et effectue une iteration de Newton sur Px0 a partir de yapp*/
d:=Derivative(Px0);
denom:=Evaluate(d,y,yapp);
num:=Evaluate(Px0,y,yapp);
ymieux:=yapp-num/denom;
return ymieux;
end function;

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

recommence_augment_prec2:=function(P,Px0,Lycomplex,C,eta,alpha0,beta0,gamma0,delta0,epsilon0,m,LXfact,k1)
/* Lycomplex est une liste de yi correspondant a une liste de bi que l'on doit decouper en 2 pour obtenir 2 facteurs, chacun de degre m */
/*ameliore prec sur les yi par une iteration de Newton et en déduit les nouveaux ai,bi
puis avec ces nouveaux bi fabriq nouvelle matrice pour LLL fait LLL et gram Schmidt.
prend 
m est le degre des facteurs a obtenir
rend nouvLX le nouveau reseau,
nouvk1 le nb de vect a garder dans le nouveau reseau*/


CCestim<i>:=ComplexField();
CCestimy<y>:=PolynomialRing(CCestim);
CCestimxy<x,y>:=PolynomialRing(CCestim,2);
nouvPx0:=CCestimy!Px0;
k:=#Lycomplex;
Amieux:=[];Bmieux:=[];
Lymieux:=[];
for i:=1 to k do
ymieux:=Newton_1_iteration(Px0,Lycomplex[i]);
Lymieux:=Lymieux cat [ymieux];
 alpha00:=Evaluate(alpha0,ymieux);
 beta00:=Evaluate(beta0,ymieux);
 gamma00:=Evaluate(gamma0,ymieux);
 delta00:=Evaluate(delta0,ymieux);
 epsilon00:=Evaluate(epsilon0,ymieux);
 a:=-alpha00/beta00;
Amieux:=Amieux cat [CCestim!a];
 b:=-(1/(2*beta00))*(gamma00+2*epsilon00*a+delta00*a^2);
 Bmieux:= Bmieux cat [CCestim!b];
end for;
M:=RowSubmatrix(LXfact,1,k1);
c:=NumberOfColumns(M);
for i:=1 to k1 do
temp1:=0; temp2:=0;
	for j:=1 to c-2 do
	temp1:=M[i,j]*Round(C*Re(Bmieux[j]))+temp1;
	temp2:=M[i,j]*Round(C*Im(Bmieux[j]))+temp2;
	end for;
M[i,c-1]:=temp1;
M[i,c]:=temp2;
end for;
F:=M*Transpose(M);
G1,T:=LLLGram(F:Sort:=false,Method:="Integral");
nouvLXfact:=T*M;
G:=nouv_nouv_GramSchmidt(nouvLXfact,k1,c);
N:=(m)+(1.5*m)^2;
nouvk1:=Norm_vec_Gram_rend_k0(G,k1,c,N);
return eta,CCestim,CCestimy,CCestimxy,Lymieux,nouvLXfact,nouvk1;
end function;

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

liste_yi_ai_bi_corresp_0_1_vect:=function(Ly,La,Lb,v)
/*prend liste des yi, ai, bi et une 0_1 matrice v avec une seule ligne qui correspond a une somme bi sans bi reels et rend liste des yi,ai,bi, corespondant a ce 0_1 vecteur*/
 
c:=NumberOfColumns(v);
Lycomplex:=[]; Lacomplex:=[]; Lbcomplex:=[];
for i:=1 to c do 
  if (v[i] eq 1) then Lycomplex:=Lycomplex cat [Ly[i]];
                        Lacomplex:=Lacomplex cat [La[i]];
                        Lbcomplex:=Lbcomplex cat [Lb[i]];
end if;
end for; 
Lycomplex:=Lycomplex cat [ComplexConjugate(Lycomplex[i]): i in [1..#Lycomplex]];
Lacomplex:=Lacomplex cat [ComplexConjugate(Lacomplex[i]): i in [1..#Lacomplex]];
Lbcomplex:=Lbcomplex cat [ComplexConjugate(Lbcomplex[i]): i in [1..#Lbcomplex]];
return Lycomplex,Lacomplex,Lbcomplex;
end function;

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

liste_yi_ai_bi_corresp_0_1_vect_apres_decoup_en_2:=function(Lycomplex,Lacomplex,Lbcomplex,w)
/*prend liste yi,ai,bi tq bi donne somme nulle a redécomposer, w une 0_1 matrice qui donne la decomposition
rend Lyi,Lai,Lbi 6 listes correspondant aux yi,ai,bi d'une somme min sur les bi*/

Ly1:=[]; Ly2:=[]; La1:=[]; La2:=[]; Lb1:=[]; Lb2:=[];
c:=NumberOfColumns(w);
for i:=1 to c do
  if (w[1,i] eq 1) then Ly1:=Ly1 cat [Lycomplex[i]];
                        La1:=La1 cat [Lacomplex[i]];
                        Lb1:=Lb1 cat [Lbcomplex[i]];
                   else Ly2:=Ly2 cat [Lycomplex[i]];
                        La2:=La2 cat [Lacomplex[i]];
                        Lb2:=Lb2 cat [Lbcomplex[i]];
  end if;
end for;

return Ly1,La1,Lb1,Ly2,La2,Lb2;
end function;

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

donne_deux_facteurs:=function(Ly1,La1,Lb1,Ly2,La2,Lb2,CCestimxy,x,y)

F:=CCestimxy!1;G:=CCestimxy!1;

for i:=1 to #Ly1 do

F:=F*(y-Ly1[i]-La1[i]*(x)-Lb1[i]*(x)^2);
G:=G*(y-Ly2[i]-La2[i]*(x)-Lb2[i]*(x)^2);
end for;
fact:=[F,G];
return fact;
end function;

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

list_fact_complex_apres_redecomp:=function(P,n,Px0,C,eta,L3,l,m,Ly,La,Lb,alpha0,beta0,gamma0,delta0,epsilon0,CCestimxy,x,y)
/* prend liste des 0-1 vecteurs sans bi reel a redecomposer
rend un booleen pour dire si tout s'est bien passe, false si k1>2 
rend la liste des facteur correspondants*/

fact:=[];

for i:=1 to #L3 do

Lycomplex,Lacomplex,Lbcomplex:=liste_yi_ai_bi_corresp_0_1_vect(Ly,La,Lb,L3[i]);
LXfact,G:=decoup_en_2_somme_bi_complex(Lbcomplex,C);
dimLXfact:=NumberOfRows(LXfact);
N:=dimLXfact+(1.5*dimLXfact)^2;
k1:=Norm_vec_Gram_rend_k0(G,dimLXfact,dimLXfact+2,N);
RREFtemp:=EchelonForm(ColumnSubmatrix(RowSubmatrix(LXfact,1,k1),dimLXfact));
test01bis:=Form_Echelon_que_des_0_1(RREFtemp,k1,dimLXfact);
	while (test01bis eq false) do
	if (test_situation_generique_pour_les_relations(RREFtemp) eq false)
		then return false,fact;
	else eta,CCestim,CCestimy,CCestimxy,Lymieux,nouvLXfact,nouvk1:=recommence_augment_prec2(P,Px0,Lycomplex,C,eta,alpha0,beta0,gamma0,delta0,epsilon0,m,LXfact,k1);
print"recommence_prec2 est utilise";
		LXfact:=nouvLXfact;
		k1:=nouvk1;
		dimLXfact:=NumberOfRows(LXfact);
		nb_col_pour_RREF:=NumberOfColumns(LXfact)-2;
		RREFtemp:=EchelonForm(ColumnSubmatrix(RowSubmatrix(LXfact,1,k1),nb_col_pour_RREF));
		test01bis:=Form_Echelon_que_des_0_1(RREFtemp,k1,nb_col_pour_RREF);
	end if;
	end while;

if (k1 gt 2) then return false,fact;
	     else 
		Ly1,La1,Lb1,Ly2,La2,Lb2:=liste_yi_ai_bi_corresp_0_1_vect_apres_decoup_en_2(Lycomplex,Lacomplex,Lbcomplex,RREFtemp);
		fact:=fact cat donne_deux_facteurs(Ly1,La1,Lb1,Ly2,La2,Lb2,CCestimxy,x,y);
end if;

end for;

return true,fact;
end function;

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

fabriq_deterministe_poly_min:=function(j,List_fact,s,CCestimt,t)
/* prend un indice j, la liste des facteurs approchées List_fact, s le nombre de facteurs,
rend le poly min approché obtenu grace au coefficient de x^0y^j des facteurs,et la liste de ses racines*/


f:=CCestimt!1;
L:=[];
for k:=1 to s do
temp:=Coefficient(UnivariatePolynomial(Coefficient(List_fact[k],1,0)),j);
	f:=f*(t-temp);
L:=L cat [temp];
end for;
f:=CCestimt!f;
return f,L;
end function;

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

fabriq_aleatoire_poly_min:=function(List_fact,m,s,CCestimt,t)
/* fabrique un poly min approche de maniere aleatoire, en prenant une combinaison lineaire des termes de deg 0 en x, et rend aussi la liste de ces racines*/

list_comb_lin:=[];
f:=CCestimt!1;
L:=[];
for i:=1 to m do 
	list_comb_lin:=list_comb_lin cat [Random(-5,5)]; 
end for;

for k:=1 to s do
	temp:=0;
	for i:=1 to m-1 do
		temp:=temp+list_comb_lin[i]*Coefficient(Coefficient(List_fact[k],1,0),2,i);
	end for;
	f:=f*(t-temp);
L:=L cat [temp];
end for;
f:=CCestimt!f;
return f,L;
end function;

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

prend_poly_min_app_rend_exact:=function(f,s,Zt,t)
g:=Zt!0;
for i:=0 to s-1 do
temp:=Zt!Round(Re(Coefficient(f,i)));
g:=g+temp*t^i;
end for;
g:=t^s+g;
g:=Zt!g;
return g;
end function;

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

trouve_poly_min:=function(List_fact,m,s,Zt,CCestimt,t)
/* rend bool qui dit si on a trouve un poly irred et le poly irred correspondant,et la liste des racines du polynome*/

test_irred:=false;
j:=0;

while (test_irred eq false) and (m gt j) do
			      f,L:=fabriq_deterministe_poly_min(j,List_fact,s,CCestimt,t);
	g:=prend_poly_min_app_rend_exact(f,s,Zt,t);
	print"f=",f;
	print"g=",g;
	test_irred:=IsIrreducible(g);
	print"test_irred=",test_irred;
	j:=j+1;
end while;
if (test_irred eq true) then return true,g,L;
			else c:=0;
				while (test_irred eq false) and (5 gt c) do
				print"je fabrique un élément primitif de maniere aleatoire";
f,L:=fabriq_aleatoire_poly_min(j,List_fact,s,CCestimt,t);
				print"f=",f;
				g:=prend_poly_min_app_rend_exact(f,s,Zt,t);
				print"g=",g;
				test_irred:=IsIrreducible(g);
				print"test_irred=",test_irred;
				c:=c+1;
				end while;
                                if (test_irred eq true) then return true,g,L;
                                                   else return false,g,L;
			        end if;
end if;

end function;

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

fabriq_vandermonde:=function(L,CCestim)
/* Prend une liste L=[x1,x2,...,x_s] 
rend la matrice deVandermonde dont la premiere ligne est egale a 
1 x1 x1^2 ... x1^(s-1)*/

coeff:=[];
s:=#L;
for i:=1 to s do
	for j:=0 to s-1 do
 	coeff:=coeff cat [L[i]^j];
	end for;
end for;

V:=Matrix(CCestim,s,s,coeff);
return V;
end function;

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

matrice_pour_coeff_app:=function(L,f,CCestim)
/*prend la liste L:=[x1,x2,...x_s] et rend la matrice M qui est le produit de l'inverse de Vandermonde par la matrice diagonale avec des f'(x_i) sur la diagonale*/

V:=fabriq_vandermonde(L,CCestim);
VV:=V^(-1);
list:=[];
g:=Derivative(f);
s:=#L;
for i:=1 to s do
	list:=list cat [Evaluate(g,L[i])];
end for;
D:=DiagonalMatrix(CCestim,s,list);
M:=VV*D;
return M;
end function;

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

trouve_coeff_exact:=function(M,s,v)
/* M est la matrice donnée par matrice_pour_coeff_app, v est une matrice s lignes 1 colonne (M est s*s),
rend la matrice 1 ligne s colonnes  M*v, ou on a pris pour chaque coeff l'entier le plus proche*/

w:=M*v;
list:=[];
for i:=1 to s do list:=list cat [Round(Re(w[i,1]))]; end for;
w:=Matrix(IntegerRing(),s,1,list);
return w;
end function;
/* ********************************************************************* */

fabriq_fact_exact_mod_x:=function(P,n,m,List_fact,f,L,CCestim,K,alpha,Ky,y,Qxy)
/* prend la liste des facteurs,
 un polynome minimal de l'etension,
 le corps K engendré par f,
 alpha est un element primitif de K de polynome minimal f, 
et y est la variable que l'on va utilisee pour fabriquer le facteur exact,
rend un booleen qui nous dit si la reconnaissance a bien fonctionne, 
le facteur exact P1 modulo x dans Ky, et le facteur Q=P2*....*Ps modulo x Q est donc un polynome de Ky
*/

  print"Je simplifie l'extension obtenue";
temps_simpl:=Cputime();
nouvK<beta>,morph:=OptimizedRepresentation(K);
temps_simpl:=Cputime()-temps_simpl;
print"Temps (en secondes) pour la simplification:",temps_simpl;
print"La simplification est faite";
/*nouvKx<x>:=PowerSeriesRing(nouvK,n+1);*/
nouvKy<y>:=PolynomialRing(nouvK);

P1:=nouvKy!0;
g:=Derivative(f);
denom:=Evaluate(g,alpha);
s:=#List_fact;
M:=matrice_pour_coeff_app(L,f,CCestim);
for i:=0 to m do
	list:=[];
	for j:=1 to s do
		list:=list cat [Coefficient(UnivariatePolynomial(Coefficient(List_fact[j],1,0)),i)];
	end for;

	/* on vient de fabriquer la liste des coefficients de y^i avec 			tous ces conjuges*/

	v:=Matrix(CCestim,s,1,list);
	coeffz:=trouve_coeff_exact(M,s,v);
	temp:=K!0;	
	for k:=0 to s-1 do
		temp2:=K!coeffz[k+1,1];
		temp:=temp+(temp2)*(alpha^k)/denom;
	end for;
	
	/*temp est l'expression exacte des coeff y^i de P1*/
	
P1:=P1+morph(temp)*y^i;
end for;


F:=nouvKy!0;
Ptemp:=Qxy!P;
/*print"Parent(Coefficient(UnivariatePolynomial(Coefficient(P,1,0)),i)",Parent(Coefficient(UnivariatePolynomial(Coefficient(P,1,0)),1));*/
for i:=0 to n do
temp:=nouvKy!Coefficient(UnivariatePolynomial(Coefficient(Ptemp,1,0)),i);
F:=F+temp*y^i;
end for;
Q,R:=Quotrem(F,P1);

bool:=true;
for i:=0 to Degree(R) do
	if (Coefficient(R,i) ne 0) then bool:=false; end if;
end for;
/*print"R=",R;*/
return bool,nouvK,P1,Q;

end function;

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

trouve_rang_avec_vect_de_norme_inf_a_M:=function(G,dim,M)
/* prend la matrice de Gram G, son nombre deligne dimX et une constante M
rend le nombre de lignes consecutives en partant de la premiere qui sont remplies avec des vecteurs de normes^2 inferieur a M*/
  test:=true;
r:=0;
nbcol:=dim+1;
while (test eq true) and (r lt dim) do
   r:=r+1;
   n:=0;
   for i:=1 to nbcol do 
       n:=n+G[r,i]^2;
   end for;
   if (n gt M)  then test:=false; end if;
end while;
r:=r-1;
return r;
end function;

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

recup_bi_de_la_1_ligne:=function(M,dimX,X)
/* M est la RREFcomplete, X est la matrice qui a servi a obtenir le premier LLL,
rend la liste eds bi correspondant a la 1 ligne de M*/
  L:=[];
for i:=1 to NumberOfColumns(M) do
  if (M[1,i] eq 1) then L:=L cat[X[i,dimX+1]]; end if;
 end for;
 return L;
 end function;


/* ******************************************************************** */
fabriq_mat_pour_LLL_opt:=function(LBopt)
/* prend la liste LBopt et fabriq la matrice pour LLL avec en derniere colonne la liste LBopt*/

temp:=#LBopt;
Xopt:=RMatrixSpace(IntegerRing(),temp,temp+1)!0;
for i:=1 to temp do Xopt[i,i]:=1;
 Xopt[i,temp+1]:=LBopt[i];end for;
 return Xopt;
 end function;

/* ********************************************************************** */
fabriq_Sylvester:=function(f1,f2,K)
/* prend deux polynomes en 1 variable y a coeff dans K 
rend la matrice de Sylvester associe*/
degf1:=Degree(f1);
degf2:=Degree(f2);
dim:=degf1+degf2;
M:=RMatrixSpace(K,dim,dim)!0;

for i:=1 to degf2 do
   for j:=1 to degf1+1 do
    M[j+i-1,i]:=Coefficient(f1,degf1-j+1);
   end for;
end for;

for i:=degf2+1 to dim do
   for j:=1 to degf2+1 do
    M[j+i-1-degf2,i]:=Coefficient(f2,degf2-j+1);
   end for;
end for;

return M,dim,degf1,degf2;

  end function;

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

mon_bezout:=function(M,dim,degf1,degf2,Ky,y,K)
/* prend la matrice de Sylvester de f1 et f2 ,
dim est la dimension de la matrice
et rend les polynomes a et b de Ky tels que 1=a*f1+b*f2;*/
  temp:=UnivariatePolynomial(y);
  a:=Ky!0;
b:=Ky!0;
  L:=[];
for i:=1 to dim-1 do
  L:=L cat [K!0];
end for;
L:=L cat [K!1];
v:=Vector(dim,L);
w:=Solution(Transpose(M),v);
for i:=1 to degf2 do
a:=a+w[i]*temp^(degf2-i);
end for;
for i:=degf2+1 to dim do
b:=b+w[i]*temp^(degf1-(i-degf2));
end for;
return a,b;
end function;


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

ma_remontee_de_hensel:=function(P,P1,Q,m,n,Kxy,x,y,Ky,K)
/* prend le polynome P afactoriser P appartient a Kxy,
P1 et Q sont deux facteurs exact mod x ils appartiennent a Ky, 
m est le degre total des facteurs,
rend deux polynomes f1 et f2 dans Kxy qui sont les releves de P1 et de Q*/
  print"Je calcule la relation de Bezout entre P1 et Q";
/*temps_bezout_defaut:=Cputime();
  c,a,b:=XGCD(P1,Q);
temps_bezout_defaut:=Cputime()-temps_bezout_defaut;
print"Temps (en secondes) pour obtenir Bezout avec version Magma:",temps_bezout_defaut;*/


temps_bezout:=Cputime();
M,dim,degf1,degf2:=fabriq_Sylvester(P1,Q,K);

a,b:=mon_bezout(M,dim,degf1,degf2,Ky,y,K);
temps_bezout:=Cputime()-temps_bezout;
print"Temps (en secondes) de mon Bezout:",temps_bezout;

/*temps_bezout:=Cputime()-temps_bezout;
print"Temps (en secondes) pour obtenir un relation de Bezout:",temps_bezout;*/
/*print"c=",c;
print"degre a",Degree(a);
print"degre b",Degree(b);*/
F1:=MultivariatePolynomial(Kxy,P1,y);
F2:=MultivariatePolynomial(Kxy,Q,y);
/*temp:=UnivariatePolynomial(y^(m+1));*/
for k:=1 to m do
  if ( (k mod 10) eq 0) then  print"j'effectue la remontee mod x^",k; end if;
c:=0;
for j:=0 to k do
c:=c+Coefficient(F1,x,j)*Coefficient(F2,x,k-j);
end for;
C:=Coefficient(P,x,k)-c;



/*q:=ExactQuotient(P-F1*F2,x^k);*/
C:=UnivariatePolynomial(C);
/*f1:=UnivariatePolynomial(Coefficient(F1,x,0));
f2:=UnivariatePolynomial(Coefficient(F2,x,0));*/
/*r2:=a*C;r2:=r2 mod Q;*/
r1:=b*C;r1:=r1 mod P1;
r2:=ExactQuotient(C-Q*r1,P1);
r2:=MultivariatePolynomial(Kxy,r2,y);
r1:=MultivariatePolynomial(Kxy,r1,y);
F1:=F1+x^k*r1;
F2:=F2+x^k*r2;
end for;
F1bon:=F1;

for k:=m+1 to n-m do
  if ( (k mod 10) eq 0) then   print"j'effectue la remontee mod x^",k;end if;
c:=0;
      for j:=0 to m do
c:=c+Coefficient(F1bon,x,j)*Coefficient(F2,x,k-j);
end for;
C:=Coefficient(P,x,k)-c;


/*q:=ExactQuotient(P-F1*F2,x^k);*/


C:=UnivariatePolynomial(C);
C:=MultivariatePolynomial(Kxy,ExactQuotient(C,P1),y);
/*f1:=UnivariatePolynomial(Evaluate(F1,1,0));*/
   /*f2:=UnivariatePolynomial(Coefficient(F2,x,0)); 
   r2:=a*C;r2:=r2 mod f2;*/
/*r1:=b*C;r1:=r1 mod f1;*/
    /*r2:=MultivariatePolynomial(Kxy,r2,y);*/
/*r1:=MultivariatePolynomial(Kxy,r1,y);*/
/*F1:=F1+x^k*r1;*/
F2:=F2+x^k*C;
end for;
return F1bon,F2;


  end function; 

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

maverif:=function(P1,Q,P,n)
/*prend P1, Q et P verifie que P1*Q=P;*/


Pprod:=P1*Q;

montest:= true;
for i:=0 to n do
    for j:=0 to n-i do
    coefP:=Coefficient(Coefficient(P,2,i),1,j);
coefPprod:=Coefficient(Coefficient(Pprod,2,i),1,j);
if (coefP ne coefPprod) then 
  print"coefP",coefP;
print"coefpprod",coefPprod;
montest:=false; end if;
end for;
end for;
return montest;
end function;

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

intrinsic Abs_Fac_Knap(p,prec) -> .
{
P must be an irreducible polynomial in Z[x,y] monic in y.
P must have no linear factor.
prec must be 1
Abs_Fac_Knap returns
K the number field generated by all the coefficients of one absolute factor P_i
P1 one absolute factor
and P.

}




  temps_global:=Cputime();
Zxy<x,y>:=PolynomialRing(IntegerRing(),2);
Qxy<x,y>:=PolynomialRing(RationalField(),2);
P:=Zxy!p;
n:=TotalDegree(p);

Z:=IntegerRing();
Zt<t>:=PolynomialRing(Z);
CCestim<i>:=ComplexField();
CCestimt<t>:=PolynomialRing(CCestim);
CCestimy<y>:=PolynomialRing(CCestim);
CCestimxy<x,y>:=PolynomialRing(CCestim,2);




if (prec eq 1) then
lambda:=0;/*Random(-10,10);*/
print "lambda=",lambda;
x0:=1;/*Random(-10,10);*/
P:=chgvar(P,x0,lambda,x,y);
print "x0=",x0;
Px0:=choixfibre(P);
Px0:=CCestimy!Px0;
/*print"Px0=",Px0;*/
C,eta:=estimprec(n);
/*print"C=",C;*/
if (n gt 200) then eta:=eta+50; end if;
/*print "eta=",eta;*/
else 
C:=10^(prec); 
eta:=prec+20;
if (n gt 200) then eta:=eta+50; end if;
print "eta=",eta;
P:=Zxy!p;
lambda:=Random(-10,10);
print "lambda=",lambda;
x0:=Random(-10,10);
P:=Zxy!chgvar(P,x0,lambda,x,y);
print "x0=",x0;
Px0:=choixfibre(P);
Px0:=CCestimy!Px0;
end if;

print"Patience je calcule les b_i....";
temps_bi:=Cputime();
L:=estimracine(Px0,eta);
A,B,alpha0,beta0,gamma0,delta0,epsilon0:=estim_ai_bi(P,CCestimy,n,L);
print"Temps (en secondes) pour obtenir les b_i",Cputime()-temps_bi;
/*print"Liste de b_i:",B;*/
test_fact_lin:=facteur_lineaire(B,n,eta);
if (test_fact_lin eq true) then
    print"Il y a de grandes chances que ce polynome n'aient que des facteurs lineaires, essayez une autre methode.";
 return RationalField(),P,P;
    end if;
l,LR,AR,BR:=estim_l_et_pretri(L,A,B,n,eta);
print "nombre de bi réels=",l;
Ly,La,Lb,Lbconj:=tribi(L,LR,A,AR,B,BR,n,eta);

temps_LLL:=Cputime();
X:=fabriq_mat_pour_LLL(Lb,n,l,C);
dimX:=NumberOfRows(X);

F:=X*Transpose(X);
print"Je commence LLL";
G1,T:=LLLGram(F:Sort:=false,Method:="Integral");
print"J'ai fini LLL";
/*print"Parent(G1)",Parent(G1);*/
LX:=T*X;
/*print"Parent(LX)",Parent(LX);*/
/*print"Lx[1]",LX[1];
print"LX[2]",LX[2];*/
print"Je calcule GramSchmidt";
temps_Gram:=Cputime();
G2:=nouv_nouv_GramSchmidt(LX,G1);
temps_Gram:=Cputime()-temps_Gram;
print"J'ai fini de calculer GramSchmidt, cela a pris (en secondes):",temps_Gram;
/*print"Parent(G2)",Parent(G2);*/

/*print"Je commence LLL";
LX:=LLL(X:Sort:=false,Method:="Integral");
print"J'ai fini LLL";
print"Je calcule GramSchmidt";
temps_Gram:=Cputime();
G:=nouv_GramSchmidt(LX,dimX,dimX+1);
temps_Gram:=Cputime()-temps_Gram;
print"J'ai fini de calculer GramSchmidt, cela a pris (en secondes):",temps_Gram;*/

M:=dimX+((1.5)*dimX)^2;
print "M^2=",M;
k0:=Norm_vec_Gram_rend_k0(G2,dimX,dimX+1,M);
print "dim=",dimX;
print "k0=", k0;

RREF:=EchelonForm(ColumnSubmatrix(RowSubmatrix(LX,1,k0),dimX));
if (k0 lt dimX) then print "RREF des k0 premiers vecteurs=",RREF;end if;
/* Ici commence l'optimization*/ 

k1:=trouve_rang_avec_vect_de_norme_inf_a_M(G2,dimX,M);
print"k1=",k1;
if (k1 lt k0) and (k1 gt 1) then 
    print"J'utilise l'optimization car k1=",k1;
RREFcomplete:=EchelonForm(ColumnSubmatrix(RowSubmatrix(LX,1,k1),dimX));
print"La RREF apres optimization est:",RREFcomplete;
    LBopt:=recup_bi_de_la_1_ligne(RREFcomplete,dimX,X);
    Xopt:=fabriq_mat_pour_LLL_opt(LBopt);
Fopt:=Xopt*Transpose(Xopt);
G1opt,Topt:=LLLGram(Fopt:Sort:=false,Method:="Integral");
LXopt:=Topt*Xopt;
G2opt:=nouv_nouv_GramSchmidt(LXopt,G1opt);
    dimXopt:=NumberOfRows(Xopt);
    Mopt:=dimXopt+((1.5)*dimXopt)^2;
    k2:=Norm_vec_Gram_rend_k0(G2opt,dimXopt,dimXopt+1,Mopt);
    if (k2 eq 1) then 
         k0:=k1;
         RREF:=EchelonForm(ColumnSubmatrix(RowSubmatrix(LX,1,k0),dimX)); end if;
end if;
if (k1 eq dimX) and (k0 eq dimX) then return $$(p,2*eta); end if;
if (k1 eq 0) and (k0 eq dimX) then return $$(p,2*eta); end if;

/* Ici se termine l'optimization*/



test01:=Form_Echelon_que_des_0_1(RREF,k0,dimX);
print "test01=", test01;



/*if (k0 eq dimX) then
return $$(p,2*eta); end if;*/ 
/*on suppose qu'il nya pas de facteurs de degre 1 donc le reseau ne pas etre de dimension = dimX  et coherent donc si on obtient ce cas c'est que nous n'avons pas éliminer de vecteurs
*/



while (test01 eq false) or (k0 eq dimX) do
  if (test_situation_generique_pour_les_relations(RREF) eq false) then 
    print "La situation semblait non generique";
return $$(p,2*eta);
else
  print"J'esaie de faire diminuer la dimension du reseau";
print"eta=",2*eta;
  eta,CCestim,CCestimy,CCestimxy,l,Ly,La,Lb,Lbconj,nouvLX,nouvk0:=recommence_augmente_prec(P,Px0,n,C^2,2*eta,LX,k0);
LX:=nouvLX;
k0:=nouvk0;
print "k0=",k0;
RREF:=EchelonForm(ColumnSubmatrix(RowSubmatrix(LX,1,k0),dimX));
test01:=Form_Echelon_que_des_0_1(RREF,k0,dimX);
print "test01=",test01;
end if;
end while;

/* a ce niveau on a une RREF qui n'a que des 0-1 vecteurs */


bool,m,L1,L2,L3:=coherence_RREF_rend_liste_0_1_vecteur_selon_type(RREF,l,n,C,Lb);
print "la coherence de la RREF est=",bool;
print "les polynomes que l'on recherche sont de degre=",m;


if (m eq n) and (l ne 0) then 
temps_global:=Cputime()-temps_global;   
print"P est absolument irreductible";  
print"Temps total (en secondes):",temps_global;
return RationalField(),P,P; end if; 


if (bool eq false) then return $$(p,2*eta);
		   else List_fact:=[];
		        List_fact:=List_fact cat list_fact_reel(L1,l,n,Ly,La,Lb,CCestimxy,x,y);
			print"list_fact_reel est ok";
			print "nb de facteurs obtenus=",#List_fact;
			List_fact:=List_fact cat list_fact_complex_bon_deg(L2,l,n,Ly,La,Lb,CCestimxy,x,y);
			print"list_fact_complex_bon_deg est ok";
			print"nb de facteurs obtenus=",#List_fact;
testbool,fact:=list_fact_complex_apres_redecomp(P,n,Px0,C,eta,L3,l,m,Ly,La,Lb,alpha0,beta0,gamma0,delta0,epsilon0,CCestimxy,x,y);
if (testbool eq false) then return $$(p,2*eta);
				else List_fact:=List_fact cat fact; 
				print"list_fact_complex_apres_redecomp est ok";
				print"nb de facteurs obtenus=",#List_fact;
				end if;

end if;
temps_LLL:=Cputime()-temps_LLL;
print"Temps (en secondes) pour obtenir les facteurs approchés mod x^3:",temps_LLL;
s:=#List_fact;



if (s eq 1)  then 
temps_global:=Cputime()-temps_global;   
print"P est absolument irreductible";  
print"Temps total (en secondes):",temps_global;
return RationalField(),P,P; end if; 

test_poly_min_irred,f,L:=trouve_poly_min(List_fact,m,s,Zt,CCestimt,t);
print"f avant simplification=",f;
K<alpha>:=NumberField(f);
/*Kx<x>:=PowerSeriesRing(K,n+1);*/
Ky<y>:=PolynomialRing(K);
test_reconnaissance,K,P1,Q:=fabriq_fact_exact_mod_x(P,n,m,List_fact,f,L,CCestim,K,alpha,Ky,y,Qxy);

print"test_reconnaissance=",test_reconnaissance;

if (test_reconnaissance eq false) then print"On arrete test_reconnaissance est faux";
return K,P1,P;
end if;

print"K=",K;
P:=Qxy!P;
Kxy<x,y>:=PolynomialRing(K,2);
P:=Kxy!P;
print"J'effectue la remontée de Hensel";
temps_hensel:=Cputime();
P1,Q:=ma_remontee_de_hensel(P,P1,Q,m,n,Kxy,x,y,Ky,K);
temps_hensel:=Cputime()-temps_hensel;
print"Temps (en secondes) pour la remontee de Hensel:",temps_hensel; 


print"Je verifie le resultat";
print"Je fais ma verification";
temps_ma_verif:=Cputime();
montest:=maverif(P1,Q,P,n);
temps_ma_verif:=Cputime()-temps_ma_verif;
print"mon_test=",montest;
print"Ma verification prend:",temps_ma_verif;





print"J'effectue le dernier changement de variable";
temps_chvar:=Cputime();
x0retour:=K!-x0;
lambdaretour:=K!-lambda;
P1:=chgvar(P1,x0retour,lambdaretour,x,y);
temps_chvar:=Cputime()-temps_chvar;
print"Temps (en secondes) pour effectuer le changement de variable:",temps_chvar;


/*print"Je fais IsdivisibleBy";
temps_verif:=Cputime();
test_final:=IsDivisibleBy(Kxy!p,P1);
temps_verif:=Cputime()-temps_verif;
print"Temps (en secondes) pour la verification avec IsDivisibleBy:",temps_verif;
print"test_final=",test_final;*/



print"Temps (en secondes) pour obtenir les facteurs approchés mod x^3:",temps_LLL;
temps_global:=Cputime()-temps_global;
print"Temps total (en secondes):",temps_global;
return K,P1,P;

end intrinsic;

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


