\\ ---------------  GP code  ---------------------------------------
\\
\\ Time-stamp: <Fri, Dec 20, 2002 - 12:40:46 - apacetti>
\\
\\ Description: Routines for arithmetic in quaternion algebras over
\\ the rationals. 
\\
\\
\\ Original Authors:    Fernando Rodriguez-Villegas 
\\                      villegas@math.utexas.edu
\\                      University of Texas at Austin
\\
\\			Ariel Martin Pacetti
\\			apacetti@math.utexas.edu
\\                      University of Texas at Austin
\\
\\ Created:             Wed Feb  12  2000
\\
\\----------------------------------------------------------------------
\\qset(a,b)
\\ Sets the vector qdef=[a,b] defining the quaternion algebra (a,b)/Q
\\ Note that almost all the routines assumes the algebra is ramified
\\ at infinity.
{
qset(a,b)=qdef=[a,b];
}

\\======================================================================
\\ Computes the product of the quaternions represented by the vectors
\\v1 and v2 in the base {1,i,j,ij}

qprod(v1,v2)=
	{local(v);
	v=[v1[1]*v2[1]+qdef[1]*(v1[2]*v2[2]-qdef[2]*v1[4]*v2[4])+
	qdef[2]*v1[3]*v2[3],v1[1]*v2[2]+v1[2]*v2[1]+qdef[2]*(v1[4]*
	v2[3]-v1[3]*v2[4]),v1[3]*v2[1]+v1[1]*v2[3]+qdef[1]*(v1[2]*
	v2[4]-v1[4]*v2[2]),v1[1]*v2[4]+v1[4]*v2[1]+v1[2]*v2[3]-v1[3]*v2[2]];
	v
}

\\======================================================================
\\ Computes the power of a quaternion.

qpower(v,n)=
	{local(ans);
	if(n==0,return([1,0,0,0]));
	ans=v;
	for(k=1,n-1,ans=qprod(ans,v));
	ans}


\\======================================================================
\\ Computes the sum of the quaternions represented by the vectors
\\v1 and v2 in the base {1,i,j,ij}

qsum (v1,v2)=v1+v2

\\======================================================================
\\ Computes the right division of the quaternions numbers v1 and v2.

qrdiv(v1,v2)=qprod(v1,qconj(v2))/qn(v2)

\\======================================================================
\\ Computes the left division of the quaternions numbers v1 and v2.

qldiv(v1,v2)=qprod(qconj(v1),v2)/qn(v1)

\\======================================================================
\\ Computes the multiplicative inverse of v.

qinv(v)=qconj(v)/qn(v)

\\======================================================================
\\ Computes the square of the quaternion number v

qsq(v)=qprod(v,v)

\\======================================================================
\\qconj(v) computes the conjugate of v.

qconj(v)=[v[1],-v[2],-v[3],-v[4]]

\\======================================================================
\\qtr(v) computes the trace of the quaternion v.

qtr(v)=2*v[1]

\\======================================================================
\\ qn(v) computes the norm of v.

qn(v)=qprod(v,qconj(v))[1]

\\======================================================================
\\ to set a lattice L, one must add its norm on the last component

qsetlat(L)=
	{local(aux1,aux2);
	aux1=qgram(concat(L,1));
	aux2=(1/2)*matdiagonal(diag(aux1));
	concat(L,abs(content(aux1-aux2)))
}

\\======================================================================
\\ gives the coordinates of the vector v is in the lattice L. For
\\ references about lattices, see qalg.txt

qcoord(v,L)=matinverseimage(qmat(L),v~)

\\======================================================================
\\ answers if the vector v is in the lattice L. For references about 
\\ lattices, see qalg.txt

qisin(v,L)=
	{local(aux);
	aux=matinverseimage(qmat(L),v~);
	if(length(aux)==0,0,denominator(aux)==1)}

\\======================================================================
\\ Computes the bilinear form in the vectors v1 and v2

qbil(v1,v2)=qtr(qprod(v1,qconj(v2)))

\\======================================================================

\\ computes the matrix of the bilinear form restricted to L, in L's basis

qgram(L)=
	{local(n);
	n=length(L)-1;
        matrix(n,n,k,j,qtr(qprod(L[k],qconj(L[j]))))/L[n+1]
}

\\======================================================================

\\ obtains a basis for the lattice L, and in some way the smallest one
\\ see qflllgram for info

qclean(L)=
	{local(mx,z,n);
	n=length(L)-1;
	mx=qmat(L);
	z=qflllgram(qgram(L));
	concat(qextract(mx*z),L[n+1])
}

