
(*------------------------------  Combined Salem/Coxeter/K3 package ---------------------------------------*)
(* 9/2011 *)

(*------------------------- Salem number routines--------------------------- *)

(* Define Tr if necessary *)

TrDefined= NameQ["Tr"]
If[Not[TrDefined],
	Tr[a_]:=Block[{i},
		Sum[a[[i,i]],{i,Length[a]}]
	]
]

(* --------------- Integers -------------------------*)
(* Factor *) 
StringFactor[n_]:=StringJoin @@ Map[FacToString,FactorInteger[n]]
FacToString[{a_,b_}]:= If[b==1,ToString[a] <> " ",
	StringJoin @@ {ToString[a],"^",ToString[b]," "}]

(* --------------- Polynomials --------------------- *)
(* Extract variable of 1-v polynomial *)
V[p_]:=Block[{vs},
	vs=Sort[Variables[p]];
	Return[If[Length[vs]>0,vs[[-1]],x]]
]

(* Degree of a polynomial *)
PolyDegree[p_]:=Length[CoefficientList[p,V[p]]]-1

(* Test for reciprocity *)
IsRecip[p_]:= Block[{as,bs},
		as=CoefficientList[p,V[p]];
		bs=Reverse[as]*as[[1]];
		Return[as==bs]
]

(* Check for irreducibility of a polynomial *)
Irr[p_]:= Block[{fs},
	fs=FactorList[p];
	Length[fs]==2 && fs[[2,2]]==1]

Irr[p_,prime_]:= Block[{fs},
	fs=FactorList[p,Modulus->prime];
	Length[fs]==2 && fs[[2,2]]==1]

(* Discriminant *)
Disc[p_]:=Det[QTr[p]]

