base_standard(rank)=
{
   standard_basis = vector(rank);
   for(i=1,rank,standard_basis[i] = vector(rank));
   j=1;
   for(i=1,rank,standard_basis[i][j] = 1;j = j++;);
   return(standard_basis);
}

addhelp(base_standard,"base_standard(rank): rank being an integer larger than or equal to one, gives the standard basis for the vector space of dimension rank.")

pairing(matrice,vecteur1,vecteur2)=
{
   return(vecteur1 * matrice * mattranspose(vecteur2));
}

addhelp(pairing,"pairing(matrice,vecteur1,vecteur2):  gives the pairing of the two vectors vecteur1 and vecteur2 for the symmetric bilinear form defined by the symmetric matrix matrice.")


projection(matrice,vecteur1,vecteur2)=
{
   return((pairing(matrice,vecteur1,vecteur2)/pairing(matrice,vecteur2,vecteur2))*vecteur2);

}


addhelp(projection,"projection(matrice,vecteur1,vecteur2):  gives the projection of vecteur1 onto vecteur2 with respect to the pairing given by the symmetric matrix matrice.")


orthogonal_basis(matrice)=
{
   ort = vector(length(matrice[1,]));

   standard_basis = vector(length(ort));
   for(i=1,length(ort),standard_basis[i] = vector(length(ort)));
   j=1;
   for(i=1,length(ort),standard_basis[i][j] = 1;j = j++;);

   ort[1] = standard_basis[1];
   for(i=2,length(ort),
      bidon = standard_basis[i];
      for(j=1,i-1,
	 bidon = bidon - projection(matrice,standard_basis[i],ort[j]);
      );
      ort[i] = bidon;
   );

   return(ort);
}


addhelp(orthogonal_basis,"orthogonal_basis(matrice):  gives an orthogonal basis for the orthogonal geometry defined by the symmetric matrix matrice.")

/* The following function is to replace the previous orthogonal_basis function in the case of D_{8}^{++}, since something funny is happening for this case when doing the Gram-Schmidt process. */

orthogonal_basis1(matrice)=
{
   ort = vector(length(matrice[1,]));

   standard_basis = vector(length(ort));
   for(i=1,length(ort),standard_basis[i] = vector(length(ort)));
   j=1;
   for(i=1,length(ort),standard_basis[i][j] = 1;j = j++;);

   ort[1] = standard_basis[1];
   for(i=2,length(ort),
      if(i==9,ort[i] = [4, 8, 6, 12, 10, 8, 6, 4, 1, 1],
      bidon = standard_basis[i];
      for(j=1,i-1,
	 bidon = bidon - projection(matrice,standard_basis[i],ort[j]);
      );
      ort[i] = bidon;
      );
   );

   return(ort);
}






apply_aut(auto,vecteur)=
{
   resultat1 = vector(length(vecteur));

   for(i=1,length(vecteur),
      resultat1 = resultat1 + vecteur[i] * auto[i];
   );

   return(resultat1);
}


addhelp(apply_aut,"apply_aut(auto,vecteur):  applies the automorphism auto to the vector vecteur.  Here, an automorphism is given by a list of vectors one for each element of the standard basis.")

reflection(matrice,vecteur)=
{
   ort = vecteur;
   standard_basis = vector(length(ort));
   for(i=1,length(ort),standard_basis[i] = vector(length(ort)));
   j=1;
   for(i=1,length(ort),standard_basis[i][j] = 1;j = j++;);

   resultat = vector(length(vecteur));
   for(i=1,length(resultat),
      resultat[i] = standard_basis[i] - (2*pairing(matrice,standard_basis[i],vecteur)/pairing(matrice,vecteur,vecteur))*vecteur;
   );

   return(resultat);

}

addhelp(reflection,"reflection(matrice,vecteur):  gives the hyperplane reflection corresponding to the vector vecteur with respect to the orthogonal geometry given by the matrix matrice.  The result is a list of vectors which defined the automorphism by specifying where the elements of the standard basis are mapped to.")




composition(autom1,autom2)=
{
   ort = autom1;
   standard_basis = vector(length(ort));
   for(i=1,length(ort),standard_basis[i] = vector(length(ort)));
   j=1;
   for(i=1,length(ort),standard_basis[i][j] = 1;j = j++;);

   resultat2 = vector(length(ort));
   for(i=1,length(ort),
      resultat2[i] = apply_aut(autom1,apply_aut(autom2,standard_basis[i]));
   );

   return(resultat2);
}


addhelp(composition,"composition(autom1,autom2):  gives the composition autom1 \circ autom2 of the two automorphisms.")


product_of_reflection(matrice,auto2)=
{
   base_orthogonal = orthogonal_basis(matrice);

   ort = auto2;
   standard_basis = vector(length(ort));
   for(i=1,length(ort),standard_basis[i] = vector(length(ort)));
   j=1;
   for(i=1,length(ort),standard_basis[i][j] = 1;j = j++;);


   record_reflections = vector(length(ort));

   auto_swap = standard_basis;

   for(k=1,length(base_orthogonal),
      auto2 = composition(auto_swap,auto2);
      if(pairing(matrice,apply_aut(auto2,base_orthogonal[k]) - base_orthogonal[k],apply_aut(auto2,base_orthogonal[k])-base_orthogonal[k]) != 0,
	 record_reflections[k] = apply_aut(auto2,base_orthogonal[k]) - base_orthogonal[k];
	 auto_swap = reflection(matrice,apply_aut(auto2,base_orthogonal[k]) - base_orthogonal[k]);,
	 record_reflections[k] = [base_orthogonal[k],apply_aut(auto2,base_orthogonal[k]) + base_orthogonal[k]];
	 auto_swap = composition(reflection(matrice,base_orthogonal[k]),reflection(matrice,apply_aut(auto2,base_orthogonal[k]) + base_orthogonal[k]));  
	); 
   );

   count = 0;
   for(i=1,length(record_reflections),if(length(record_reflections[i])==2,;count = count+2,count = count++));

   res = vector(count);
   k=1;
   for(i=1,length(record_reflections),if(length(record_reflections[i])==2,res[k] = record_reflections[i][1];res[k+1] = record_reflections[i][2];k = k+ 2,res[k] = record_reflections[i];k = k++ ));

   return(res);
}

addhelp(product_of_reflection,"product_of_reflection(matrice,auto2): given the automorphism auto2, returns a list of vectors.  The automorphism can be written as a product of hyperplane reflections corresponding to those vectors in the list.  auto2 has to be an orthogonal map, even though the code does not do this sanity check.  It auto2 is not orthogonal, it will probably return garbage.")



spinor_norm(matrice,vecteur_long)=
{
   hip = 1;
   for(i=1,length(vecteur_long),hip = hip * pairing(matrice,vecteur_long[i],vecteur_long[i]));

   return(hip);
}


addhelp(spinor_norm,"spinor_norm(matrice,vecteur_long):  given a list of vectors, as returned by the function product_of_reflection for instance, returns the products of quadratic form applied to each vector in the list.  Hence, it gives the spinor norm (up to a rational square number).")