\\======================================================================
\\ Makes something like qclean, but it's just for matrices, not for
\\ lattices, so it's get a "minimal" basis in some way

clean(m)=
	{local(z);
        z=qflllgram(m);
        z~*m*z
}

\\======================================================================
\\ Computes just the diagonal of the bilinear form in L basis.

qgramdiag(L)=
	{local(u);
	u=[];
	L=qclean(L);
	for(k=1,4,u=concat(u,qn(L[k])));
	u=vecsort(u/L[5]);u
}

\\======================================================================
\\ Gives the diagonal of the matrix

diag(m)=
        {local(u);
	u=[];
	for(k=1,length(m),u=concat(u,m[k,k]));
	u
}
\\======================================================================
\\ internal, writtes a matrix with the values of L as columns

qmat(L)=
	{local(n,mx);
	n=length(L)-1;
	mx=Mat(L[1]~);
	for(k=2,n,mx=concat(mx,L[k]~));mx
}
\\======================================================================
\\ it's an internal routine to convert a matrix into a lattice

qextract(mat)=
	{local(mx);
	mx=[];
		for(k=1,length(mat),
		mx=concat(mx,[mat[,k]~]));
	mx
}

\\======================================================================
\\ Computes a basis for the lattice spanned by the lattices L1, and L2

qextendlat(L1,L2)=
	{local(mx,n1,n2,aux);
	mx=Mat([]);
	n1=length(L1);n2=length(L2);
        mx=concat(qmat(L1),qmat(L2));
	aux=qextract(mx*qflll(denominator(mx)*mx,1));
	qclean(qsetlat(aux))
}

\\======================================================================
\\ Computes a basis for the lattice spanned by L1, and the vectors in B

qextendvec(L1,B)=
	{local(mx,n1,n2,aux);
	mx=Mat([]);
	n1=length(L1);n2=length(B);
	mx=qmat(L1);
        for(k=1,n2,mx=concat(mx,B[k]~));
	aux=qextract(mx*qflll(denominator(mx)*mx,1));
	qclean(qsetlat(aux))
}

\\======================================================================
\\ computes the conjugate of the lattice m
    
qlatconj(L,k)=concat(vector(length(L)-1,x,qconj(L[x])),L[length(L)])


\\======================================================================
\\ computes the lattice obtained from the product of 2 lattices

qlatprod(L1,L2)=
	{local(mx,n1,n2);
	mx=Mat();
	n1=length(L1)-1;n2=length(L2)-1;
	for(k=1,n1,for(j=1,n2,mx=concat(mx,Mat(qprod(L1[k],L2[j])~))));
	qclean(qsetlat(qextract(mx*qflll(denominator(mx)*mx,1))))
}

\\======================================================================
\\ determines if two lattices are the same

qaresamelat(L1,L2)=
	{local(aux1);
	aux1=(qmat(L1)^(-1)*qmat(L2));
	(denominator(aux1)==1)&&(denominator(aux1^(-1))==1)
}

\\======================================================================
\\ Computes given a invertible matrix m with Q entries, the lattice of
\\ vectors in Q  such that m*v is in Z

qcond(m)=
	{local(n1,L);
	n1=length(m);
	m=1/m;
	L=vector(n1,k,m[,k]~);
	L}

\\======================================================================
\\ Computes given a invertible matrix m with Q entries, the lattice of
\\ vectors in Q  such that v*m is in Z

qrcond(m)=
	{local(n1);
	n1=length(m);
	m=1/m;
	vector(n1,k,m[k,])
	}

\\======================================================================
\\ this is different in the sence that computes which vectors in Z
\\ satisfy the same condition.

qcondz(m)=
	{local(aux1,aux2,n1,L);
	L=[];
	n1=length(m);
	aux1=concat(m,-matid(n1));
	for(k=1,n1,aux1=concat(aux1,vector(2*n1)));
	for(k=1,n1,aux1[k,]=aux1[k,]*(content(aux1[k,])^(-1)));
	aux2=matkerint(aux1);
	aux2=matrix(n1,n1,j,k,aux2[j,k]);
	for(k=1,n1,L=concat(L,[aux2[,k]~]));
	L

}

\\======================================================================
\\ Given a lattice L1, and a vector v, computes the lattice L2 such that 
\\ L2*v \in L1

qlveccond(v,L1)=
	{local(L2,aux1);
	aux1=[v[1],qdef[1]*v[2],v[3]*qdef[2],-v[4]*qdef[1]*qdef[2];
	v[2],v[1],-v[4]*qdef[2],v[3]*qdef[2];v[3],v[4]*qdef[1],v[1],
	-qdef[1]*v[2];v[4],v[3],-v[2],v[1]];
	L2=qsetlat(qcond(qmat(L1)^(-1)*aux1));
	L2
}

