\\ ---------------  GP code  ---------------------------------------
\\
\\ Time-stamp: <Mon, Jan 13, 2003 - 14:09:59 - apacetti>
\\
\\ Description: Computes change of
\\ basis to achieve symplectic basis of a skew-symmetric matrix
\\
\\
\\ Original Author:     Fernando Rodriguez-Villegas 
\\                      villegas@math.utexas.edu
\\                      University of Texas at Austin
\\
\\
\\ With the assistance of:  Ariel Pacetti
\\			    apacetti@math.utexas.edu
\\			    University of Texas at Austin
\\
\\ Created:             Mon Jan 11 1999
\\
\\-----------------------------------------------------------------

\\ Given skew-symmetric matrix a find matrix m so that m~*a*m is in
\\ standard form (like J); in terms of symplectic basis re-writes a in
\\ a basis e_1, ..., e_n,f_1,...,f_n (n needs to be even)

symbas(a)=
      {local(ax,n,d,ux,m,w,e,e1,u,u1,jj,v);
      n=length(a);
      v=vector(n,k,matid(n)[k,]);
      forstep(j=1,n-3,2,
      u=vector(n-j+1,k,v[j]*a*v[k+j-1]~);
      e=content(u);jj=j;
      for(l=1,n-j,
        u1=vector(n-j+1,k,v[j+l]*a*v[k+j-1]~);
        e1=content(u1);
          if(e1<e,
            u=u1;e=e1;jj=j+l
            ,)
        );
          if(jj-j,
             u1=v[jj];v[jj]=v[j];v[j]=u1
             ,);
      u=vector(n-j+1,k,v[j]*a*v[k+j-1]~);
      ax=mbezout(u);
      d=ax[2];ax=ax[1];
      w=sum(k=1,n-j+1,ax[k]*v[k+j-1],0);
      ux=vector(n-j,k,w*a*v[k+j]~);
         for(k=j+1,n,
             v[k]=v[k]-u[k-j+1]/d*w+ux[k-j]/d*v[j];
             );
      m=matrix(n,n-j,s,t,v[t+j][s]);
      m=m*qflll(m,1);
      v[j+1]=w;
      for(k=0,length(m)-1,
               v[j+2+k]=Vec(vecextract(m,2^k))[1]~
        );
           );
      if(v[n-1]*a*v[n]~<0,v[n]=-v[n],);
      m=Mat([]);\
      forstep(k=1,n,2,m=concat(m,v[k]~));
      forstep(k=2,n,2,m=concat(m,v[k]~));
      m
}

\\======================================================================
\\ Computes the gcd of all the entries of the vector v, and gives the
\\ "smallest" linear combination.

mbezout(r)=
       {local(a,d,m,z,n);
       d=r[1];a=[1];n=length(r);
  if(n>1,
   if(norml2(r),
    for(k=2,n,
       z=bezout(d,r[k]);
       a=concat(z[1]*a,z[2]);
       d=z[3]
       );
     m=matkerint(Mat(r));
     sign(d)*[(m*round(matinverseimage(m,(d*r/(r*r~)-a)~)))~+a,d],
      [r,0]
      ),
     [[1],1]
    )
}

\\-----------------------------------------------------------------
\\ Computes the symmetric form, not the change of base, and the length
\\ of a need not be even. 

sym(a)=
	{local(n1,n2,m,aux);
	n1=length(a);
	n2=length(matkerint(a));
	m=[;];
	for(k=1,n2,m=concat(m,Mat(vector(n1,l,0)~)));
	aux=matsnf(a);
	forstep(k=1,(n1-n2),2,
	d=aux[n1-k+1];
	m=concat(m,-d*matid(n1)[,(n1+n2)/2+(k+1)/2]));
	for(k=1,(n1-n2)/2,
		d=-content(m[,n2+k]);
		m=concat(m,d*matid(n1)[,n2+k]));
	m
}

{helpskewsymbas()=
print("-------------------------------------------------------------------");
print("      		Computational number Theory                        "); 
print("-------------------------------------------------------------------");
print("");
print("Description: Computes change of basis to achieve symplectic basis of a skew-symmetric matrix");
print("");
print("Original Author:     Fernando Rodriguez-Villegas");
print("                     villegas@math.utexas.edu");
print("                     University of Texas at Austin");
print("");
print("With the assistance of:  Ariel Pacetti");
print("                         apacetti@math.utexas.edu");
print("                         University of Texas at Austin");
print("-------------------------------------------------------------------");
print("");
print("List of routines' helps");
print("");
print("helpsymbas    helpsym    helpmbezout");
print("");
}


addhelp(symbas,"Given skew-symmetric matrix a find matrix m so that m~*a*m is in standard form (like J); in terms of symplectic basis re-writes a in a basis e_1, ..., e_n,f_1,...,f_n (n needs to be even)")

addhelp(sym,"Computes the symmetric form, not the change of base, and the length of a need not be even. ")

addhelp(mbezout,"Computes the gcd of all the entries of the vector v, and gives the 'smaller' linear combination.")