(* Square root of a square polynomial *)
PolySqrt[p_]:=
	Times @@ Map[#[[1]]^(#[[2]]/2)&,FactorList[p]]

(* Ramification of quad extension to real field *)
Ramif[r_]:=Disc[RealToCx[r]]/Disc[r]^2

(* Height *)
Height[p_]:=Max[Abs[CoefficientList[p,V[p]]]]

(* Reduce p modulo q *)
PolyRed[p_,q_]:=PolynomialRemainder[p,q,V[q]]

(* --------------- Matrices ------------------------- *)
IdMat[p_]:=IdentityMatrix[PolyDegree[p]]

Zero[p_]:=0*IdMat[p]

(* Apply a polynomial to a matrix *)
PofM[x_Real,m_]:= x * Id[m]
PofM[x_Integer,m_]:= x * Id[m]
PofM[p_,m_]:= CoefficientList[p,V[p]].Table[MatrixPower[m,i],
		{i,0,PolyDegree[p]}]


(* Apply a rational function to a matrix *)
RofM[r_,m_]:= PofM[Numerator[Together[r]],m].Inverse[PofM[Denominator[Together[r]],m]]

(* Express a as p(b) *)
MtoP[a_,b_]:=Block[{cs,d,i,j,t,ts,p},
                d = Length[a];
                t = Table[Tr[a.MatrixPower[b,i]],{i,0,d-1}];
                ts = Table[Tr[MatrixPower[b,i+j]],{i,0,d-1},{j,0,d-1}];
                cs = Inverse[ts].t;
                p = cs.Table[x^i,{i,0,d-1}];
                Return[p]
]

(* A matrix with a given monic p as its characteristic polynomial *)
Mat[p_]:= Block[{i,m,n,q},
		q=CoefficientList[p,V[p]];
		n=Length[q]-1;
		m=Table[0,{n},{n}];
		Do[m[[i+1,i]]=1,{i,n-1}];
		Do[m[[i,n]]=-q[[i]],{i,n}];
		Return[m]
]

(* For reciprocal polys, action of Z/2 Galois group in same basis *)
AntiDiag[n_]:=Block[{i,j},Table[If[i+j==n+1,1,0],{i,1,n},{j,1,n}]]

Gal[f_]:=Block[{n},
                n = Length[f];
                MatrixPower[f,1-n].AntiDiag[n]
]


(* Turn a row of matrices into a single matrix *)
UnTensor1[r_]:=Transpose[Join @@ Map[Transpose,r]]

(* Turn a matrix of matrices into a single matrix *)
UnTensor[m_]:=Join @@ Map[UnTensor1,m]

(* A matrix with eigenvalues y, where y+1/y =x is a root of p *)
UMat[p_]:= UnTensor[{{Zero[p],IdMat[p]},{-IdMat[p],Mat[p]}}]

(* Short hand for characteristic polynomial *)
CharPoly[m_,v_]:=Det[v IdentityMatrix[Length[m]] - m]

(* -------------- Real/Complex -------------------- *)
(* Find the real polynomial associated to a separable recip poly *)
CxToReal[p_]:=MakeSep[CharPoly[Mat[p] + Inverse[Mat[p]],V[p]]]

(* Make polynomial separable (remove multiple roots) *)
MakeSep[p_]:=Times @@ FactorList[p][[All,1]]

(* Find the Salem polynomial associated to a real polynomial *)
RealToCx[p_]:= CharPoly[UMat[p],V[p]]

(* Test for separability *)
SeparableQ[p_]:=PolyDegree[p]==PolyDegree[MakeSep[p]]

(* ---------------- Quadratic Forms -------------- *)
(* Signature *)
Sig[q_]:=Block[{e,n,p},
		e = Eigenvalues[N[q]];
		e = Chop[e];
		p = Length[Select[e,#>0&]];
		n = Length[Select[e,#<0&]];
		{p,n}]

(* Parity -- is q even or odd? *)
Parity[q_]:=Max[Table[Mod[q[[i,i]],2],{i,1,Length[q]}]]

(* Quadratic forms on the lattice L of algebraic integers spanned 
   by {1,x,x^2,...}, x a root of p.  *)
(* Tr(x.y) *)
QTr[p_]:=Block[{i,j,m,n,q,r,t},
		m=Mat[p];
		n=Length[m];
		r=IdMat[p];
		Do[t[i]=Tr[r];r=r.m,{i,0,2n-2}];
		q=Table[t[i+j],{i,0,n-1},{j,0,n-1}];
		Return[q]
]

(* Tr(x.s(y) *)
(* If p is symmetric, s(x)=1/x is a field automorphism. *)
QHerm[p_]:=Block[{i,j,m,n,q,r,t},
		m=Mat[p];
		n=Length[m];
		r=IdMat[p];
		Do[t[i]=t[-i]=Tr[r];r=r.m,{i,0,n-1}];
		q=Table[t[i-j],{i,0,n-1},{j,0,n-1}];
		Return[q]
]

(* Hermitian form modified by canonical different *)
QEuler[s_]:= Block[{dif,mr,q},
	r=CxToReal[s];
	mr=Mat[s] + Inverse[Mat[s]];
	dif=D[r,V[r]]; 
	q=QHerm[s].Inverse[PofM[dif,mr]];
	Return[q]
]

(* Hermitian form modified by unit and canonical different *)
QDiff[r_,u_]:= Block[{dif,mr,ms,qh,q,s},
	s=RealToCx[r];
	ms=Mat[s];
	mr=ms+Inverse[ms];
	qh=QHerm[s];
	dif=u*D[r,V[r]]; 
	q=qh.Inverse[PofM[dif,mr]];
	Return[q]
]

(*-------------------Euler's and Yang's forms -------------------------*)
(* Symmetric unimodular form, for any p separable *)
(* Make Mat[p] self-adjoint *)
QUnim[p_]:= Block[{dif,mp,q},
	mp=Mat[p];
	dif=D[p,V[p]]; 
	q=QTr[p].Inverse[PofM[dif,mp]];
	Return[q]
]

(* Unimodular form *invariant* under Mat[p] *)
(* p symmetric => q alternating, p anti => q symmetric *)
SUnim[p_]:=Block[{i,j,m,n,q,r,t},
	mp=Mat[p];
	dif=D[p,V[p]]; 
	g=Length[mp]/2;
	q = QHerm[p];
	q= q.MatrixPower[mp,g-1].Inverse[PofM[dif,mp]];
	Return[q]
]

(* Test orthogonality *)
IsOrtho[m_,q_]:=0==Max[Abs[Transpose[m].q.m-q]]

(*-------------------Automorphisms of H3-------------------------*)
(* Basis for Wedge^2(Z^4) is 12, 13, 14, 23, 24, 34.
   Inner products with signs (e.g. 12.34=1, 12.23=0, 13.24=-1) below *)

H3wed=
   {{ 0, 0, 0, 0, 0, 1},
    { 0, 0, 0, 0,-1, 0},
    { 0, 0, 0, 1, 0, 0},
    { 0, 0, 1, 0, 0, 0},
    { 0,-1, 0, 0, 0, 0},
    { 1, 0, 0, 0, 0, 0}}

H3=
   {{ 0, 1, 0, 0, 0, 0},
    { 1, 0, 0, 0, 0, 0},
    { 0, 0, 0, 1, 0, 0},
    { 0, 0, 1, 0, 0, 0},
    { 0, 0, 0, 0, 0, 1},
    { 0, 0, 0, 0, 1, 0}}

(* Change of basis between H3wed and H3 *)
(* H3wed ==Transpose[Ch].H3.Ch *)
BasisCh={ { 1, 0, 0, 0, 0, 0},
     { 0, 0, 0, 0, 0, 1},
     { 0,-1, 0, 0, 0, 0},
     { 0, 0, 0, 0, 1, 0},
     { 0, 0, 1, 0, 0, 0},
     { 0, 0, 0, 1, 0, 0}}

(* From SL_2 to SO(H3) *)
MMat[p_]:= BasisCh.Minors[Mat[p],2].Transpose[BasisCh]
Deg6[p_]:= CharPoly[MMat[p],x]
Deg3[p_]:= CxToReal[Deg6[p]]

(* ------------- Sturm tests for Salem polynomials --------------*)
Salem[p_]:= If[!Irr[p] || !IsRecip[p], False, 
		RealSalem[CxToReal[p]]]

SalemCyc[p_]:= If[!IsRecip[p], False, 
		RealSalemCyc[CxToReal[p]]]

(* Note:  this version is absolutely reliable *)
SlowRealSalem[p_]:= Irr[p] && RealSalemCyc[p] 

(* Note:  this version is for speed *)
RealSalem[p_]:= Irr[p] && RealSalemCyc[N[p]] 

(* Number of roots of p in (a,b) and (c,d) *)
NumRoots[p_,{a_,b_},{c_,d_}]:=Block[{st},
		st=Sturms[p];
		Return[{SturmCount[st,a]-SturmCount[st,b],
		        SturmCount[st,c]-SturmCount[st,d]}]
]

(* Sturm sequence of an irreducible polynomial *)
Sturms[p_]:=Block[{r,st,x},
		x=V[p];
		st={p,D[p,x]};
		While[r= -PolynomialRemainder[st[[-2]],st[[-1]],x];
			!TrueQ[0 == r],
			st=Append[st,r]];
		Return[st]
]

(* Count the number of sign changes in a sequence of +/- 1 *)
SignChgs[s_]:= Sum[1-s[[i]] s[[i+1]],{i,1,Length[s]-1}]/2

(* Sturm count at a given point *)
SturmCount[st_,x_]:=SignChgs[Sign[st/.V[st[[1]]]->x]]

(* Test for real Salem/cyclotomic numbers *)
RealSalemCyc[p_]:= {PolyDegree[p]-1,1}== NumRoots[p,{-2,2},{2,1000}]

(* ------------- Extract Salem factor --------------*)
SalemFactor[p_]:=Select[FactorList[p][[All,1]],Salem][[1]]

(* ------------- Roots --------------*)
(* Numerical roots *)
RtsDigs=20
Rts[p_]:=V[p]/.Chop[N[Solve[p==0],RtsDigs]]

(* Find the root > 1 of a Salem polynomial *)
SalemRoot[p_]:=Block[{r},
		r=RealSalemRoot[CxToReal[p]];
		Return[(r+Sqrt[r^2-4])/2]
]

(* Find the root > 2 of a real Salem polynomial *)
RealSalemRoot[p_]:=Max[Rts[p]]

(*-------------------Lefschetz / Atiyah-Bott-----------------*)
(* Find eigenvalues (a,b) for a given determinant d *)
(* 2Re[d] must be in the range [1-2 Sqrt[2],2] for  *)
(* d, a, b to *all* be on the unit circle.          *)  

(* Positive eigenvalue = action on H^{2,0} *)
HodgeEv[ms_,q_]:=
	Block[{evs,i,norms,vs},
	{evs,vs}=Eigensystem[N[ms,30]]//Chop;
	norms=Map[#.q.Conjugate[#]&,vs]//Chop;
	i=Flatten[Position[Map[#>0&,norms],True]][[1]];
	d=evs[[i]];
	Return[d]
]
	
EVS[d_]:=Block[{a,b}, {a,b}/.
	Solve[	{a*b==d,a+b==(1+d+d^2)/(1+d)},{a,b}][[1]] ]

(*-------------------Check example---------------------------*)
TestEx[{r_,u_}]:=Block[{a,b,d,ms,q},
	If[BadPlaces[r] == 0,
		Print["* Alg Independence May Fail *"];
		Return[False]];

	ms=Mat[RealToCx[r]];
	q=QDiff[r,u];
	If[Abs[Det[q]] != 1 || Sig[q] != {3,19} || 
	   Parity[q] != 0 || !IsOrtho[ms,q] || Max[Denominator[q]] != 1,
	   Printf["* Quad Form Fails *"];
	   Return[False]];

	d=HodgeEv[ms,q];
	{a,b}=EVS[d];
	If[Abs[eigenvalues] != {1,1},
	Print["* Eigenvalues Fail *"];
	Return[False]];

	Return[True]
]

(* Determine the value of tau = d(f)+1/d(f) for a given example *)
Tau[{r_,u_}]:=Block[{ps},
	ps=Places[r];
	ts=Select[ps,(D[r,x] u /.x->#)>0&];
	Return[ts[[1]]]
]

(* ------------------DIALOG WITH GP ---------------------- *)
(* Revised 9/2013 for latest version of GP *)
(* Basis for the group of units in field Q[r] *)
(* Don't forget -1! *)
Units[r_]:=Block[{in,out,units},
	out=OpenWrite["/tmp/gp.in",
		FormatType->OutputForm,
		PageWidth->Infinity];
        Write[out,"field = bnfinit(" <> ToString[r,InputForm] <> ");"];
        Write[out,"field.fu"];
	Write[out,"quit"];
	Close[out];
	Run["gp -q < /tmp/gp.in > /tmp/gp.out"];
	Run["sed 's/\[/{/;s/\]/}/;s/(//g;s/)//g;s/Mod//g' < /tmp/gp.out > /tmp/gp.math"];
	units=<<"/tmp/gp.math";
	units=units[[-1+2*Range[Length[units]/2]]];
	units=Join[{-1},units];
	Return[units]
]


(* ------------------Field theory ---------------------- *)
(* Discriminant of maximal order in number field *)
FieldDisc[p_]:=NumberFieldDiscriminant[V[p]/.Solve[p==0,V[p]][[1]]]

(* Given index of Z[x]/p(x) in its integral closure *)
IntClosure[p_]:=Sqrt[Abs[Disc[p]/FieldDisc[p]]]

(* Class number of field Q[p] *)
(* Mathematica version:  cannot yet do ClassNo[salem10]! *)
(*
ClassNo[p_]:=  NumberFieldClassNumber[Root[p/.V[p]->#&,1]]
*)

(* GP version, 2013 *)
ClassNo[p_]:=Block[{answer,out},
        out=OpenWrite["/tmp/gp.in",
                FormatType->OutputForm,
                PageWidth->Infinity];
        Write[out,"field = bnfinit(" <> ToString[p,InputForm] <> ");"];
        Write[out,"field.clgp.no"];
        Write[out,"quit"];
        Close[out];
        Run["gp -q -s 160000000 < /tmp/gp.in > /tmp/gp.out"];
        answer=<<"/tmp/gp.out";
        Return[answer];
]

(*----------------Info on Salem field and real subfield  -------*)
(* {Z[x]/S[x] : O_K; h(K); same 2 for k; Rel disc of K/k} *)
SalemFieldInfo[salem_]:= Block[{real},
		If[!Irr[salem],Return[{}]];
		real = CxToReal[salem];
		Return[{IntClosure[salem],
			ClassNo[salem],
			IntClosure[real],
			ClassNo[real],
			Ramif[real]}
		]
]


(* FactorOf:  what is the simplest n such that ab=nI ?*)
FactorOf[a_]:=LCM @@ Flatten[Denominator[Inverse[a]]]
(* ------------------DIALOG WITH SAGE ---------------------- *)
(* Factorization of prime *)
FactorIdeal[p_,ring_]:=Block[{out,ans},
        out=OpenWrite["/tmp/sage.in",
                FormatType->OutputForm,
                PageWidth->Infinity];
        Write[out,"sage.rings.number_field.number_field_ideal.SMALL_DISC = 10^40"];
        Write[out,"K.<x> = NumberField(" <> ToString[ring,InputForm] <>")"];
        Write[out,"A = K.order(x)"];
        Write[out,"factor(A.ideal(",p,"))"];
        Write[out,"quit"];
        Close[out];
        Run["sage < /tmp/sage.in | grep Fractional > /tmp/sage.raw"];
	Run["sed 's/^.*sage: /{/;s/(//g;s/))^./))/g;s/$/}/;s/Fractional ideal /{/g;s/))/}/g;s/)/}/g;s/ \* /,/g' < /tmp/sage.raw > /tmp/sage.out"];
	ans = <</tmp/sage.out;
(*	Run["rm -f /tmp/sage.in /tmp/sage.raw /tmp/sage.out"];   *)
	ans = Flatten[ans];
	Return[ans]
]

(* ------------------Field :  Units ---------------------- *)
(* Units + Different => Invariant Forms                  *)
(* Usage:  units=FindUnits[{r,us},{3,19}}]               *)
(*         q = QDiff[r,units[[i]]]                       *) 
(*	   q is even, unimodular, sig (3,19) and         *)
(*	   invariant under Mat[s], s=RealToCx[r]	 *)

Field[r_]:={r,Units[r]}

(* Sort by Salem *)
SortFields[fs_]:=Block[{gs},	
	gs=Map[Join[{RealSalemRoot[#[[1]]]},#]&,fs];
	gs=Sort[gs];
	gs=Map[Drop[#,1]&,gs];
	Return[gs]
]

VectorSpace[n_,p_]:=Table[Mod[Floor[i/p^j],p],{i,0,p^n-1},{j,0,n-1}]

Places[r_]:=Drop[Sort[Rts[r]],-1];

(* Places with d on S^1 but a,b not on S^1 *)
SmallPlaces[r_]:=Length[Select[Rts[r],#<1-2 Sqrt[2]&]]

Signs[p_,places_]:=Map[(1-Sign[p/.V[p]->#])/2&, places];

UnitSigns[{r_,us_}]:= Map[Signs[#,Places[r]]&,us]

DiffSigns[r_]:=Signs[D[r,V[r]],Places[r]]

(* Unit specified by 0s and 1s *)
Unit[{r_,us_},v_]:=
	PolynomialRemainder[Times @@ 
	us[[Flatten[Position[v,1]]]], r,V[r] ]

(* Convert signs to signature *)
SigSigns[s_]:={1,1} + Plus @@ {{2,0},{0,2}}[[1+s]]

(*--------------------------Modular Units-------------------*)
(* Units in real subfield of Z[x], x^m = 1, expressed as polynomials in x. *)
(* Only correct for a odd *)
ModUnit[a_]:=Together[((1-x^a)/(1-x)) x^((1-a)/2)]
ModUnits[m_]:=Map[ModUnit[#]&,Select[Range[3,m-1,2],GCD[#,m]==1&]]

(*----------------------Select Unit--------------------------*)
(* Find all units, up to squares, giving desired signature *)
FindUnits[{r_,us_},{a_,b_}]:=Block[{found,sigs,sd,su,vs},
	su=UnitSigns[{r,us}];
	sd=DiffSigns[r];
	vs=VectorSpace[Length[su],2];
	sigs=Map[SigSigns[Mod[sd+#.su,2]]&,vs];
	vs=vs[[Flatten[Position[sigs,{a,b}]]]];
	found=Map[Unit[{r,us},#]&,vs];
	Return[found]
]

(* Units only in Q[x]! guarantee in Z[x] *)
IntUnits[us_]:=Select[us,1==Max[Denominator[CoefficientList[#,V[#]]]]&]

(* Guarantee Siegel disks *)
SiegelUnits[r_,us_]:=Block[{d,places,tus},
		places=Places[r];
		places=Select[places,#>1-2 Sqrt[2]&];
		d = D[r,V[r]];
		tus=Select[us,Max[# d/.V[r]->places] > 0&];
		Return[tus]
]

(* Select unit*)
PickUnit[{r_,us_}]:=Block[{units},
	units=FindUnits[{r,us},{3,19}];
	units=IntUnits[units];
	units=SiegelUnits[r,units];
	If[units=={},Print["* Unit Search Failed *"];
		Return[]];
	units=SortUnits[Expand[units],Height];
	Return[units[[1]]]
]

SortUnits[u_,f_]:= Transpose[Sort[Transpose[{Map[f,u],u}]]][[2]]

(* -----------------------Search -------------------*)
(* SearchSalem[{a_1,...,a_{deg/2}},{b_1,...b_{deg/2}] *)
(* SearchUn[] picks unramified Salem numbers: |disc(s)/disc(r)^2| = 1 *)
(* a1=-Trace *)
(* Type "Verbose" to each polynomials as found; "UnVerbose" to stop *)
UnVerbose:=Remove["Verbose"]

(* Prepare for a given degree *)
(* Defines RealNo, As and alist *)
SearchSetup[deg_]:=Block[
	{ailist,m,elist,eqs,ram,r1,r2,s,silist},
	alist= Take[{a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14},deg];
	ailist=Take[{a1_,a2_,a3_,a4_,a5_,a6_,a7_,a8_,a9_,a10_,a11_,a11_,a12_,a13_,a14_},deg];
	silist=Drop[ailist,-2];
	elist=Take[alist,-2];

(* Salem numbers with given coefficients *)
	s=Table[x^i,{i,0,2*deg}].Join[{1},alist,Drop[Reverse[alist],1],{1}];
	m=Mat[s];
RealNo[ailist]= Part[Factor[CharPoly[m+Inverse[m],x]],1];

(* Solve for ramification +1/-1 *)
	ram=Factor[Det[m-Inverse[m]]];
	{r1,r2}={Part[ram,2],Part[ram,3]};
	eqs={{r1==1,r2==1},{r1==-1,r2==-1},{r1==1,r2==-1},{r1==-1,r2==1}};
	solns=Table[alist/.Solve[eqs[[i]],elist][[1]],{i,1,4}];
	solns=Select[solns,Max[Denominator[#]]==1&];
As[silist]=solns
]

(* Advance sequence *)
NextSeq[seqi_,low_,hi_]:=Block[{i,seq=seqi},
	seq[[1]]=seq[[1]]+1;
	i=1;
	While[seq[[i]] > hi[[i]],
		If[i==Length[seq],Return[{}]];
		seq[[i]]=low[[i]];
		seq[[i+1]]=seq[[i+1]]+1;
		i=i+1
	];
	Return[seq]
]

(* Find Salem numbers in given region *)
SearchSalem[low_,hi_]:=SearchGen[low,hi,0]
SearchUn[low_,hi_]:=SearchGen[low,hi,1]

(* Flag=0 for regular Salems, 1 for unramified *)
SearchGen[low_,hi_,flag_]:= Block[
	{As,RealNo,alist,deg,new,r,rs,seq},
	deg = Length[low]+2*flag;
	Print["Search region size: ",Times @@ (1+hi-low)];
	SearchSetup[deg];
	Run["echo Starting `date`"];
	rs={}; 
	seq=low;
	While[seq != {},
		rs=Join[rs,TestSalem[seq,flag]];
		seq=NextSeq[seq,low,hi]
	];
	Run["echo Finished `date`"];
	Print["Found ",Length[rs]];
	Return[rs]
]

(* Test and echo result *)
TestSalem[seq_,flag_]:=Block[{as,rs},
	If[flag==0,rs={RealNo[seq]},rs=Map[RealNo,As[seq]]];
	rs=Select[rs,RealSalem];
	If[NameQ["verbose"] && rs != {},
		Print[rs]];
	Return[rs]
]

(*----------------------Generate Salem Numbers by Degree-----------------------*)
(* Only with 0,1,-1 as coefficients *)

NtoDig[n_,i_]:={1,-1,0}[[1+Mod[Floor[n/3^i],3]]]
NtoDigs[n_,m_]:=(NtoDig[n,#]&) /@ Range[0,m-1]
PCoefs[n_,d_]:=Join[{1},NtoDigs[n,d/2], Reverse[NtoDigs[n,d/2-1]],{1}]
Pofn[n_,d_]:=PCoefs[n,d].Table[x^i,{i,0,d}]

PolyGen[deg_]:= Pofn[#,deg]& /@ Range[0,3^(deg/2)-1];
SalemGen[deg_]:=Block[{dim,ss},
	ss = PolyGen[deg];
	ss = Select[ss,Salem];
	Return[ss]
]

(*----------------------Small Salem Numbers by Degree-----------------------*)
salem2:=1 -3x + x^2
salem4:=1 - x - x^2 - x^3 + x^4
salem6:=1 - x^2 - x^3 - x^4 + x^6
salem8:=1 - x^3 - x^4 - x^5 + x^8
salem10:=1 + x - x^3 - x^4 - x^5 - x^6 - x^7 + x^9 + x^10
salem12:=  1 - x + x^2 - x^3 - x^6 - x^9 + x^10 - x^11 + x^12
salem14:= 1 - x^3 - x^4 + x^7 - x^10 - x^11 + x^14
salem16:=  1 - x - x^8 - x^15 + x^16
salem18:=1 - x + x^2 - x^3 - x^6 + x^7 - x^8 + x^9 - x^10 + x^11 - x^12 - x^15 + x^16 - x^17 + x^18
salem20:=1 - x - x^5 + x^6 - x^9 + x^10 - x^11 + x^14 - x^15 - x^19 + x^20
salem22:=1 - x^2 - x^3 + x^7 + x^8 - x^10 - x^11 - x^12 + x^14 + x^15 - x^19 - x^20 + x^22
salem24:= 1 - x - x^4 + x^5 - x^7 + x^8 - x^9 + x^11 - x^12 + x^13 - x^15 + x^16 - x^17 + x^19 - x^20 - x^23 + x^24
salem26:= 1 - x^2 - x^5 - x^8 + x^10 + x^13 + x^16 - x^18 - x^21 - x^24 + x^26

(*----------------------List of Salem Numbers < 1.3--------------------------*)
(* Based on http://www.math.ucla.edu/~mjm/lc/lists/SalemList.html *)
(* Sign corrections have been applied to that list! *)
SmallSalems=
{1 + x - x^3 - x^4 - x^5 - x^6 - x^7 + x^9 + x^10, 
 1 - x + x^2 - x^3 - x^6 + x^7 - x^8 + x^9 - x^10 + x^11 - x^12 - x^15 + 
  x^16 - x^17 + x^18,
  1 - x^3 - x^4 + x^7 - x^10 - x^11 + x^14, 
 1 - x^2 - x^7 - x^12 + x^14,
 1 - x^4 - x^5 - x^6 + x^10,
 1 - x - x^8 + x^9 - x^10 - x^17 + x^18,
 1 - x^3 - x^5 - x^7 + x^10, 
 1 - x - x^5 + x^6 - x^9 + x^10 - x^11 + x^14 - x^15 - x^19 + x^20, 
 1 - x^2 - x^3 + x^7 + x^8 - x^10 - x^11 - x^12 + x^14 + x^15 - x^19 - x^20 + x^22,
 1 - x - x^8 - x^15 + x^16,
 1 - x^2 - x^5 - x^8 + x^10 + x^13 + x^16 - x^18 - x^21 - x^24 + x^26,
 1 - x + x^2 - x^3 - x^6 - x^9 + x^10 - x^11 + x^12,
 1 - x^6 - x^7 - x^8 - x^9 - x^10 - x^11 - x^12 + x^18, 
 1 - x^2 - x^5 - x^15 - x^18 + x^20, 1 - x^2 - x^3 + x^5 - x^7 + x^9 - x^11 - 
  x^12 + x^14, 1 - x - x^4 + x^5 - x^9 + x^13 - x^14 - x^17 + x^18, 
 1 - x - x^4 + x^5 - x^7 + x^8 - x^9 + x^11 - x^12 + x^13 - x^15 + x^16 - 
  x^17 + x^19 - x^20 - x^23 + x^24,
  1 - x - x^3 + x^4 - x^8 + x^9 - x^10 + 
  x^11 - x^12 + x^13 - x^14 + x^18 - x^19 - x^21 + x^22, 
 1 - x^2 - x^5 - x^8 + x^10, 1 - x - x^6 + x^13 - x^20 - x^25 + x^26, 
 1 - x - x^6 + x^7 - x^8 - x^13 + x^14, 1 - x - x^2 + x^3 - x^9 + x^11 - 
  x^13 + x^19 - x^20 - x^21 + x^22, 1 - x^3 - x^4 - x^5 + x^8, 
 1 - x^6 - x^7 - x^8 - x^9 - x^10 - x^11 - x^12 - x^13 - x^14 - x^15 - x^16 - 
  x^17 - x^18 - x^19 - x^20 + x^26, 1 - 2*x + 2*x^2 - 2*x^3 + 2*x^4 - 2*x^5 + 
  x^6 - x^8 + x^9 - x^10 + x^11 - x^12 + x^14 - 2*x^15 + 2*x^16 - 2*x^17 + 
  2*x^18 - 2*x^19 + x^20, 1 - x^4 - x^6 - x^7 - x^9 - x^11 - x^12 - x^14 + 
  x^18, 1 - 2*x + x^2 + x^3 - 2*x^4 + x^5 - x^8 + x^9 - x^11 + x^12 - x^13 + 
  x^14 - x^15 + x^17 - x^18 + x^21 - 2*x^22 + x^23 + x^24 - 2*x^25 + x^26, 
 1 - x^5 - x^6 - x^7 - x^8 - x^9 - x^10 + x^15 - x^20 - x^21 - x^22 - x^23 - 
  x^24 - x^25 + x^30, 1 - 2*x + 2*x^2 - 2*x^3 + x^4 - x^6 + 2*x^7 - 2*x^8 + 
  x^9 - x^11 + x^12 - x^13 + x^14 - x^15 + x^16 - x^17 + x^18 - x^19 + x^21 - 
  2*x^22 + 2*x^23 - x^24 + x^26 - 2*x^27 + 2*x^28 - 2*x^29 + x^30, 
 1 - x - x^8 - x^12 - x^15 - x^18 - x^22 - x^29 + x^30, 
 1 - x^2 - x^3 + x^7 - x^9 - x^10 + x^12 + x^13 + x^14 - x^16 - x^17 + x^19 - 
  x^23 - x^24 + x^26, 1 - x - x^7 - x^11 + x^19 + x^22 + x^25 - x^33 - x^37 - 
  x^43 + x^44, 1 - x^2 - x^5 - x^6 + x^10 + x^13 - x^15 + x^17 + x^20 - 
  x^24 - x^25 - x^28 + x^30, 1 - x - x^4 + x^5 - x^6 + x^8 - x^9 + x^10 - 
  x^12 + x^13 - x^14 + x^16 - x^17 + x^18 - x^20 + x^21 - x^22 + x^24 - 
  x^25 + x^26 - x^28 + x^29 - x^30 - x^33 + x^34, 
 1 - 2*x + 2*x^2 - 2*x^3 + 2*x^4 - 2*x^5 + 2*x^6 - 3*x^7 + 3*x^8 - 3*x^9 + 
  3*x^10 - 3*x^11 + 2*x^12 - 2*x^13 + 2*x^14 - 2*x^15 + 2*x^16 - 2*x^17 + 
  x^18, 1 - x - x^4 + x^5 - x^6 + x^8 - x^9 + x^10 - x^12 + x^13 - x^14 + 
  x^16 - x^17 + x^18 - x^20 + x^21 - x^22 - x^25 + x^26, 
 1 - x - x^6 - x^18 - x^23 + x^24, 1 - x^2 - x^5 - x^8 + x^10 - x^12 - x^15 - 
  x^18 + x^20, 1 - x^3 - x^5 - x^7 - x^9 - x^11 + x^14 + x^16 + x^18 + x^20 + 
  x^22 + x^24 + x^26 - x^29 - x^31 - x^33 - x^35 - x^37 + x^40, 
 1 - x^4 - x^5 - x^6 - x^7 + x^21 + x^22 + x^23 + x^24 + x^25 - x^39 - x^40 - 
  x^41 - x^42 + x^46, 1 - x^2 - x^3 + x^5 - x^7 - x^8 + x^10, 
 1 - x - x^4 + x^5 - x^6 + x^8 - x^9 + x^10 - x^12 + x^13 - x^14 - x^17 + 
  x^18, 1 - x - x^3 + x^5 + x^7 - 2*x^8 + x^11 + x^12 - x^13 - x^14 - x^15 + 
  x^16 + x^17 + x^18 - x^19 - x^20 - x^21 + x^22 + x^23 - 2*x^26 + x^27 + 
  x^29 - x^31 - x^33 + x^34, 1 - x - x^5 + x^11 - x^17 - x^21 + x^22, 
 1 - x^4 - x^5 - x^6 - x^7 - x^8 + x^12 + x^13 + x^14 + x^15 + x^16 - x^20 - 
  x^21 - x^22 - x^23 - x^24 + x^28, 1 + x - x^3 - 2*x^4 - 2*x^5 - x^6 + x^8 + 
  x^9 - x^11 - x^12 + x^14 + x^15 - x^17 - x^18 - x^19 + x^21 + x^22 - x^24 - 
  x^25 + x^27 + x^28 - x^30 - 2*x^31 - 2*x^32 - x^33 + x^35 + x^36, 
 1 - x - x^2 + 2*x^4 - 2*x^6 - x^7 + 2*x^8 + 2*x^9 - 2*x^10 - 2*x^11 + 
  3*x^13 - 2*x^15 - 2*x^16 + 2*x^17 + 2*x^18 - x^19 - 2*x^20 + 2*x^22 - 
  x^24 - x^25 + x^26};

(*-------------------------------Small Salems of Small Degree---------------*)
(* Salem numbers < golden mean, by degree *)
SmallSalems6 = 
	{1-x^2-x^3-x^4+x^6,
	1-x-x^3-x^5+x^6,
	1-x-x^2+x^3-x^4-x^5+x^6,
	1-x^2-2*x^3-x^4+x^6}

SmallSalems8 = 
	{1-x^3-x^4-x^5+x^8,
	1-x+x^2-2*x^3+x^4-2*x^5+x^6-x^7+x^8,
	1-x-x^3+x^4-x^5-x^7+x^8,
	1-x^2-x^3-x^5-x^6+x^8,
	1-x-x^2+x^4-x^6-x^7+x^8,
	1-2*x+2*x^2-3*x^3+3*x^4-3*x^5+2*x^6-2*x^7+x^8,
	1-2*x+x^2-x^4+x^6-2*x^7+x^8}

SmallSalems10 = {
	1+x-x^3-x^4-x^5-x^6-x^7+x^9+x^10,
	1-x^4-x^5-x^6+x^10,
	1-x^3-x^5-x^7+x^10,
	1-x^2-x^5-x^8+x^10,
	1-x^2-x^3+x^5-x^7-x^8+x^10,
	1-x-x^5-x^9+x^10,
	1-x-x^4+x^5-x^6-x^9+x^10,
	1-x-x^3+x^4-x^5+x^6-x^7-x^9+x^10,
	1-x-x^2+x^3-x^5+x^7-x^8-x^9+x^10,
	1-2*x+2*x^2-2*x^3+x^4-x^5+x^6-2*x^7+2*x^8-2*x^9+x^10,
	1-x-x^4-x^6-x^9+x^10,
	1-2*x^2-2*x^3+x^4+3*x^5+x^6-2*x^7-2*x^8+x^10,
	1-x^3-2*x^4-x^5-2*x^6-x^7+x^10,
	1-x-x^2+x^5-x^8-x^9+x^10,
	1-2*x+x^2-2*x^4+3*x^5-2*x^6+x^8-2*x^9+x^10,
	1-x^2-x^3-x^4-x^5-x^6-x^7-x^8+x^10,
	1-2*x+x^2-x^3+2*x^4-3*x^5+2*x^6-x^7+x^8-2*x^9+x^10}

SmallSalems12 = {
	1-x+x^2-x^3-x^6-x^9+x^10-x^11+x^12,
	1-x-x^5+x^6-x^7-x^11+x^12,
	1-x^4-x^5-x^6-x^7-x^8+x^12,
	1-x^2-x^3-x^9-x^10+x^12,
	1-x^2-x^5-x^6-x^7-x^10+x^12,
	1-2*x+2*x^2-2*x^3+x^4-x^5+x^6-x^7+x^8-2*x^9+2*x^10-2*x^11+x^12,
	1-x^3-x^4-x^5-x^6-x^7-x^8-x^9+x^12,
	1-x^2-x^3-x^4+x^6-x^8-x^9-x^10+x^12,
	1-2*x+x^2-x^5+x^6-x^7+x^10-2*x^11+x^12,
	1-2*x+2*x^2-3*x^3+3*x^4-3*x^5+3*x^6-3*x^7+3*x^8-3*x^9+2*x^10-2*x^11+x^12,
	1+x-x^3-2*x^4-3*x^5-3*x^6-3*x^7-2*x^8-x^9+x^11+x^12,
	1-x-x^3-x^9-x^11+x^12,
	1-x-x^4-x^6-x^8-x^11+x^12}

SmallSalems14:= {
	1-x^3-x^4+x^7-x^10-x^11+x^14,
	1-x^2-x^7-x^12+x^14,
	1-x^2-x^3+x^5-x^7+x^9-x^11-x^12+x^14,
	1-x-x^6+x^7-x^8-x^13+x^14,
	1-x-x^3+x^4-x^7+x^10-x^11-x^13+x^14,
	1-x-x^4+x^6-x^7+x^8-x^10-x^13+x^14,
	1-x^2-x^3-x^4+x^6+x^7+x^8-x^10-x^11-x^12+x^14,
	1-x-x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10-x^11-x^13+x^14,
	1-x^2-x^3-x^11-x^12+x^14,
	1-2*x+2*x^2-2*x^3+x^4-x^6+x^7-x^8+x^10-2*x^11+2*x^12-2*x^13+x^14,
	1-x-x^3+x^5-x^6+x^7-x^8+x^9-x^11-x^13+x^14,
	1-x-x^2+x^3-x^5+x^7-x^9+x^11-x^12-x^13+x^14,
	1-x-x^4-x^10-x^13+x^14,
	1-x-x^3+x^7-x^11-x^13+x^14,
	1-x-x^2+x^4+x^5-x^6-x^7-x^8+x^9+x^10-x^12-x^13+x^14,
	1-x^2-x^5-x^6-x^7-x^8-x^9-x^12+x^14,
	1-2*x^2-x^3+x^4+x^7+x^10-x^11-2*x^12+x^14,
	1-2*x+x^2-2*x^6+3*x^7-2*x^8+x^12-2*x^13+x^14,
	1-x-x^2+x^3-x^4+x^6-x^7+x^8-x^10+x^11-x^12-x^13+x^14,
	1-x-x^2+x^4-x^7+x^10-x^12-x^13+x^14,
	1-x^3-x^4-x^5-x^6-x^7-x^8-x^9-x^10-x^11+x^14,
	1-2*x+x^2-x^5+x^6-x^7+x^8-x^9+x^12-2*x^13+x^14,
	1-x^2-x^3-x^6-x^7-x^8-x^11-x^12+x^14}

SmallSalems16 = {
	 1-x-x^8-x^15+x^16,
	1+x-x^3-x^4-x^5-x^6-x^7-x^8-x^9-x^10-x^11-x^12-x^13+x^15+x^16,
	1-x^3-x^5-x^7-x^9-x^11-x^13+x^16,
	1-x-x^6-x^10-x^15+x^16,
	1-x^2-x^3-x^13-x^14+x^16,
	1-x^4-x^5-x^6-x^7-x^8-x^9-x^10-x^11-x^12+x^16,
	1-x-x^2+x^3-x^5+x^6-x^8+x^10-x^11+x^13-x^14-x^15+x^16,
	1-2*x+2*x^2-2*x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10-x^11+x^12-2*x^13+2*x^14-2*x^15+x^16,
	1-x-x^5-x^8-x^11-x^15+x^16,
	1-x-x^3+x^7-x^8+x^9-x^13-x^15+x^16,
	1-x^3-x^4-x^5-x^6-x^8-x^10-x^11-x^12-x^13+x^16,
	1-x-x^5-x^7-x^9-x^11-x^15+x^16,
	1-x-x^4-x^7+x^8-x^9-x^12-x^15+x^16,
	1-2*x+x^2-x^5+x^6-x^8+x^10-x^11+x^14-2*x^15+x^16,
	1-x+x^2-2*x^3+x^4-2*x^5+x^6-2*x^7+x^8-2*x^9+x^10-2*x^11+x^12-2*x^13+x^14-x^15+x^16,
	1-x-x^3-x^5+x^6+x^8+x^10-x^11-x^13-x^15+x^16,
	1-x-x^2+x^4-x^6+x^8-x^10+x^12-x^14-x^15+x^16}

SmallSalems18={
	1-x+x^2-x^3-x^6+x^7-x^8+x^9-x^10+x^11-x^12-x^15+x^16-x^17+x^18,
	1-x-x^8+x^9-x^10-x^17+x^18,
	1-x^6-x^7-x^8-x^9-x^10-x^11-x^12+x^18,
	1-x-x^4+x^5-x^9+x^13-x^14-x^17+x^18,
	1-x^4-x^6-x^7-x^9-x^11-x^12-x^14+x^18,
	1-2*x+2*x^2-2*x^3+2*x^4-2*x^5+2*x^6-3*x^7+3*x^8-3*x^9+3*x^10-3*x^11+2*x^12-2*x^13+2*x^14-2*x^15+2*x^16-2*x^17+x^18,
	1-x-x^4+x^5-x^6+x^8-x^9+x^10-x^12+x^13-x^14-x^17+x^18,
	1-x-x^2+x^3-x^7+x^9-x^11+x^15-x^16-x^17+x^18,
	1-x^5-x^6-x^7-x^8-2*x^9-x^10-x^11-x^12-x^13+x^18,
	1-x^3-x^4-x^8-x^9-x^10-x^14-x^15+x^18,
	1-x-x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10-x^11+x^12-x^13+x^14-x^15-x^17+x^18,
	1-2*x+2*x^2-2*x^3+x^4-x^6+2*x^7-3*x^8+3*x^9-3*x^10+2*x^11-x^12+x^14-2*x^15+2*x^16-2*x^17+x^18,
	1-x-x^4+x^9-x^14-x^17+x^18,
	1-x^2-x^5-x^6-x^9-x^12-x^13-x^16+x^18,
	1-x^4-x^5-x^6-x^7-x^8-x^9-x^10-x^11-x^12-x^13-x^14+x^18,
	1-x-x^2+x^3-x^5+x^7-x^9+x^11-x^13+x^15-x^16-x^17+x^18,
	1-2*x+2*x^2-2*x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10-x^11+x^12-x^13+x^14-2*x^15+2*x^16-2*x^17+x^18,
    	1-x-x^4-x^14-x^17+x^18}

SmallSalems20={
	1-x-x^5+x^6-x^9+x^10-x^11+x^14-x^15-x^19+x^20,
	1-x^2-x^5-x^15-x^18+x^20,
	1-2*x+2*x^2-2*x^3+2*x^4-2*x^5+x^6-x^8+x^9-x^10+x^11-x^12+x^14-2*x^15+2*x^16-2*x^17+2*x^18-2*x^19+x^20,
	1-x^2-x^5-x^8+x^10-x^12-x^15-x^18+x^20,
	1+x-x^3-x^4-x^5-x^6-x^7-x^8-x^9-x^10-x^11-x^12-x^13-x^14-x^15-x^16-x^17+x^19+x^20,
	1-x-x^3+x^5-x^9+x^10-x^11+x^15-x^17-x^19+x^20,
	1-x+x^2-2*x^3+x^4-2*x^5+x^6-x^7+x^8+x^10+x^12-x^13+x^14-2*x^15+x^16-2*x^17+x^18-x^19+x^20,
	1-x-x^3+x^4-x^5+x^6-x^7+x^8-x^9+x^10-x^11+x^12-x^13+x^14-x^15+x^16-x^17-x^19+x^20,
	1-x^2-x^3-x^17-x^18+x^20,
	1-x-x^6-x^10-x^14-x^19+x^20,
	1-2*x+2*x^2-2*x^3+x^4-x^6+x^7-x^8+x^9-x^10+x^11-x^12+x^13-x^14+x^16-2*x^17+2*x^18-2*x^19+x^20,
	1-x-x^2+x^3-x^5+x^6-x^8+x^10-x^12+x^14-x^15+x^17-x^18-x^19+x^20,
	1-x-x^6-x^9-x^11-x^14-x^19+x^20,
	1-x-x^3+x^6-x^10+x^14-x^17-x^19+x^20,
	1-x^4-x^5-x^6-x^7-x^8-x^9-x^10-x^11-x^12-x^13-x^14-x^15-x^16+x^20,
	1-x^3-x^4-x^5-x^9-x^10-x^11-x^15-x^16-x^17+x^20,
	1-x-x^2+x^4+x^5-x^6-x^7+x^10-x^13-x^14+x^15+x^16-x^18-x^19+x^20}

SmallSalems22={
	1-z^2-z^3+z^7+z^8-z^10-z^11-z^12+z^14+z^15-z^19-z^20+z^22,
	1-z-z^3+z^4-z^8+z^9-z^10+z^11-z^12+z^13-z^14+z^18-z^19-z^21+z^22,
	1-z-z^2+z^3-z^9+z^11-z^13+z^19-z^20-z^21+z^22,
	1-z-z^5+z^11-z^17-z^21+z^22,
	1-z-z^4+z^7-z^11+z^15-z^18-z^21+z^22,
	1-z^3-z^4-z^5-z^6+z^9+z^10+z^11+z^12+z^13-z^16-z^17-z^18-z^19+z^22,
	1-z-z^2+z^3-z^7+z^9-z^11+z^13-z^15+z^19-z^20-z^21+z^22,
	1+z-z^3-z^4-z^5-z^6-z^7-z^8-z^9-z^10-z^11-z^12-z^13-z^14-z^15-z^16-z^17-z^18-z^19+z^21+z^22,
	1-z^3-z^5-z^7-z^9-z^11-z^13-z^15-z^17-z^19+z^22,
	1-z^2-z^3-z^4+z^6+z^7+z^8-z^10-z^11-z^12+z^14+z^15+z^16-z^18-z^19-z^20+z^22,
	1-z-z^3+z^4-z^5+z^6-z^7+z^8-z^9+z^10-z^11+z^12-z^13+z^14-z^15+z^16-z^17+z^18-z^19-z^21+z^22,
	1-z^2-z^3-z^19-z^20+z^22,
	1-z-z^3+z^5-z^8+z^11-z^14+z^17-z^19-z^21+z^22,
	1-z^2-z^3-z^4+z^6+z^7-z^11+z^15+z^16-z^18-z^19-z^20+z^22,
	1-z^3-z^4-z^5-z^6-z^7+z^9+z^10+2*z^11+z^12+z^13-z^15-z^16-z^17-z^18-z^19+z^22,
	1-z-z^4+z^10-z^11+z^12-z^18-z^21+z^22,
	1-z-z^2+z^4+z^5-z^6-z^7+z^9-z^11+z^13-z^15-z^16+z^17+z^18-z^20-z^21+z^22,
	1+z-z^3-2*z^4-3*z^5-3*z^6-2*z^7+2*z^9+4*z^10+5*z^11+4*z^12+2*z^13-2*z^15-3*z^16-3*z^17-2*z^18-z^19+z^21+z^22,
	1-z^4-z^5-z^6-z^7-z^8-z^9-z^10-z^11-z^12-z^13-z^14-z^15-z^16-z^17-z^18+z^22,
	1-z-z^5-z^9-z^13-z^17-z^21+z^22,
	1-z-z^4-z^7+z^8+z^11+z^14-z^15-z^18-z^21+z^22,
	1-z-z^2+z^3-z^5+z^7-z^9+z^11-z^13+z^15-z^17+z^19-z^20-z^21+z^22,
	1-2*z+2*z^2-2*z^3+z^4-z^5+z^6-z^7+z^8-z^9+z^10-z^11+z^12-z^13+z^14-z^15+z^16-z^17+z^18-2*z^19+2*z^20-2*z^21+z^22,
	1-z-z^4-z^18-z^21+z^22,
	1-2*z+z^2-z^6+z^7-z^8+z^9-z^11+z^13-z^14+z^15-z^16+z^20-2*z^21+z^22,
	1-z^2-z^3-z^8-z^9+z^11-z^13-z^14-z^19-z^20+z^22,
	1-z-z^2+z^4-z^9+z^11-z^13+z^18-z^20-z^21+z^22,
	1-z-z^3+z^4-z^5-z^7+z^8+z^10-z^11+z^12+z^14-z^15-z^17+z^18-z^19-z^21+z^22,
	1-2*z+z^2-z^5+z^6-z^9+z^10-z^11+z^12-z^13+z^16-z^17+z^20-2*z^21+z^22,
	1-z-z^4+z^5-z^6-z^9-z^11-z^13-z^16+z^17-z^18-z^21+z^22,
	1-z-z^5-z^8-z^10-z^12-z^14-z^17-z^21+z^22,
	1-2*z+z^2-z^4+2*z^5-2*z^6+z^7-z^9+z^10-z^11+z^12-z^13+z^15-2*z^16+2*z^17-z^18+z^20-2*z^21+z^22,
	1-2*z+2*z^2-2*z^3+z^4-z^5+z^6-2*z^7+2*z^8-z^9+z^11-z^13+2*z^14-2*z^15+z^16-z^17+z^18-2*z^19+2*z^20-2*z^21+z^22,
	1-z^3-z^4-z^5-z^6-z^7-z^15-z^16-z^17-z^18-z^19+z^22,
	1-z-z^2+z^4-z^7+z^8-z^11+z^14-z^15+z^18-z^20-z^21+z^22,
	1-2*z+z^2-z^5+z^6-z^8+z^9-z^10+z^11-z^12+z^13-z^14+z^16-z^17+z^20-2*z^21+z^22,
	1-z-z^3+z^7-z^10-z^12+z^15-z^19-z^21+z^22,
	1-z^2-z^3-z^5-z^6+z^8+z^10+z^11+z^12+z^14-z^16-z^17-z^19-z^20+z^22,
	1-z-z^2+z^4+z^5-2*z^6-z^7+z^8+2*z^9-3*z^11+2*z^13+z^14-z^15-2*z^16+z^17+z^18-z^20-z^21+z^22,
	1-z^2-z^3-z^6-z^7+z^9-z^11+z^13-z^15-z^16-z^19-z^20+z^22
}


(* -------------------------------Coxeter routines---------------------------------*)
(*
	DiCox[Q] - to display properties of the quadratic form Q
	and its Coxeter element.

	Spherical forms:  QA[n], QB[n], QC[n], QD[n], QE[n].
	Euclidean and hyperbolic:  QAE[n], QAH[n], etc.
	Y-graphs:  QY[{n1,...,nk}]

	CoxPoly[{i1,...,ik},q] gives the characteristic polynomial
	of the given product of reflections.
			
*)

(* -----------------Examples ---------------- *)
(* Smallest Salem number *)
Ex1:=DiCox[QE[10],"II_9,1"]

(*--------------------------------Cubes and Boxes--------------------------*)
(* Given {d1...dn}, return all {a1,...,an} where ai is in {di1,di2} or {0,di-1} *)
IntegerBox[{d_Integer}]:=Map[{#}&,Range[0,d-1]]
IntegerBox[{{d_,e_}}] := Map[{#}&,Range[d,e]]
IntegerBox[ds_]:=Flatten[Outer[Join,IntegerBox[Take[ds,1]],IntegerBox[Drop[ds,1]],1],1]
IntegerCube[side_,dim_]:=IntegerBox[Table[{0,side-1},{dim}]]

RandomBox[ds_,num_]:=Table[Map[Random[Integer,#]&,ds],{num}]

(* -----------------Matrices---------------- *)
Id[n_Integer]:=IdentityMatrix[n]
Id[v_]:=Id[Length[v]]

CharPoly[m_,v_]:=Det[v Id[m] - m]

(* Turn a row of matrices into a single matrix *)
UnTensor1[r_]:=Transpose[Join @@ Map[Transpose,r]]

(* Turn a matrix of matrices into a single matrix *)
UnTensor[m_]:=Join @@ Map[UnTensor1,m]

(* p x q zero matrix *)
Zero[p_,q_]:=Table[0,{p},{q}] 

(* Direct sum of 2 square matrices *)
DirSum[m_,q_]:=UnTensor[ {{m,Zero[Length[m],Length[q]]},
			{Zero[Length[q],Length[m]],q}}]

(* Boost size of matrix *)
Boost[q_,n_]:=Dirsum[q,Zero[n-Length[q],n-Length[q]]]

(* {Mahler measure,rank} *)
Mahler[F_]:=Block[{as},
		as=Abs[Eigenvalues[N[F]]];
		as=1+Chop[as-1];
		as=Select[as,#>1+10^-6&];
		Return[{Times @@ as,Length[as]}]
]

(* Leading eigenvalue and eigenvector *)
LeadingEV[F_]:=Max[Abs[Eigenvalues[N[F]]]]
LeadingV[F_]:=Block[{es,i,vs,v},
		{es,vs}=Eigensystem[N[F]];
		es = Abs[es];
		i = Position[es,Max[es]][[1,1]];
		v = vs[[i]];
		v = v*Sign[Re[v[[1]]]];
		Return[v]
]

Diagof[q_]:=Map[q[[#,#]]&,Range[Length[q]]]

(* Make diagonal = 2 if necessary *)
NormalizeDiag[q_]:=Block[{m},
		m=DiagonalMatrix[Sqrt[2/Diagof[q]]];
		Return[m.q.m]
]

(* -----------------Reflections/Projections ---------------- *)

ProjV[v_,Q_]:=Outer[Times,v,Q.v]/(v.Q.v)

RefV[v_Integer,Q_]:=RefV[Subspace[v,Q],Q]

RefV[v_,Q_]:=Id[v] - 2*ProjV[v,Q]

Subspace[is_,Q_]:=Id[Q][[is]]

(* Project to kernel of a(f) *)
(* Splits space as Ker(a)+Ker(b); I=ac+bd; then bd = proj to Ker(a) *)
FProj[f_,a_]:=Block[{poly},
		poly = MakeSep[CharPoly[f,x]];
		b = Simplify[poly/a];
		{c,d} = PolynomialExtendedGCD[a,b][[2]];
		proj = PofM[b*d,f];
		Return[proj]
]

(* -----------------Coxeter Element and eigenvalues---------------- *)
(* General form:  Cox = R_A R_B *)
(* If no 2 coloring exists, use s_1...s_n *)

Cox[is_,Q_]:=Dot @@ Map[RefV[#,Q]&,is]

Cox[Q_]:=Cox[Join @@ AnyColor[Q],Q]

CoxPoly[Q_]:= CharPoly[Cox[Q],x]

CoxPoly[is_,Q_]:= CharPoly[Cox[is,Q],x]

(* List all Coxeter matrices *)
Coxs[q_]:=Block[{fs,n,ps},
	n = Length[q];
	ps=Permutations[Range[n-1]];
	ps=Map[Join[#,{n}]&,ps];
	fs=Map[Cox[#,q]&,ps];
	Return[Union[fs]]
]

CoxEV[q_]:=LeadingEV[Cox[q]]
CoxEV[is_,q_]:=LeadingEV[Cox[is,q]]

CoxEVs[q_]:=Map[LeadingEV,Coxs[q]]

(* Random product of n simple reflections *)
FRand[n_,q_]:=Dot @@ Map[RefV[#,q]&,
	Table[Random[Integer,{1,Length[q]}],{n}]]

(* ----------------Colorings---------------------- *)
TwoColor[Q_]:=TwoColor[Q,Null]
TwoColor[Q_,alt_]:=Block[{M,n,ns},
	n=Length[Q];
	M=Sign[Abs[Q]]-IdentityMatrix[n];
	M=M.M;
	M=Sign[MatrixPower[M,n]];
	ns={Select[Range[n],M[[1,#]] != 0&], 
	    Select[Range[n],M[[1,#]] == 0&]};
	If[ns[[1]]=={} || ns[[2]]=={},Return[alt]];
	Return[ns]
]

EvenOdd[Q_]:={2*Range[Length[Q]/2],
		2*Range[Length[Q]/2 + 1/2]-1}

(* Returns 2 coloring or n coloring *)
AnyColor[Q_]:=TwoColor[Q,EvenOdd[Q]]

(* -----------------Adjacency eigensystem---------------- *)
Adj[q_]:=2*Id[q] - NormalizeDiag[q]

AdjV[q_]:=Block[{es,i,v},
		es=Eigensystem[Adj[q]//N];
		i=Position[Re[es[[1]]],Max[Re[es[[1]]]],1][[1,1]];
		v=es[[2,i]];
		Return[v/v[[1]]]
]

AdjEV[q_]:=Max[Eigenvalues[Adj[q]//N]]

(* Bound on Coxeter eigenvalues based on eigenvalues of A *)
BiEV[q_]:=Block[{a,s,t},
		a=AdjEV[q];
		t=2-a^2;
		s=(-t + Sqrt[t^2-4])/2;
		Return[Abs[s]]
]

(* ------------------Generate roots of a root lattice -------------------*)
GenRoots[q_,gens_]:=Block[{i,n,refs,rs},
		n = Length[q];
		rs = Map[Subspace[#,q]&,Range[n]];
		refs = Map[RefV[#,q]&,Range[n]];
		i = gens;
		While[i-- >0, rs = Union[rs,Flatten[Outer[Dot,refs,rs,1],1]]];
		Return[rs]
]

(*---------------Find roots crossed by Coxeter element -----------------*)
CoxRoots[is_,Q_]:=Block[{i,ps,ss,rs,roots},
		rs = Map[Subspace[#,Q]&,is];
		ss = Map[RefV[#,Q]&,rs];
		ps = FoldList[Dot,Id[Q],ss];
		ps = Drop[ps,-1];
		roots = Table[ps[[i]].rs[[i]],{i,1,Length[ps]}];
		Return[roots]
]

CoxRoots[Q_]:=CoxRoots[Flatten[AnyColor[Q]],Q]

(* -----------------Diagram operations---------------- *)
(* Inner product for wt. n in Coxeter diagram *)
Wt[n_]:=-2 Cos[Pi/n]

QDrop[q_,j_]:=Drop[Transpose[Drop[q,{j}]],{j}]

(* Glue node i to j in two Coxeter diagrams *)
(* Default is to glue FIRST nodes together; result remains first *)
QJoin[u_,v_]:=QJoin[u,1,v,1]
QJoin[u_,i_,v_,j_]:=Block[{q},
		q=DirSum[u,QDrop[v,j]];
		q[[Length[u]+Range[Length[v]-1],i]] =
		q[[i,Length[u]+Range[Length[v]-1]]] =
			Drop[v[[j]],{j}];
		Return[q]
]

(* Multiple joins - glue 1st nodes together *)
QJoins[{u_,v_}]:=QJoin[u,v]
QJoins[qs_]:=QJoin[qs[[1]],QJoins[Drop[qs,1]]]

(* Linear diagram *)
(*   a   b   c     z   *)
(* O - O - O - ... - O *)
QL[{}]:={{2}}
QL[{n_}]:={{2,Wt[n]},{Wt[n],2}}
QL[ns_]:=QJoin[QL[Take[ns,1]],2,QL[Drop[ns,1]],1]

(* Join ends of diagram *)
QClose[u_]:=Block[{q=u},
		q[[Length[q],1]]=q[[1,Length[q]]]=-1;
		Return[q]
]

QEdge={{2,-1},{-1,2}}
AddNode[q_,i_]:=QJoin[q,i,QEdge,1]
AddNode[q_]:=AddNode[q,Length[q]]

QDouble[q_]:=QJoin[q,Length[q],q,Length[q]]

(* -----------------Spherical Diagrams---------------- *)
(* A_n O--O--O--O--...--O *)
QA[1]:={{2}}
QA[n_]:=AddNode[QA[n-1]]

(* B_n O==O--O--O--O--...--O *)
QB[2]:={{1,-1},{-1,2}}
QB[n_]:=QJoin[QB[2],2,QA[n-1],1]

(* C_n O==O--O--O--O--...--O *)
QC[2]:={{4,-2},{-2,2}}
QC[n_]:=QJoin[QC[2],2,QA[n-1],1]

(* D_n O--O--O...O--O *)
(*        |           *)
(*        O           *)
QD[3]:={{ 2, 0,-1},
	{ 0, 2,-1},
	{-1,-1, 2}}
QD[n_]:=QJoin[QD[3],3,QA[n-2],1]

(* E_n: O--O--O--O--...--O    *)
(*            |               *)
(*            O               *)
QE[1]:=QA[1]
QE[2]:=QA[2]
QE[3]:=DirSum[QA[2],QA[1]]
QE[4]:=QA[4]
QE[5]:=QD[5]
QE[6]:=QY[{3,3,2}]
QE[7]:=QY[{2,4,3}]
QE[8]:=QY[{2,3,5}]
QE[n_]:=QY[{2,3,n-3}]
QY[ns_]:=QJoins[Map[QA,ns]]

(* F_4 O--O=>=O--O *)
QF[3]:={{ 4,-2, 0},
	{-2, 4,-2},
	{ 0,-2, 2}}
QF[4]:= AddNode[QF[3]]

(* G_2:  O-6-O *)
QG[2]:={{6,-3},{-3,2}}

(* H_3, H_4: O-5-O--O--O *)
QH[3]:=QL[{5,3}]
QH[4]:=QL[{5,3,3}]

(* -----------------Euclidean Diagrams---------------- *)
QAE[n_]:=QClose[QA[n]]
QBE[n_]:=QJoin[QB[n-2],n-2,QD[3],3]
QCE[n_]:=QJoin[QC[n-1],n-1,QC[2],2]
QDE[n_]:=QJoin[QD[n-2],n-2,QD[3],3]
QEE[n_]:=AddNode[QE[n-1]]
QFE[5]:=AddNode[QF[4]]
QGE[3]:=AddNode[QG[2]]

(* -----------------Hyperbolic Diagrams---------------- *)
QAH[n_]:=AddNode[QAE[n-1]]
QBH[n_]:=AddNode[QBE[n-1]]
QDH[n_]:=AddNode[QDE[n-1]]
QEH[n_]:=AddNode[QEE[n-1]]
QK[ns_]:=AddNode[QL[ns],Length[ns]]

(* Q_n:  Loop with 1 edge O==O *)
QQ[n_]:=QClose[QL[Join[{4},Table[3,{n-2}]]]]

(* X_5,6 *) 
QX[5]:=QJoins[{QA[2],QA[2],QA[2],QL[{4}]}]
QX[6]:=QY[{2,2,2,2,2}]

(*----------- Unimodular Lattices --------------*)
(* I_{n,1} *)
QSL[2,1]=QL[{3,Infinity}]

(* II_{n,1} *)
QII[1,1]:={{0,1},{1,0}}
QII[9,1]:= QE[10]
QII[17,1] := QDouble[QE[10]]

(* Complete bipartite graphs *)
(* The first m nodes join the last n nodes *)
QCB[m_,n_]:=Table[If[Min[i,j]<=m && Max[i,j]>m,-1,0],
			{i,m+n},{j,m+n}] +
		2*IdentityMatrix[m+n]

(* The complete graph *)
QCG[n_]:=Table[If[i != j,-1,0],{i,n},{j,n}] +
		2*IdentityMatrix[m+n]

(* -----------------Doubles and Joins---------------- *)
QED[n_]:=QDouble[QE[n]]

(* -----------------Random Graphs---------------- *)
(* d=density of fill; 0<d<1 *)
QRand[n_,d_]:=Block[{fill,q},
		q=2*IdentityMatrix[n];
		fill = n;
		While[fill < d*n^2,
			i=Random[Integer,{1,n}];
			j=Random[Integer,{1,n}];
			If[q[[i,j]]==0,
				fill = fill+2;
				q[[i,j]]=q[[j,i]]=-1
			];
		];
		Return[q]
]

(* -----------------Additional Symmetries ---------------- *)
(* Product of these 2 gives small Salem number sometimes *)
Eq[x_,y_]:=If[x==y,1,0]

Sym[Q_]:=Block[{n,perm,q,S},
	q = Length[Q];
	n = Floor[q/2];
	perm=Join[Range[q-n+1,q],Range[n+1,q-n],Range[1,n]];
	S=Outer[Eq,Range[q],perm];
	Return[S]
]

SubCox[Q_]:=Block[{F,n,U},
	n=Ceiling[Length[Q]/2];
	U=Q[[Range[n],Range[n]]];
	F=Cox[Join @@ TwoColor[U],Q];
	Return[F]
]

SymSub[Q_]:=Sym[Q].SubCox[Q]

(* -----------------Signature of Form ---------------- *)
(* Returns {p,q} *)
Sig[Q_]:=Block[{p,q,s,x},
	s=Eigenvalues[Q];
	s=Chop[s];
	p=Length[Select[s,Re[#]>0&]];
	q=Length[Select[s,Re[#]<0&]];
	Return[{p,q}]
]

(* -----------------Test for even-ness ---------------- *)
EvenForm[q_]:=And @@ EvenQ[Diagonal[q]]

(* -----------------Sample value on the glue group ---------------- *)
FF[Q_]:=Block[{ds,ns},
		If[Det[Q]==0,Return[0]];
		ns = Mod[Diagonal[Inverse[Q]],1];
		ns = Select[ns,Denominator[#]==Max[Denominator[ns]]&];
		Return[Min[ns]]
]
		

(* -----------------Display ---------------- *)
Pad[x_,n_]:= Block[{s},
	s=ToString[x];
	s= s <> StringJoin[Table[" ",{n-StringLength[s]}]];
	s=StringTake[s,n];
	Return[s]
]

DiQ[Q_]:=DiQ[Q,""]
DiQ[Q_,header_]:= Block[{gp},
		Print[
		Pad[header,14],
		" Det ",Pad[Det[Q],6],
		" Sig ",Pad[Sig[Q],8],
		" Glue ",gp = AbGp[Q],
		" FF ",InputForm[FF[Q]]
		]
]

DiCox[Q_]:=DiCox[Q,""]
DiCox[Q_,header_]:=DiCox[Q,header,Cox[Q]]
DiCox[Q_,header_,F_]:=Block[{ev,x},
	DiQ[Q,header];
        ev=LeadingEV[F];
	Print["EV ",Pad[ev,8]," ",InputForm[Factor[CharPoly[F,x]]]]
]

Double[Q_]:=Double[Q,""]
Double[Q_,header_]:= DiCox[QDouble[Q],header,SymSub[QDouble[Q]]]

(* -------GLUING THEORY )--------------*)
(* -------Determine abelian group L/m(L)--------------*)
(* Ans is Z/c[1] + Z/c[2] + ... ; where Z/0  Z *)
AbGp[m_]:=Block[{a,b,c},
	If[! (And @@ Map[IntegerQ,Flatten[m]]),Return[{-1}]];
	{a,b} = IntegerDiagonalize[m];
	c = Abs[Diagonal[a.m.b]];
	c = Join[Table[0,{Length[m]-Length[c]}],c];
	c = Select[c, # != 1&];
	Return[Primary[c]]
]

AbGp[m_,prime_]:=Select[AbGp[m],Mod[#,prime]==0&]

(* Replace numbers by prime power factors *)
Primary[n_Integer]:=Primary[{n}]
Primary[n_List]:=Map[#[[1]]^#[[2]]&,Sort[Join @@ FactorInteger[n]]]

(* ----------------Integer Diagonalize --------------*)
(* Return {a,b} such that a.m.b is upper diagonal *)
(* In Hermite Decomp, a.m = b *)
(* Several iterations may be required! *)
IntegerDiagonalize[m_]:=Block[{a,b,a2,b2,c},
		{a,b} = HermiteDecomposition[m];
		{b,c} = HermiteDecomposition[Transpose[b]];
		b = Transpose[b];
		If[!DiagQ[c],{a2,b2}=IntegerDiagonalize[Transpose[c]];
				a = a2.a; b = b.b2];
		Return[{a,b}]
]

DiagQ[m_]:=Block[{ans,d1,d2,i,j},
		ans = True;
		{d1,d2} = Dimensions[m];
		Do[If[(m[[i,j]] != 0) && (i != j), ans=False],
			{i,1,d1},{j,1,d2}];
		Return[ans]
]

(*-----------Compute ideals----------------*)
(* Ideal[{f,g,...}] computes the Abelian group *)
(* Z^n mod <Im(f), Im(g)> *)

Ideal[ms_]:=AbGp[Transpose[Join @@ Map[Transpose,ms]]]

MaxPeriod = 10000
(* Check period of map *)
Period[f_]:=Block[{evs},
		If[f=={},Return[1]];
		evs = Union[Eigenvalues[N[f]]];
		If[Chop[Max[Abs[evs]]-1] > 0,Return[MaxPeriod]];
		evs = Rationalize[Log[evs]/(2 Pi I)];
		Return[Denom[evs]]
]

(* Check period by brute force in finite group case *)
Period[f_,mod_]:=Block[{g=f,per=1},
		If[f=={},Return[1]];
		While[(g != Id[g]) && (per < MaxPeriod),per++; g = Mod[g.f,mod]];
		Return[per]
]


(*--------------------Glue group routines-------------------*)
(*  We think of L* as lying in L tensor Q *) 
(*  GlueBasis returns {e1,...,en},{d1,...,dn} generators for prod (Z/di) in Q^n/Z^n *)
(*  GlueVectors gives just the ei with di>1 *)
(*  The ei are the columns of (1/Q)B where A(1/Q)B = diag(di) *)
(*  GlueMod gives {d1,...,dn} so group is prod Z/di *)
(*  GlueQ gives fractional form on this group *)
(*  GlueF gives map on this group *)
(*  All can be localized at a prime *)

(* Take just the p part of the glue group *)
Content[0,p_]:=0
Content[n_,1]:=n
Content[n_,p_]:=If[!IntegerQ[n/p],1,p Content[n/p,p]]

GlueBasis[q_]:=GlueBasis[q,1]

GlueBasis[q_,prime_]:=Block[{a,b,cs,ds,es},
		{a,b} = IntegerDiagonalize[q];
		ds = Diagonal[a.q.b];
		es = Transpose[Inverse[a.q]];
		cs = Map[Content[#,prime]&,ds];
		es = es * (ds/cs);
		ds = cs;
		{ds,es} = GlueSort[{ds,es}];
		Return[{ds,es}]
]

(* Sort so ds are ascending *)
GlueSort[{ds_,es_}]:=Block[{is,cs},
			cs = Transpose[{ds,Range[Length[ds]]}];
			cs = Sort[cs];
			is = cs[[All,2]];
			Return[{ds[[is]],es[[is]]}]
]

GlueVectors[q_]:=GlueVectors[q,1]

GlueVectors[q_,prime_]:=Block[{ds,es,n},
			{ds,es} = GlueBasis[q,prime];
			n = Length[Select[ds,#>1&]];
			es = Take[es,-n];
			Return[es]
]

GlueMod[q_]:=GlueMod[q,1]

GlueMod[q_,prime_]:=Select[GlueBasis[q,prime][[1]], #>1&]


GlueQ[q_]:=GlueQ[q,1]

GlueQ[q_,prime_]:=Block[{ds,es,m,n},
		{ds,es} = GlueBasis[q,prime];
		n = Length[Select[ds,#>1&]];
		m = es.q.Transpose[es];
		m = Take[m,-n,-n];
		m = Mod[m,1];
		Return[m]
]

GlueF[q_,f_]:=GlueF[q,f,1]

GlueF[q_,f_,prime_]:=Block[{ds,es,g,n},
		{ds,es} = GlueBasis[q,prime];
		es = Transpose[es];
		n = Length[Select[ds,#>1&]];
		g = Inverse[es].f.es;
		g = Mod[g,ds];
		g = Take[g,-n,-n];
		Return[g]
]

GlueChar[q_,f_,prime_]:=Factor[CharPoly[GlueF[q,f,prime],x],Modulus->prime]

GluePeriod[q_,f_]:=Period[GlueF[q,f],GlueMod[q]]

GluePeriod[q_,f_,prime_]:=Period[GlueF[q,f,prime],GlueMod[q,prime]]

(* Determine minimum dimension with a rotation of order n preserving a lattice *)
GlueDim[n_]:=If[Mod[n,4]==2,GlueDim[n/2],Plus @@ Map[EulerPhi[#[[1]]^#[[2]]]&,FactorInteger[n]]]

(* Find how Ker(p(f)) is glued int L *)
GlueGp[q_,f_,poly_]:=Block[{b,b1,l,l1,q1,r,r1},
		b1 = Kernel[PofM[poly,f]];
		r1 = Length[b1];
		b  = Basis[b1];
		r  = Length[b];
		q1 = b.q.Transpose[b];
		l1 = Take[q1,r1,r1];
		l  = Take[q1,r,r1];
		l =LatticeReduce[l];
		m = l1.Inverse[l];
		Return[AbGp[m]]
]

(*------------------ Integer basis for kernel --------------------*)
(* Note:  NullSpace does not work! *)
Kernel[f_]:=Block[{a,b,n},
		{a,b} = IntegerDiagonalize[f];
		b = Transpose[b];
		n = Select[b,Max[Abs[f.#]]==0&];
		Return[n]
]

Kernel[f_,prime_]:=NullSpace[f,Modulus->prime]

(*------------------ Add glue vectors to get new lattice and auto.----------*)
(* Returns new q,f *)
Glue[qa_,fa_,qb_,fb_,vs_]:= Glue[DirSum[qa,qb],DirSum[fa,fb],vs]

Glue[q_,f_,vs_]:=Block[{f1,es,n,q1},
		n = Length[q];
		es = Id[n];
		es = Join[es,vs];
		es = LatticeReduce[es];
		es = Transpose[es];
		q1 = Transpose[es].q.es;
		f1 = Inverse[es].f.es;
		Return[{q1,f1}]
]

Glue[qa_,fa_,qb_,fb_,vs_]:= Glue[DirSum[qa,qb],DirSum[fa,fb],vs]

(*-----------------Gluing with a map--------------*)
(* phi : G(qa) -> G(qb), anti-isometry *)
Glue[qa_,fa_,qb_,fb_,prime_,phi_]:=Block[{vas,vbs,vs},
	If[phi == {},Return[{}]];
	vas = GlueVectors[qa,prime];
	vbs = Transpose[phi].GlueVectors[qb,prime];
	vs = Map[Join @@ #&, Transpose[{vas,vbs}]];
	{q,f} = Glue[qa,fa,qb,fb,vs];
	Return[{q,f}]
]

(*-----------------Gluing without a map--------------*)
GlueAuto[qa_,fa_,qb_,fb_,prime_]:=
	Glue[qa,fa,qb,fb,prime, FindPhi[qa,fa,qb,fb,prime]]

(* ------------------------Find gluing map-------------------------------*)
(* Work over F_p *)
(* phi.fa = fb.phi *)
Orbit[fs_,v_]:=Map[#.v&,fs]
NewQ[q_,fs_,v_]:=Orbit[fs,v].q.Transpose[Orbit[fs,v]]

FindPhi[qa_,fa_,qb_,fb_,prime_]:=Block[{ga,gb,ma,mb},
	ga = GlueF[qa,fa,prime];
	ma = prime GlueQ[qa,prime];
	gb = GlueF[qb,fb,prime];
	mb = prime GlueQ[qb,prime];
	Return[FindPhiLoc[ma,ga,mb,gb,prime]]
]

FindPhiLoc[ma_,ga_,mb_,gb_,prime_]:=Block[{dim,i,mod,fas,fbs,phi,va,vb,vs},
	dim = Length[ma];
	vs = IntegerCube[prime,dim];
	gas = Table[Mod[MatrixPower[ga,i],prime],{i,0,dim-1}];
	gbs = Table[Mod[MatrixPower[gb,i],prime],{i,0,dim-1}];
	va = Select[vs,MatrixRank[Orbit[gas,#],Modulus->prime]==dim&];
	If[va == {},Print["No cyclic vector found!"];Return[{}]];
	va = va[[1]];
	vb = Select[vs,Max[Mod[NewQ[mb,gbs,#] + NewQ[ma,gas,va],prime]] == 0&];
	If[vb == {},Print["No isometry found!"];Return[{}]];
	vb = vb[[1]];
	phi = Transpose[Inverse[Orbit[gas,va],Modulus->prime].Orbit[gbs,vb]];
	phi = Mod[phi,prime];
	If[Max[Mod[phi.ga-gb.phi,prime]] > 0,
		Print["Failed to conjugate!"]; Return[{}]];
	If[Max[Mod[Transpose[phi].mb.phi+ma,prime]] > 0,
			Print["Not anti-isometric!"]; Return[{}]];
	Return[phi]
]

(*-----------------1D Self-glue --------------*)
GlueSelf[q_,f_,prime_]:=Block[{g,m,mod,vs,v},
		m = GlueQ[q,prime];
		mod = GlueMod[q,prime];
		If[mod != {prime,prime},Print["CANNOT SELF GLUE"];
					Return[{q,f}]];
		vs = IntegerCube[prime,2];
		vs = Drop[vs,1];
		v = Select[vs,Mod[#.m.#,prime]==0&][[1]];
		v = v.GlueVectors[q,prime];
		Return[Glue[q,f,{v}]]
]

(*-----------------One-dimensional gluing --------------*)
Glue1[qa_,fa_,qb_,fb_,prime_Integer]:= Glue1[qa,fa,qb,fb,{prime}]

(* Primes is a list of primes with 1D glue groups *)
Glue1[qa_,fa_,qb_,fb_,primes_]:=Block[{vs},
	vs = Map[Glue1V[qa,fa,qb,fb,#]&,primes];
	{q,f} = Glue[qa,fa,qb,fb,vs];
	Return[{q,f}]
]
		
Glue1V[qa_,fa_,qb_,fb_,prime_]:=Block[{f,ma,mb,mod,phi,q,v},
	ma = GlueQ[qa,prime];
	mb = GlueQ[qb,prime];
	If[Length[ma]+Length[mb] > 2,
		Print["Gluing at ",prime," is not 1D"];
		Return[{}]];
	mod = GlueMod[qa,prime][[1]];
	ma = mod ma[[1,1]];
	mb = mod mb[[1,1]];
	phi = Select[Range[mod],Mod[#*#*ma + mb,mod]==0&];
	If[phi=={},Print["Gluing at ",prime," impossible"]; Return[{}]];
	phi = phi[[1]];
	v = Join[phi*GlueVectors[qa,prime][[1]], GlueVectors[qb,prime][[1]]];
	Return[v]
]

(*-------------------------Restrict (q,f) to Ker poly[f] ---------------------*)
(* Returns (qnew,fnew) for form/map restricted to ker(poly(f)) *)
UnGlue[q_,f_,poly_]:=Take[ UnGlueProj[q,f,poly], 2]

(* Return (qnew,fnew,qproj), qproj = q on projection of original lattice *)
UnGlueProj[q_,f_,poly_]:= GlueSlice[q,f,Kernel[PofM[poly,f]]]

(* Restrict (q,f) to Ker poly[f]^perp *)
UnGluePerp[q_,f_,poly_]:=Take[ UnGluePerpProj[q,f,poly], 2]

UnGluePerpProj[{},f_,poly_]:={{},{},{}}
UnGluePerpProj[q_,f_,poly_]:= Block[{ker},
		ker = Kernel[PofM[poly,f]];
		If[ker=={},Return[{q,f,q}]];
		Return[GlueSlice[q,f,Kernel[ker.q]]]
]

(* Restrict/project to given subspace *)
GlueSlice[{},f_,vs_]:={{},{},{}}
GlueSlice[q_,f_,{}]:={{},{},{}}
GlueSlice[q_,f_,vs_]:= Block[{a,fnew,n,proj,qnew,qproj}, 
		n = Length[vs];
		qnew = vs.q.Transpose[vs];
		a = Basis[vs];
		a = Transpose[a];
		fnew = Inverse[a].f.a;
		fnew = Take[fnew,n,n];
		proj = LatticeReduce[Transpose[vs.q]].Inverse[qnew];
		qproj=proj.qnew.Transpose[proj];
		Return[{qnew,fnew,qproj}]
]

GlueSlicePerp[q_,f_,vs_]:=GlueSlice[q,f,Kernel[vs.q]]

(*-------------------------Extend to a basis---------------------*)
Basis[vs_]:=Block[{basis=vs,es},
		es = Id[Dimensions[vs][[2]]];
		Map[(basis=Basis[basis,#])&,es];
		Return[basis]
]

Basis[vs_,v_]:=Block[{ws},
		ws = Join[vs,{v}];
		If[MatrixRank[ws]>MatrixRank[vs],Return[ws],Return[vs]]
]

(*-------------------------Pull out Pic /Transc ---------------*)
Pic[q_,f_]:=Block[{ps,qpic,fpic},
		ps = GlueFacs[f];
		ps = Select[ps,IsPic[q,f,#]&];
		ps = Times @@ ps;
		{qpic,fpic} = UnGlue[q,f,ps];
		Return[{qpic,fpic}]
]

Transc[q_,f_]:=Block[{ps,qt,ft},
		ps = GlueFacs[f];
		ps = Select[ps,!IsPic[q,f,#]&];
		ps = Times @@ ps;
		{qt,ft} = UnGlue[q,f,ps];
		Return[{qt,ft}]
]

IsPic[q_,f_,p_]:=(Salem[p] && Sig[UnGlue[q,f,p][[1]]][[2]] <= 1) || (!Salem[p] && Sig[UnGlue[q,f,p][[1]]][[2]] == 0)

(*-------------------------Find minimum of a quadratic form ---------------*)
GpRow[v_,end_]:=(StringJoin @@ Map[ToString[#] <> ", "&,Drop[v,-1]]) <> ToString[v[[-1]]] <> end

(* Use gp to find the minimum of a positive definite integral form *)
QMin[q0_]:=Block[{d,i,n,out,q=q0},
	(* gp will only do dim > 1 *)
	If[Length[q0]==0,Return[Infinity]];
	If[Length[q0]==1,Return[q0[[1,1]]]];

	(* Make integral *)
	d = LCM @@ Denominator[Flatten[q]];  
	q = d*q;
	n = Length[q];
	out=OpenWrite["/tmp/gp.in",
                FormatType->OutputForm,
                PageWidth->Infinity];
        Write[out,"m = [\\"];
	Do[Write[out,GpRow[q[[i]],";\\"]],{i,1,n-1}];
	Write[out,GpRow[q[[-1]],"]"]];
	Write[out,"qfminim(m,,1,0)"]; 
        Close[out];
        Run["(gp < /tmp/gp.in | grep \%2 > /tmp/gp.out) >& /dev/null"];
	Run["sed 's/\%2 = \[[0-9]*,//;s/,.*//' < /tmp/gp.out > /tmp/gp.ans"];
	ans = <</tmp/gp.ans;
	Return[ans/d]
]

(* Use gp to find the minimum of a positive definite integral form *)
(* This version works with VERY large integers; returns real number *)
SuperQMin[q0_]:=Block[{d,i,n,out,q=q0},
	(* gp will only do dim > 1 *)
	If[Length[q0]==0,Return[Infinity]];
	If[Length[q0]==1,Return[q0[[1,1]]]];

	(* Make integral *)
	d = LCM @@ Denominator[Flatten[q]];  
	q = d*q;
	n = Length[q];
	out=OpenWrite["/tmp/gp.in",
                FormatType->OutputForm,
                PageWidth->Infinity];
        Write[out,"m = [\\"];
	Do[Write[out,GpRow[q[[i]],";\\"]],{i,1,n-1}];
	Write[out,GpRow[q[[-1]],"]"]];
	(* Use lots of digits to avoid E notation *)
	Write[out,"\\p 100"];
	Write[out,"qfminim(m,,1,2)"]; 
        Close[out];
        Run["(gp < /tmp/gp.in | grep \%2 > /tmp/gp.out) >& /dev/null"];
	Run["sed 's/\%2 = \[[0-9]*,//;s/,.*//' < /tmp/gp.out > /tmp/gp.ans"];
	ans = <</tmp/gp.ans;
	Return[ans/d]
]

(*-------------------------Basic twisting ---------------*)
QTwist[q_,f_,poly_]:=q.PofM[poly,f + Inverse[f]]

QSub[q_,poly_]:=QTwist[q,Cox[q],AdjPoly[poly]]

(* Find q so p(x) p(x^-1) = q(x+1/x) *)
(* This allows one to twist down to the sublattice P(f) *)

AdjPoly[p_Integer]:=p
AdjPoly[p_Real]:=p
AdjPoly[p_]:=Block[{c,f,q},
                c = CoefficientList[p,V[p]][[-1]];
                f= Mat[p/c];
                q = Det[Id[f] - x f + f.f];
                q = Expand[c^2 q];
                Return[q]
]

(*-------------------------Twisting by units to fix signature---------------*)
(* Change twist by units to get signature (p,1) *)
FixSig[q_,f_,poly_]:=FixSig[q,f,poly,1]

(* Change twist by units to get signature (p,k) *)
FixSig[q_,f_,poly_,k_]:=Block[{r,us},
                r = RealPoly[f];
                us = AllUnits[r];
                ps = poly us;
                ps = Select[ps,Sig[QTwist[q,f,#]][[2]] == k&];
		ps = Map[PolyRed[#,r]&,ps];
                Return[ps]
]

(* For all products of distinct elements of xs *)
(* For units this gives reps of U/U^2 *)
Products[{}]:={1}
Products[{x_}]:={1,x}
Products[xs_]:=Join[Products[Drop[xs,1]],xs[[1]] Products[Drop[xs,1]]]

AllUnits[p_]:=Products[Units[p]]

(* --------------Twists for given prime-----------------*)
PrimeTwists[q_,f_,p_]:= FactorIdeal[p,RealPoly[f]]

(* With given signature *)
PrimeTwists[q_,f_,p_,k_]:=
        Join @@ Map[FixSig[q,f,#,k]&,PrimeTwists[q,f,p]]

GlueFacs[f_]:=Drop[FactorList[CharPoly[f,x]],1][[All,1]]

(*-------------------------POSITIVITY--------------------------*)
(* Return:  [code, info] *)
(* code {0,{k}} - signature (p,k), k>= 2 *)
(* code {1,r} - obstructing cyclic root *)
(* code {2,{minr,mini}} - finite order, pos *)
(* code {3,{minq,minr,mini}} - passed lin. alg test *)
(* code {4,{minp,mins,minr,mini}} - passed projection test *)
(* code {5,{minp,mins,minr,mini}} - passed geometric test *)
(* code {6,r} - geometric obstr. root found *)
(* code {7,r} - root found but not obstructive *)

(* mini - minimum fixed vector *)
(* minp - minimum in projection *)
(* minr - minimum rotating vector *)
(* minq - minimum crossing vector *)
(* mins - minimum in Salem factor *)

PosReport[q_,f_]:= Block[{k,fc,fi,fr,fs,minc,mini,minp,minr,mins,qc,qi,qr,qs},
		If[k==Sig[q][[2]] >= 2,Return[{0,{k}}]];

		{qi,fi} = UnGlue[q,f,x-1];
		mini = QMin[qi];
		salem = Select[GlueFacs[f],Salem];
		salem = If[salem=={},1,salem[[1]]];
		{qs,fs} = UnGlue[q,f,salem];
		{qc,fc,qp} = UnGluePerpProj[q,f,salem];
		{qr,fr} = UnGluePerp[qc,fc,x-1];
		minr= QMin[qr];
		If[minr<=2,Return[{1,GetQMin}]];
		If[qs=={},Return[{2,{minr,mini}}]];

		(* Lower bound for norm of crossing obstructions *)
		minq = SalemOb[q,f];
		If[minq > 2,Return[{3,{minq,minr,mini}}]];

		minp = QMin[qp];
		mins = If[q != qs, SalemOb[qs,fs], minq];
		If[minp >= 2 && mins > 2, Return[{4,{minp,mins,minr,mini}}]];

		r = GeomOb[q,f];
		If[r == {}, Return[{5,{minp,mins,minr,mini}}]];
		If[VerifyOb[q,f,r][[2]], Return[{6,r}]];
		Return[{7,r}]
]

(* Return TRUE if positive, FALSE if not or unverified or wrong signature *)
PosQ[q_,f_]:= MemberQ[{2,3,4,5}, PosReport[q,f][[1]]]

(*----------------------Positivity Obstruction:  Salem case----------------*)
(* Give minimum of x.x over obstructing hyperplanes for q,f *)
(* LP Minimizes c.x subject to m.x >= b and x >= l and x are integers *)
(* b: Allow numerically slightly negative numbers except at 1st place where strict positivity is required *)
SalemOb[q_,f_]:=Block[{b,c,embed,l,m,pos,r},
		If[q=={},Return[Infinity]];
                embed = Transpose[QDual[q,f]];
                pos = RealPos[RealPoly[f]];
                c = embed[[1]];
                m = pos.embed;
		m = Chop[m];
                b = -0.0001 + 0*c;
                b[[1]] = 0.0001;
                l = -1000;
                r = embed.LinearProgramming[c,m,b,l,Integers];
                Return[r[[1]]]
]

(* Find separable char polynomial associated to a = f+1/f *)
RealPoly[f_]:=MakeSep[CharPoly[f+Inverse[f],x]]

(* RealPos[p].v gives component of Q(v,v) coming from leading eigenspace of f + 1/f *)
(* Allows positivity to be checked *)
RealPos[p_]:=Block[{d,ps,rs,v,x},
		If[!SeparableQ[p],Print["* * * f IS NOT SEPARABLE * * *"]];
		d = PolyDegree[p];
		x = V[p];
		rs = x/.NSolve[p==0];
		rs = Reverse[Sort[rs]];
		If[Max[rs] <= 2, Print["* * * f IS NOT POSITIVE * * *"]];
		(* Vandermonde matrix : R[x] -> R^d *)
		(* Need special function so 0^0 = 1*)
		v = Vandermonde[rs];
		Return[Transpose[Inverse[v]]]
]

Vandermonde[rs_]:=Outer[If[#2 != 0,#1^#2,1]&,rs,Range[0,d-1]]

(* Pair an element of the lattice with a polynomial in a = f + 1/f *)
(* I.e. compute B(p(a)x,x), where q gives matrix of B *)
Pair[v_,p_,q_,f_]:=v.q.PofM[p,f+Inverse[f]].v

Pair[v_,q_,f_]:=Block[{d,i,n,m,ps,x,xs},
		d = PolyDegree[RealPoly[f]];
		ps = Table[x^i,{i,0,d-1}];
		m = Map[Pair[v,#,q,f]&,ps];
		Return[m]
]

(* Give vectors that span lattice in Z[a]^* *)
Pair[q_,f_]:=Block[{d,i,n,m,ps,x,xs},
		n = Length[q];
		xs = Table[If[i==j,1/2,1]*(Subspace[i,q]+Subspace[j,q]),{i,n},{j,i,n}];
		xs = Flatten[xs,1];
		d = PolyDegree[RealPoly[f]];
		ps = Table[x^i,{i,0,d-1}];
		m = Outer[Pair[#1,#2,q,f]&,xs,ps,1];
		Return[m]
]

(* Maps Z^d to dual, d = Degree of RealPoly[f] *) 
QDual[q_,f_]:=LatticeReduce[Pair[q,f]]

(*--------- Coxeter test for obstructing roots ------------*)
(* Lattice is poly[Cox[q]](L) *)
(* Returns list of roots which still cross *)

CoxOb[q_,poly_]:=Block[{rs,u},
	rs = CoxRoots[q];
	u = PofM[poly,Cox[q]];
	rs = Select[rs,RangeQ[#,u]&];
	Return[rs]
]

CoxRoots[is_,Q_]:=Block[{i,ps,ss,rs,roots},
		rs = Map[Subspace[#,Q]&,is];
		ss = Map[RefV[#,Q]&,rs];
		ps = FoldList[Dot,Id[Q],ss];
		ps = Drop[ps,-1];
		roots = Table[ps[[i]].rs[[i]],{i,1,Length[ps]}];
		Return[roots]
]

CoxRoots[Q_]:=CoxRoots[Flatten[AnyColor[Q]],Q]

Denom[m_]:=LCM @@ Denominator[Flatten[m]]
RangeQ[v_,u_]:= Denom[Inverse[u].v] <= 1

(*--------- Search test for obstructing roots ------------*)
(* {v1,v2} in positive eigenspace, (v1+v2)^2 < 0 *)
(* YES, < 0 is correct for signature (n,1) *)
PosVs[q_,f_]:=Block[{es,vs},
                {es,vs} = Eigensystem[N[f]];
                vs = Sort[Transpose[{Abs[es],vs}]];
                {v1,v2} = vs[[{-1,1},2]];
                If[v1.q.v2 > 0, v2 = -v2];
		Return[{v1,v2}]
]

(* If this is < 0, OBSTRUCTED *)
VerifyOb[q_,f_,r_]:=Block[{t,v1,v2},
		If[r=={},Return[{False,0}]];
		{v1,v2} = PosVs[q,f];
		t = Chop[(v1.q.r)(v2.q.r)];
		Return[{t,t<0}]
]

(*--------- Search for obstructing roots ------------*)
(* 0 = deterministic *)
CubeOb[q_,f_]:= CubeOb[q,f,1,0]

CubeOb[q_,f_,side_]:= CubeOb[q,f,side,0]

(* num = try num times at random *)
CubeOb[q_,f_,side_,num_]:=Block[{dim,ds,rs,ts,v1,v2},
                dim = Length[q];
		ds = Join[{{0,side}},Table[{-side,side},{dim-1}]];
		If[num==0, ts = IntegerBox[ds]];
		If[num >0, ts = RandomBox[ds,num]];
                rs = Select[ts,#.q.# == 2&];
		{v1,v2} = PosVs[q,f];
                rs = Select[rs,Chop[(#.q.v1)(#.q.v2)] < 0&];
                Return[rs]
]

(*--------- Random search ------------*)
(* Note:  we *can* have roots whose planes contain the axis of f *)
(* Thus we need Chop below *)
MaxRand = 10^6
RandOb[q_,f_]:=RandOb[q,f,1]

RandOb[q_,f_,n_]:=Block[{nr,nt,r,v1,v2},
	{v1,v2}=PosVs[q,f];
	nt = nr = 0;
	While[nt < MaxRand,
		nt++;
		r = Table[Random[Integer,{-n,n}],{Length[q]}];
		If[r.q.r == 2, nr++];
	        If[(r.q.r == 2) && (Chop[(r.q.v1)(r.q.v2)] < 0),
			Break[]];
	];
	If[nt<MaxRand,Print["Obstruction found after ",nt," trials"], 
			r={}];
	Return[r]
]

(*------------- Verify all twists are obstructed ------------------*)
(* Return {T/F, {twists, obstr. roots}} *)
AllOb[q_,f_,poly_]:=Block[{fixes,ob},
		fixes = FixSig[q,f,poly,1];
		qs = Map[QTwist[q,f,#]&,fixes];
		rs = Map[RandOb[#,f]&,qs];
		ob = And @@ Map[# != {}&,rs];
		Return[{ob,{fixes,rs}}]
]

AllOb[salem_]:=AllOb[QEuler[salem],Mat[salem],1]

(*------------- Geodesic test for obstructing roots ------------------*)
(* Warning:  will founder on roots with x^perp containing geodesic *)
(* Amount to blowup point on geodesic *)
GeomFactor = 4
GeomErr    = 10^-5

(* Distance in Minkowski model *)
HypDist[a_,b_,q_]:= ArcCosh[ -a.q.b/Sqrt[(a.q.a)(b.q.b)]]

(* Return lower bound for distance from v0 to nearest (root)^perp *)
(* 2 Sinh d(x^perp,a) = <x,a> if x^2=2, a^2 = -1 *)
GeomOb[q_,f_,v0_]:=Block[{d1,d2,form,minform,qv,v=v0},
	v = v/Sqrt[-v.q.v] //N;
	v = Rationalize[v,10^-5]; 
	qv =  Outer[Times,q.v,q.v];
	form = q + GeomFactor * qv;
	minform = SuperQMin[form];
	If[minform <= 2, Return[0]];
	d1 = ArcSinh[Sqrt[(minform-2.0)/(2 GeomFactor)]];
	d2 = HypDist[v,v0,q];
	Return[d1-d2]
]

(* Return {} or obstructing root *)
GeomOb[q_,f_]:=Block[{d,ob,s,t,v,v1,v2},
		{v1,v2} = PosVs[q,f];
		t=0;
		d= Log[Max[Abs[Eigenvalues[f]//N]]];
		While[t <= d,
			v = Exp[t]*v1 + Exp[-t]*v2;
			s = GeomOb[q,f,v];
			t =t + s;
			If[s <= 0, t=Infinity];
		];
		If[t != Infinity,Return[{}],Return[GetQMin]];
]

(* Retrieves vector found by minimizing length *)
GetQMin:= (Run["sed 's/.*\[//;s/^/{/;s/;/,/g;s/\]\]/}/' < /tmp/gp.out > /tmp/gp.m"];
		Return[<</tmp/gp.m])

(* ------------------------Cyclotomic Cousins of Salem polynomials---------------------------*)
(* All primes s.t. salem can have a common factor, mod p, *)
(* with a cyclotomic poly of degree <= dim *)
CycloPrimes[salem_,dim_]:=Block[{cs,gs,ns,primes},
                ns = Range[1000];
                ns = Select[ns,EulerPhi[#]<=dim&];
                cs = Map[Cyclotomic[#,x]&,ns];
                gs = Map[PGCD[salem,#]&,cs];
                gs = Union[Select[gs,#>1&]];
                primes = Union[Flatten[FactorInteger[gs],1][[All,1]]];
                Return[primes]
]

(* Default:  dim = space remaining in K3 *)
CycloPrimes[salem_]:=CycloPrimes[salem,22-PolyDegree[salem]]

(* PGCD[p_,q_]:=Denom[CoefficientList[PolynomialExtendedGCD[p,q],x]] *)
PGCD[p_,q_]:=Abs[Resultant[p,q,x]]

(* Smallest dimension realizing a linear map of period n > 1*)
CycloDim[n_]:=If[Mod[n,4]==2,CycloDim[n/2],Plus @@ EulerPhi[Map[#[[1]]^#[[2]]&,FactorInteger[n]]]]

(* --------------All possible twists (using primes found above)-----------------*)
(* Up to units *)
(* Removes those whose period is too large *)
(* Or whose signature cannot be fixed *)
AllPrimeTwists[salem_]:=Block[{f,polys,primes,q},
	primes = CycloPrimes[salem];
	q = QEuler[salem];
	f = Mat[salem];
	polys = Join @@ Map[PrimeTwists[q,f,#]&,primes];
	polys = Select[polys,DimNeeded[q,f,#] <= 22&];
	Return[polys]
]

(* Allow products of primes, but require correct signature *)
AllTwists[salem_]:=Block[{f,new,polys,q},
	q =QEuler[salem];	
	f = Mat[salem];
	polys = AllPrimeTwists[salem];
	While[{} != (new =TwistEvol[q,f,polys]),
		polys = Join[polys,new]];
	Print["All twists:  those with dim'n small enough: ",Length[polys]];
	polys = Select[polys,GluePeriods[QTwist[q,f,#],f] != {}&];
	Print["All twists:  those with period small enough: ",Length[polys]];
	polys = Select[polys,FixSig[q,f,#,1] != {}&];
	Print["All twists:  those with right signature: ",Length[polys]];
	Return[polys]
]

(* Fast dimension test, but sometimes more dim is needed *)
DimNeeded[q_,f_,poly_]:=
	Length[q] + CycloDim[GluePeriod[QTwist[q,f,poly],f]] 

TwistN[f_,poly_]:=Denom[Inverse[PofM[poly,f+Inverse[f]]]]

TwistEvol[q_,f_,polys_]:=Block[{new},
	new = Flatten[Outer[Times,polys,polys,1],1];
	new = Union[new];
	new = Complement[new,polys];
	new = Select[new,DimNeeded[q,f,#] <= 22&];
	Return[new]
]

(* ----------------------Cyclotomic Cousin of (Q,F) -----------------*)
(* GluePeriods[q,f]  --- lists all possible {n1,...,nk} such that *)
(* {q,f} can be mated with these cyclotomic periods and still have dim <= 22 *)

(* Glue map g must satisfy P(g) = 0 *)
(* P(g) = product of *distinct* cyclotomic polynomials C(n_i) *)
(* sum phi(n_i) = dimension *)

(* Returns lists of all increasing sequence s.t. sum phi(n_i) = dim *)
CycloNs[dim_]:=CycloNs[dim] = Block[{ns},
	ns = Range[1000];
	ns = Select[ns,EulerPhi[#]<=dim&];
	ns = Map[{#}&,ns];
	ns = Join @@ Map[CycloX[#,dim]&,ns];
	Return[ns]
]

CycloNs[0]={}
CycloNs[1]={{1,2}}
CycloNs[2]={{1,1},{3},{4},{6}}

CycloX[ns_,dim_]:=Block[{ks,rem},
		rem = dim - (Plus @@ EulerPhi[ns]);
		If[rem==0,Return[{ns}]];
		ks = CycloNs[rem];
		ks = Select[ks,Min[#] > Max[ns]&];
		ks = Map[Join[ns,#]&,ks];
		Return[ks]
]

(* All s.t. sum phi(n_i) <= dim *)
AllCycloNs[dim_]:=Join @@ Map[CycloNs[#]&,Range[dim]]

(* Find possible periodicity *)
GluePeriods[q_,f_]:=Block[{dim,g,mod,ns},
		g = GlueF[q,f];
		mod = GlueMod[q];
		dim = 22 - Length[q];
		ns = AllCycloNs[dim];
		ns = Select[ns,GlueKilled[g,mod,#]&];
		ns = MinNs[ns];
		Return[ns]
]

(* See if these periods kill the glue map *)
GlueKilled[g_,mod_,ns_]:=
	Max[Mod[PofM[Times @@ Map[Cyclotomic[#,x]&,ns],g],mod]] == 0

(* Test if a is a subset of b *)
SubQ[a_,b_]:=Length[Union[a,b]]==Length[b];

(* Test if a is a minimal subset among the as *)
MinQ[a_,as_]:=And @@ Map[(# == a) || !SubQ[#,a]&,as]

(* Extract minimal elements *)
MinNs[ns_]:=Select[ns,MinQ[#,ns]&]

(* Compute rotation in transcendental lattice *)
K3Rot[q_,f_]:=Block[{a,es,ns,r,u,vs},
		a = f + Inverse[f];
		{es,vs} = Eigensystem[a];
		ns = Sign[Map[#.q.#&,vs]];
		i = Position[ns,-1][[1,1]];
		r = es[[i]];
		r = u/.Solve[u + 1/u == r,{u}][[1]];
		r = Log[r]/(2Pi I);
		r = Chop[N[r]];
		r = Rationalize[r];
		Return[Abs[r]]
]

(*-------------------------Display characteristic polynomials ---------------*)
DiChar[q_,f_]:=Block[{per,poly,primes},
		per = Period[f];
		poly = Factor[CharPoly[f,x]];
		If[per < MaxPeriod,
			Print["Char poly  ",poly,"  Period ",per],
			Print["Char poly  ",poly]
		];
		primes = FactorInteger[Det[q]][[All,1]];
		primes = Select[primes, #>1&];
		Map[Print["Localization at ",#,": ",GlueChar[q,f,#]]&,primes];
]

(*-------------------------Display gluing info ---------------*)
DiF[q_,f_]:=Block[{null,sig},
		null = Length[NullSpace[q]];
		sig = Sig[q][[2]];
		If[(null > 0) || (sig > 1), DiF2[q,f]; Return[]];
		If[sig == 0, DiF0[q,f],DiF1[q,f]]
]

(* For weird signatures or nullities *)
DiF2[q_,f_]:=Block[{lambda},
		lambda = Max[Abs[Eigenvalues[N[f]]]];
		Print["Glue group ",AbGp[q],", Signature ",Sig[q],", lambda = ",lambda];
		DiChar[q,f]
]

(* Definite form *)
DiF0[q_,f_]:= (Print["Glue group ",AbGp[q],", Signature ",Sig[q]]; DiPos[q,f]; DiChar[q,f])

(* Hyperbolic form *)
DiF1[q_,f_]:=Block[{lambda},
	lambda = Max[Abs[Eigenvalues[N[f]]]];
	If[lambda < 1.001,DiF2[q,f]; Return[]];
	Print["Glue group ",AbGp[q],", Signature ",Sig[q],", lambda ",lambda];
	DiPos[q,f];
	DiChar[q,f]
]

DiGlue[q_,f_,p_]:=Block[{d,g,h,m,mod,per,sig},
                mod = GlueMod[q,p];
                g = GlueF[q,f,p];
                m = GlueQ[q,p];
		per=Period[g,mod];
                d = LCM @@ Flatten[Denominator[m]];
                Print["At prime ",p," : Glue group ",mod,", Period ",per,", Char poly ",GlueChar[q,f,p]];
                Print["                 q = 1/",d," * ",MatrixForm[d m],"   f = ",MatrixForm[g]]
]

DiGlue[q_,f_]:=Block[{p,ps},
		DiF[q,f];
                ps = Transpose[FactorInteger[Abs[Det[q]]]][[1]];
                If[ps=={1},Print["Unimodular"];ps={}];
                Map[DiGlue[q,f,#]&,ps];
]       

(* Pass to sublattice defined by poly(f) *)
DiSubCox[q_,poly_]:=Block[{f,rs},
		Print["               Sub Coxeter System ",poly];
		f = Cox[q];
		DiF[QTwist[q,f,AdjPoly[poly]],f];
		Print["Number of obstructing roots: ",Length[CoxOb[q,poly]]]
]

(* ------------------------------------Display full K3 lattice---------------------------------*)

DiK3[q_,f_]:=DiK3[q,f,1]

DiK3[q_,f_,polys_]:=Block[{fpic,ft,qpic,qt},
		(* Check basic features *)
		If[Denom[Join[q,f]] > 1,	Print["Q/F ARE NOT INTEGRAL"];Return[]];
		If[!EvenForm[q],           	Print["Q IS NOT EVEN"];Return[]];
		If[AbGp[q] != {},           	Print["Q IS NOT UNIMODULAR"];Return[]];
		If[Sig[q] != {19,3},		Print["Q IS NOT SIG (19,3)"];Return[]];
		If[Transpose[f].q.f != q,	Print["Q IS NOT F-INVARIANT"];Return[]];
		If[Abs[Det[f]] != 1,		Print["F IS NOT INVERTIBLE"];Return[]];

		{qpic,fpic} = Pic[q,f];
		{qt,ft} = Transc[q,f];

		Print["\n		Picard group "];
		DiF[qpic,fpic];

		Print["\n		Transcendental cycles"];
		DiF[qt,ft];
		Print["Rotation of T(X): ",InputForm[K3Rot[qt,ft]]];

		DiPic[qpic,fpic,polys];
]

DiPos[q_,f_]:=Block[{code,r},
		{code,r} = PosReport[q,f];
		Print[""];
		Switch[code,
		0, Null,
		1, Print["* * * * POSITIVITY FAILED:  Cyclic root ",r],
		2, {minr,mini} = r; 
			Print["Positivity verified.  Cyclic: ",minr," Fixed: ",mini],
		3, {minq,minr,mini} = r; 
			Print["Positivity verified.  Crossing: ",minq," Cyclic: ",minr," Fixed: ",mini],
		4, {minp,mins,minr,mini} = r; 
			Print["Positivity verified.  Proj: ",minp," Salem: ",mins," Cyclic: ",minr," Fixed: ",mini],
		5, {minp,mins,minr,mini} = r; 
			Print["Positivity verified geometrically.  Proj: ",minp," Salem: ",mins," Cyclic: ",minr," Fixed: ",mini],
		6, Print["* * * * POSITIVITY FAILED:  Obstructing root ",r],
		7, Print["* * * * POSITIVITY UNVERIFIED."]
		]
]

(*-------------------Picard group report-----------------*)
DiPic[qpic_,fpic_]:= DiPic[qpic,fpic,1]

DiPic[qpic_,fpic_,polys_]:=Block[{ps,fc,fs,qc,qs,salem},
		If[Sig[qpic][[2]] == 0, Return[]];
		ps = GlueFacs[fpic];
		salem = Select[ps,Salem];
		If[salem=={},Return[]];

		salem = salem[[1]];
		{qs,fs} = UnGlue[qpic,fpic,salem];
		{qc,fc} = UnGluePerp[qpic,fpic,salem];
		If[qc != {},
			Print["\n		Periodic factor of Pic"];
			DiF[qc,fc]];
		Print["\n		Salem factor of Pic"];
			DiF[qs,fs];
		Print["\n		Field Information"];
		DiField[salem];
		DiIdeal[salem,polys];
]

(*-------------------Field and ideal properties -----------------*)
(* Work for reciprocal polynomials *)

DiField[salem_]:= Print["A/Z[x], h(K), B/Z[y], h(k), Ramif: ",
		SalemFieldInfo[salem]];

DiIdeal[salem_,poly_]:=Block[{facs,idealr,ideals, prime,real},
		If[poly==1,Return[]];
		If[!Irr[salem],Print["Not a field"]; Return[]];
		Print["Twisted by ",poly];
		real = CxToReal[salem];
		facs = FactorList[poly][[All,1]];
		facs = Select[facs, PolyDegree[#]>0&];
		Map[DiPrimeIdeal[real,#]&,facs];
]

DiPrimeIdeal[real_,poly_]:=Block[{a,deg,prime,idealr,ideals,x,y},
		a = PofM[poly,Mat[real]];
		{prime,deg}=FactorInteger[Abs[Det[a]]][[1]];
		If[prime==1,Return[]];
		idealr = PolynomialGCD[real,poly,Modulus->prime];
		ideals = Expand[y^PolyDegree[idealr] (idealr/.x->y+1/y)];
		ideals = Factor[ideals,Modulus->prime];
		If[Irr[ideals,prime],
			Print["Prime ",prime,"^",deg,"  <-  (",idealr,"), inert in K"],
			Print["Prime ",prime,"^",deg,"  <-  (",idealr,"), splits in K"]];
]

DiPrimes[salem_]:=Block[{polys,primes,real},
	real = CxToReal[salem];
	primes = CycloPrimes[salem];
	Print[""];
	Print["		Primes and Twists"];
	Print[""];
	Print["Feasible Primes: ",primes];
	polys = Join @@ Map[FactorIdeal[#,real]&,primes];
	Map[DiPrimeIdeal[real,#]&,polys]
]

(*-------------------Groundwork for realizing a given Salem number -----------------*)
DiK3Salem[]:=Map[DiStrat,Select[SmallSalems,PolyDegree[#]<=20&]]

DiStrat[salem_]:=(
		Print[""];
		Print["		Salem number ",SalemRoot[salem]," of degree ",PolyDegree[salem]];
		Print[""];
		DiField[salem];
		DiK3Salem[salem]
)

DiK3Salem[salem_]:=Block[{f,polys,q},
		DiPrimes[salem];

		polys = AllTwists[salem];
		Print["Twists found: ",Length[polys]];
		If[polys=={},Print["NOT REALIZABLE"];
				Return[]];

		q = QEuler[salem];
		f = Mat[salem];
		Map[DiTwist[q,f,#]&,polys];
]

DiTwist[q_,f_,poly_]:=Block[{fix,n,per,primes,rs,qt,str},
		qt = QTwist[q,f,poly];
		per = GluePeriods[qt,f][[1]];
		n = TwistN[f,poly];
		fix = FixSig[qt,f,1][[1]];
		qt = QTwist[qt,f,fix];
		min = SalemOb[qt,f];
		str = "";
		If[min <= 2,
			If[RandOb[qt,f] != {}, str = "OBSTRUCTED"]];
		Print["Ideal ", Pad[StringFactor[n],10],
			" Periods ", Pad[per,14],
			" Min ",Pad[min,2]," ",str];
]