\\======================================================================
\\ Given a lattice L1, and a vector v, computes the lattice L2 such that 
\\ v*L2 \in L1

qrveccond(v,L1)=
	{local(L2,aux1);
	aux1=[v[1],qdef[1]*v[2],v[3]*qdef[2],-v[4]*qdef[1]*qdef[2];
	v[2],v[1],v[4]*qdef[2],-v[3]*qdef[2];v[3],-v[4]*qdef[1],v[1],
	qdef[1]*v[2];v[4],-v[3],v[2],v[1]];
	L2=qsetlat(qcond(qmat(L1)^(-1)*aux1));
	L2
}
\\======================================================================
\\ Computes the intersection of two full rank lattices

qlatint(L1,L2)=
	{local(mtx1,mtx2);
	mtx1=qmat(L2);
	mtx2=qmat(qsetlat(qcondz(qmat(L1)^(-1)*mtx1)));
	qclean(qsetlat(qextract(mtx1*mtx2)))
}
	
\\======================================================================
\\ Given a lattice, computes the conjugate of it

qconjlat(m)=
	{local(ans);
	ans=vector(length(m)-1,x,qconj(m[x]));
	ans=concat(ans,m[length(m)]);
	ans}

\\======================================================================
\\ Given a lattice, computes the conjugate of it

qinvlat(m)=
	{local(ans);
	ans=qconjlat(m)/m[length(m)];
	ans[5]=ans[5]/m[length(m)];
	ans}

\\======================================================================
\\ Computes the right order of the lattice, in any case

qrorder(L)=
	{local(L2,k);
	L2=[];
	L2=qrveccond(L[1],L);
	for(k=2,4,L2=qlatint(L2,qrveccond(L[k],L)));
	L2
}

\\======================================================================
\\ Computes the right order of the lattice, in the case the lattice is
\\ the ideal of an Eichler order

qrorder2(L)=
	{local(ans);
	ans=qlatprod(qconjlat(L),L);
	ans=ans/sqrtint(ans[length(ans)]);
	ans[5]=1;
	ans}

\\======================================================================
\\ Computes the left order of the lattice, in any case

qlorder(L)=
	{local(L2);
	L2=[];
	L2=qlveccond(L[1],L);
	for(k=2,4,L2=qlatint(L2,qlveccond(L[k],L)));
	L2
}

\\======================================================================
\\ Computes the right order of the lattice, in the case the lattice is
\\ the ideal of an Eichler order

qlorder2(L)=
	{local(ans);
	ans=qlatprod(L,qconjlat(L));
	ans=ans/sqrtint(ans[length(ans)]);
	ans[5]=1;
	ans}

\\======================================================================
\\ Just changes the scale of the lattice, multiplying the elements of
\\ the basis by the scalar r, and fixing the lattice norm

qscale(L,r)=
	{local(n);
	n=length(L)-1;
	concat(vector(n,k,L[k]/r),L[n+1]/r^2)
}

\\======================================================================
\\ Answers if the lattice L1 is in L2

qlatisin(L1,L2)=
	{local(ans,n1);
	n1=length(L1)-1;
	ans=1;
	for(k=1,n1,ans=ans&&(abs(content(qisin(L1[k],L2)))>=1));
	ans
}

\\======================================================================
\\ Answers if the lattice is an ideal or not

qlatisideal(L)=
	{local(ans,n1);
	ans=1;
	n1=length(L)-1;
	for(j=1,n1,
		for(k=1,n1,ans=ans&&(qisin(qprod(L[j],L[k]),L))));
	ans
}

\\======================================================================
\\ Answers if the lattice is an order or not

qlatisorder(L)=if(qisin([1,0,0,0],L),qlatisideal(L),0)

\\======================================================================
\\ Computes if two lattices m1, m2, are in the same class number. Note
\\ that the order itself is not necesary, because all orders have norm 1.

qareequiv(m1,m2)=
	{local(m);
	m=qlatprod(qconjlat(m1),m2);
	if(qn(m[1])==m[5],1,0)
	}

\\======================================================================
\\ Computes if two lattices m1, m2, are in the same class number. Note
\\ that the order itself is not necesary, because all orders have norm 1.

qareequivr(m1,m2)=
	{local(m);
	m=qlatprod(m2,qconjlat(m1));
	if(qn(m[1])==m[5],1,0)
	}
	

