w1 = quadgen(-4,'w1);
w2 = quadgen(-8,'w2);
w3 = quadgen(-3,'w3);
w7 = quadgen(-7,'w7);
w11 = quadgen(-11,'w11);

/* Hermitian scalar product */
h(u,v) = {
	if(matsize(u)[2] == 1,
		u~*conj(v),
		u*conj(v)~
  )
};


M={[22 + 49*w11 , 58 + 87*w11 , 86 + 22*w11 ,  6 + 62*w11 , 91 + 91*w11 , 83 + 47*w11 , 34 + 13*w11 , 79 + 95*w11 , 16 + 68*w11 , 97 + 88*w11;
    59 + 84*w11 , 53 + 96*w11 ,  34 + 2*w11 , 44 + 36*w11 ,  56 + 6*w11 , 50 + 68*w11 , 28 + 96*w11 , 47 + 10*w11 , 76 + 81*w11 ,  39 + 2*w11;
    58 + 30*w11 ,  42 + 6*w11 , 74 + 97*w11 ,  61 + 5*w11 , 40 + 32*w11 ,  3 + 12*w11 ,  96 + 2*w11 , 98 + 21*w11 , 22 + 61*w11 , 80 + 75*w11;
    24 + 98*w11 , 51 + 45*w11 , 80 + 36*w11 ,  7 + 82*w11 , 17 + 75*w11 , 99 + 24*w11 , 40 + 69*w11 , 88 + 77*w11 , 12 + 91*w11 , 82 + 81*w11;
    72 + 43*w11 ,  4 + 37*w11 , 73 + 46*w11 , 98 + 96*w11 , 43 + 70*w11 , 81 + 77*w11 ,  98 + 4*w11 , 47 + 78*w11 ,  9 + 73*w11 , 46 + 15*w11;
    84 + 26*w11 ,  2 + 31*w11 , 54 + 46*w11 , 38 + 36*w11 , 99          , 84 + 45*w11 , 43 + 70*w11 , 82 + 38*w11 ,  7 + 50*w11 , 91 + 74*w11;
    64 + 37*w11 , 78 + 13*w11 , 47 + 91*w11 , 85 + 49*w11 , 20 + 33*w11 , 82 + 49*w11 , 33 + 51*w11 , 18 + 77*w11 , 86 + 79*w11 , 20 + 90*w11;
    39 + 43*w11 , 17 + 39*w11 , 26 + 33*w11 , 69 + 66*w11 , 19 + 79*w11 , 58 + 95*w11 ,  58 + 2*w11 , 71 + 26*w11 , 53 + 33*w11 , 11 + 34*w11;
    72 + 86*w11 , 51 + 49*w11 ,  1 + 21*w11 , 83 + 45*w11 , 95 + 82*w11 , 64 + 98*w11 , 23 + 24*w11 , 23 + 54*w11 , 41 + 57*w11 , 93 + 57*w11;
    94 + 89*w11 , 78 + 44*w11 , 40 + 94*w11 , 93 + 35*w11 , 78 + 93*w11 , 35 + 42*w11 , 37 + 97*w11 , 13 + 31*w11 ,  96 + 2*w11 , 54 + 11*w11]};

M_2 = {[-5 + 2*w1 ,  5 - 2*w1 ,  -2;
w1  ,  -3 - w1  ,  -1 - 4*w1;
3*w1  ,  -3 + 4*w1  ,  1 + 5*w1]};

M_2J = {[-5 + 2*w3 ,  5 - 2*w3   ,  -2;
        w3        ,  -3 - w3    ,  -1 - 4*w3;
        3*w3      ,  -3 + 4*w3  ,  1 + 5*w3]};



mapgen = Map([-4 , w1 ; -8 , w2 ; -3 , w3 ; -7 , w7 ; -11 , w11]);

/* Take u in Q[w] and return the closest element in the integer ring Z[w] with respect to the algebraïc norm */
round_QF(u) = {
  if(u.disc % 4 == 0, 
    round(real(u)) + round(imag(u))*mapget(mapgen,u.disc),

    /* else : On se ramène au cas du parallèlogramme fondamentale */
    my(w = mapget(mapgen,u.disc));
    my(a = floor(real(u)) + floor(imag(u))*w);
    my(b = u - a);

    my(x = real(b) + imag(b)/2); /* b = x + iy */

    /* On différencie les cas */
    if(x <= 0.5,
      if((-b.disc + 1)*2*imag(b) <= -b.disc + 1 - 4*real(b),
        a,
        a + w),
      if(x <= 1,
        if((-b.disc - 1)*2*imag(b) <= -b.disc - 3 + 4*real(b),
          a + 1,
          a + w),
        if((-b.disc + 1)*2*imag(b) <= -b.disc + 5 - 4*real(b),
          a + 1,
          a + 1 + w
        )
      )
    )
  )
}




gso_norm(M) = {
	my(N = vector(matsize(M)[2],i,h(M[,i],M[,i])));
	my(mu = matid(matsize(M)[1]));
	for(j=1,matsize(M)[2],
		for(i=j+1,matsize(M)[2],
			mu[i,j] = (1/N[j])*(h(M[,i],M[,j]) - sum(k=1,j-1,conj(mu[j,k])*mu[i,k]*N[k]));
			N[i] -= norm(mu[i,j])*N[j];
		)
	);
	[N,mu]
};


/* i>=2 */
lll_swap(M,N,mu,i) = {
	my([M1,N1,mu1] = [M,N,mu]);
	
	M1[,i] = M[,i-1];
	M1[,i-1] = M[,i];

	N1[i-1] = N[i] + norm(mu[i,i-1])*N[i-1];
	N1[i] = N[i]*N[i-1]/N1[i-1];

	mu1[i,i-1] = conj(mu[i,i-1])*N[i-1]/N1[i-1];

	for(k=1,i-2,
		mu1[i,k] = mu[i-1,k];
		mu1[i-1,k] = mu[i,k];
	);

	for(k=i+1,matsize(M)[2],
		mu1[k,i-1] = mu[k,i]*N[i]/N1[i-1] + mu[k,i-1]*mu1[i,i-1];
		mu1[k,i] = mu[k,i-1] - mu[k,i]*mu[i,i-1];
	);
	[M1,N1,mu1]
};

/**
	* M est la base
	* d est le delta (proche de 1, exemple 0.99)
	* l est mk (0.5 dans Z)
  **/
lll_QF_col(M,d,l) = {
	
	my([N,mu] = gso_norm(M));
	my(M1 = M);
	i=2;
	
	while(i <= matsize(M)[2],	
		forstep(j=i-1,1,-1,
			if(norm(mu[i,j])>l,
				/*size reduction*/
				c = round_QF(mu[i,j]);
				M1[,i] -= c*M1[,j];
				for(k=1,j,
					mu[i,k] -= c*mu[j,k];
				);
			)
		);
		if(i>1 && N[i] < (d - norm(mu[i,i-1]))*N[i-1],
			[M1,N,mu] = lll_swap(M1,N,mu,i);
			i--,
			i++
		)
	);
	M1
};


islllreduced(M,d,l) = {
  my([N,mu] = gso_norm(M));
  a = sum(i = 1 , matsize(M)[2] , sum(j = 1, i-1 , norm(mu[i,j]) > l));
  b = sum(i = 2 , matsize(M)[2] , N[i] < (d - norm(mu[i,i-1]))*N[i-1]);
  a+b == 0;
};






gram_mat_row(M) = {  
  matrix(matsize(M)[1],matsize(M)[1],i,j,h(M[i,],M[j,]))
}


gso_norm_row(M) = {
  my(N = vector(matsize(M)[1],i,real(h(M[i,],M[i,]))));
  my(mu = matid(matsize(M)[1]));
  for(j=1,matsize(M)[1],
    for(i=j+1,matsize(M)[1],
      mu[i,j] = (1/N[j])*(h(M[i,],M[j,]) - sum(k=1,j-1,conj(mu[j,k])*mu[i,k]*N[k]));
      N[i] -= real(norm(mu[i,j])*N[j]);
    )
  );
  [N,mu]
};

/* i>=2 */
lll_swap_row(M,N,mu,i) = {
  my([M1,N1,mu1] = [M,N,mu]);
  
  M1[i,] = M[i-1,];
  M1[i-1,] = M[i,];

  N1[i-1] = N[i] + norm(mu[i,i-1])*N[i-1];
  N1[i] = N[i]*N[i-1]/N1[i-1];

  mu1[i,i-1] = conj(mu[i,i-1])*N[i-1]/N1[i-1];

  for(k=1,i-2,
    mu1[i,k] = mu[i-1,k];
    mu1[i-1,k] = mu[i,k];
  );

  for(k=i+1,matsize(M)[1],
    mu1[k,i-1] = mu[k,i]*N[i]/N1[i-1] + mu[k,i-1]*mu1[i,i-1];
    mu1[k,i] = mu[k,i-1] - mu[k,i]*mu[i,i-1];
  );
  [M1,N1,mu1]
};

/**
  * M est la base
  * d est le delta (proche de 1, exemple 0.99)
  * l est mk (0.5 dans Z[i])
  **/
lll_QF_row(M,d=0.99,l=0.51) = {
  
  my([N,mu] = gso_norm_row(M));
  my(M1 = M);
  i=2;
  
  while(i <= matsize(M)[1], 
    forstep(j=i-1,1,-1,
      if(norm(mu[i,j])>l,
        /*size reduction*/
        c = round_QF(mu[i,j]);
        M1[i,] -= c*M1[j,];
        for(k=1,j,
          mu[i,k] -= c*mu[j,k];
        );
      )
    );
    if(i>1 && N[i] < (d - norm(mu[i,i-1]))*N[i-1],
      [M1,N,mu] = lll_swap_row(M1,N,mu,i);
      i--,
      i++
    )
  );
  M1
};


islllreduced_row(M,d=0.99,l=0.51) = {
  my([N,mu] = gso_norm_row(M));
  a = sum(i = 1 , matsize(M)[1] , sum(j = 1, i-1 , norm(mu[i,j]) > l));
  b = sum(i = 2 , matsize(M)[1] , N[i] < (d - norm(mu[i,i-1]))*N[i-1]);
  a+b == 0;
};





/* i>=2 */
lll_swap_row_u(M,U,N,mu,i) = {
  my([M1,U1,N1,mu1] = [M,U,N,mu]);
  
  M1[i,] = M[i-1,];
  M1[i-1,] = M[i,];

  U1[i,] = U[i-1,];
  U1[i-1,] = U[i,];

  N1[i-1] = N[i] + norm(mu[i,i-1])*N[i-1];
  N1[i] = N[i]*N[i-1]/N1[i-1];

  mu1[i,i-1] = conj(mu[i,i-1])*N[i-1]/N1[i-1];

  for(k=1,i-2,
    mu1[i,k] = mu[i-1,k];
    mu1[i-1,k] = mu[i,k];
  );

  for(k=i+1,matsize(M)[1],
    mu1[k,i-1] = mu[k,i]*N[i]/N1[i-1] + mu[k,i-1]*mu1[i,i-1];
    mu1[k,i] = mu[k,i-1] - mu[k,i]*mu[i,i-1];
  );
  [M1,U1,N1,mu1]
};


/**
  * M est la base
  * d est le delta (proche de 1, exemple 0.99)
  * l est mk (0.5 dans Z[i])
  * return the tranformation matrix
  **/
lll_QF_row_u(M,d=0.99,l=0.51) = {
  
  my(M1 = M);
  my([N,mu] = gso_norm_row(M1));
  my(U = matid(matsize(M1)[1]));

  i=2;
  
  while(i <= matsize(M1)[1], 
    forstep(j=i-1,1,-1,
      if(norm(mu[i,j])>l,
        /*size reduction*/
        c = round_QF(mu[i,j]);
        M1[i,] -= c*M1[j,];
        U[i,] -= c*U[j,];
        for(k=1,j,
          mu[i,k] -= c*mu[j,k];
        );
      )
    );
    if(i>1 && N[i] < (d - norm(mu[i,i-1]))*N[i-1],
      [M1,U,N,mu] = lll_swap_row_u(M1,U,N,mu,i);
      i--,
      i++
    )
  );
  [M1,U]
};