
(* 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=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]

(* 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]]]]

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

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

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

(* A matrix with 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]
]

(* 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_]:=CharacteristicPolynomial[m,v]

(* -------------- Real/Complex -------------------- *)
(* Find the real polynomial associated to an irreducible recip poly *)
CxToReal[p_]:=Block[{m,r},
		m=Mat[p];
		r=CharPoly[m+Inverse[m],V[p]];
		r=If[PolyDegree[r]>1
			,r=Expand[PolySqrt[Factor[r]]], r];
		r=If[CoefficientList[r,V[r]][[-1]]>0,r,-r];
		Return[r]
		]

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

(* ---------------- Quadratic Forms -------------- *)
(* Signature *)
Sig[q_]:=Block[{e,n,p},
		e = Eigenvalues[N[q]];
		p = Length[Select[e,#>0.0&]];
		n = Length[Select[e,#<0.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]
]

(* 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 *)
Ch={ { 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_]:= Ch.Minors[Mat[p],2].Transpose[Ch]
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]]]

(* 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}]

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

(* 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 ---------------------- *)
(* 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,"r = " <> ToString[r,InputForm]];
	Write[out,"bnfunit(bnfinit(r))"];
	Write[out,"quit"];
	Close[out];
	Run["gp -q < /tmp/gp.in > /tmp/gp.out"];
	in=ReadList["/tmp/gp.out",String];
	in=in[[2]];
	in=StringTake[in,{3,StringPosition[in,"]"][[1,1]]-1}];
	units=ToExpression["{" <> in <> "}"];
	units=Join[{-1},units];
	Return[units]
]

(* Class number of field Q[p] *)
ClassNo[p_]:=Block[{in,out},
        out=OpenWrite["/tmp/gp.in",
                FormatType->OutputForm,
                PageWidth->Infinity];
        Write[out,"p = " <> ToString[p,InputForm]];
        Write[out,"bnfclgp(p)"];
        Write[out,"quit"];
        Close[out];
        Run["gp -q < /tmp/gp.in > /tmp/gp.out"];
        in=ReadList["/tmp/gp.out",String];
        in=in[[2]];
        in=StringTake[in,{2,StringPosition[in,","][[1,1]]-1}];
        Return[ToExpression[in]];
]


(* ------------------Field theory ---------------------- *)
(* 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]
]

(* 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]
]

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]]

(*----------------------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=SortBy[Expand[units],Height];
	Return[units[[1]]]
]

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

(* -----------------------Search -------------------*)
(* Search[{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 *)
Next[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 *)
Search[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=Next[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]
]

(* -----------------------Search -------------------*)
(* Search[{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 *)
Next[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 *)
Search[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=Next[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]
]

(*----------------------Small Salem Numbers--------------------------*)
SmallSalem[2]:=1 -3x + x^2
SmallSalem[4]:=1 - x - x^2 - x^3 + x^4
SmallSalem[6]:=1 - x^2 - x^3 - x^4 + x^6
SmallSalem[8]:=1 - x^3 - x^4 - x^5 + x^8
SmallSalem[10]:=1 + x - x^3 - x^4 - x^5 - x^6 - x^7 + x^9 + x^10
SmallSalem[18]:=1 + x - x^2 - 3x^3 - 2x^4 + x^6 + x^7 + x^8 + x^9 + 
         x^10 + x^11 + x^12 - 2*x^14 - 3*x^15 - x^16 + x^17 + x^18