{helpqalg()=
print("-------------------------------------------------------------------");
print("      		Computational number Theory                        "); 
print("-------------------------------------------------------------------");
print("");
print("Description: Routines for arithmetic in quaternion");
print("algebras over the rationals");
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:");
print("");
print("qset	qprod	qsum	qrdiv");
print("qinv	qsq		qconj	qtr");
print("qsetlat	qcoord	qisin	qldiv");
print("qclean	clean	qgramdiag	qextendlat");
print("qextendvec	qlatconj	qlatprod	qaresamelat");
print("qinvlat	qcond	qcondz		qlveccond");
print("qrveccond	qlatint	qlorder	qrorder");
print("qlatscale	qbil	qgram	qlatisin");
print("qn		qlatisideal	qareequiv");
print("");
}

addhelp(qset,"qset(a,b) Sets the vector qdef=[a,b] defining the quaternion  algebra (a,b)/Q.");

addhelp(qprod,"qprod (v1,v2) Computes the product of the quaternions represented  by the vectors v1 and v2 in the base {1,i,j,ij}");

addhelp(qsum,"Computes the sum of the quaternions represented by the vectors v1 and v2 in the base {1,i,j,ij}");

addhelp(qrdiv,"qrdiv (v1,v2) Computes the right division of the quaternions  numbers v1 and v2.");

addhelp(qldiv,"qldiv (v1,v2) Computes the left division of the quaternions  numbers v1 and v2.");

addhelp(qinv,"qinv(v) Computes the multiplicative inverse of v.");

addhelp(qsq,"qsq(v) Computes the square of the quaternion number v");

addhelp(qconj,"qconj(v) computes the conjugate of v.");

addhelp(qtr,"qtr(v) computes the trace of the quaternion v.");

addhelp(qn,"qn(v) computes the norm of v.");

addhelp(qsetlat,"qsetlet (L) adds the lattice norm on the last component");

addhelp(qcoord,"qcoord (v,L)gives the coordinates of the vector v is in the lattice L. For references about lattices, see qalg.txt");

addhelp(qisin,"qisin (v,L)decides if the vector v is in the lattice L. For  references about lattices, see qalg.txt");

addhelp(qbil,"qbil (v1,v2) Computes the bilinear form in the vectors v1 and v2");

addhelp(qgram,"qgram (L) computes the matrix of the bilinear form restricted to L, in L's basis");

addhelp(qclean,"qclean(L) obtains a basis for the lattice L, and in some way the smallest one see qflllgram for info");

addhelp(clean,"clean (M) Makes something like qclean, but it's just for matrices,  not for lattices, so it's get a minimal basis in some way");

addhelp(qgramdiag,"qgramdiag(L) Computes just the diagonal of the bilinear form in L  basis.");

addhelp(qextendlat,"qextendlat (L1,L2) Computes a basis for the lattice spanned by the lattices L1, and L2");

addhelp(qextendvec,"qextend(L1,B)Computes a basis for the lattice spanned by L1,  and the vectors in B");

addhelp(qlatconj,"qlatconj(L) computes the conjugate of the lattice m");

addhelp(qlatprod,"qlatprod(L1,L2) computes the lattice obtained from the product  of the two lattices");

addhelp(qaresamelat,"qaresamelat(L1,L2) determines if two lattices are the same");

addhelp(qcond,"qcond(m) Computes given a invertible matrix m with Q entries,  the lattice of vectors in Q  such that m*v is in Z");

addhelp(qcondz,"qcondz(m) computes which vectors in Z satisfy that m*v is in Z");

addhelp(qlveccond,"qlveccond(v,L1) Given a lattice L1, and a vector v, computes the  lattice L2 such that L2*v \in L1");

addhelp(qrveccond,"qrveccond(v,L1) Given a lattice L1, and a vector v, computes the lattice L2 such that v*L2 \in L1");

addhelp(qlatint,"qlatint(L1,L2) Computes the intersection of two full rank lattices");

addhelp(qrorder,"qrorder(L) Computes the right order of the lattice");

addhelp(qlorder,"qlorder(L) Computes the left order of the lattice");

addhelp(qscale,"qscale(L,r) Just changes the scale of the lattice, multiplying the elements of the basis by the scalar r, and fixing the lattice norm");

addhelp(qlatisin,"qlatisin(L1,L2) Answers if the lattice L1 is in L2");

addhelp(qlatisideal,"qlatisideal(L) Answers if the lattice is an ideal or not");

addhelp(qinvlat,"Given a lattice L, computes it's inverse (which is the conjugate over the norm)");

addhelp(qareequiv,"Given two lattices with the same left order, computes if they are equivalent, i.e. differ from a principal ideal");