
(* -------------------------Candidate for Bedford-Kim example---------------------------*)
Print["\n\n             CANDIDATE FOR BEDFORD-KIM EXAMPLE\n"]

<<cox.m


(*-------------------Realize lambda_8 -------------------*)
salem = salem8
cyclo1 = Cyclotomic[1,x]
cyclo2 = Cyclotomic[2,x]
cyclo6 = Cyclotomic[6,x]
cyclo18 = Cyclotomic[18,x]

p = 3
(*-------------------Find viable extra twists -------------------*)
Extras[q_,f_]:=Block[{a,ps},
	a = f + Inverse[f];
	ps = HalfCube[5,5].{1,x,x^2,x^3,x^4}; 
	ps = Union[ps,-ps];
	ps = Select[ps,0<Abs[Det[PofM[#,a]]] <= 100&];
	ps = Select[ps,Sig[QTwist[q,f,#]]=={7,2}&];
	ps = Select[ps,EvenForm[QTwist[q,f,#]]&];
	ps = Select[ps,TestGlue[QTwist[q,f,#],f]==2&];
	Return[ps]
]

(* Find matrix conjugating action on Z/9 x Z/27 to standard form *)
(* h = {{0,1},{3,0}} and m = +/-{{3,0},{0,1}} *)

FindGlue[m_,g_]:=Block[{a,b,e1,e2,es,h},
		h = g + 2 Id[g];
		If[Max[Mod[h.h-3Id[h],mod]] > 0,	
			Print[" * * * Incorrect module * * *"]];
		e2 = {0,1};
		e1 = h.e2;
		es = {e1,e2};
		es = Transpose[es];
		a = 27 (Transpose[es].m.es)[[2,2]];
		b=Select[Range[0,26],Mod[#*#*a-1,27]==0 || Mod[#*#*a+1,27]==0 &][[1]];
		es = Mod[b es,mod];
		Return[es]
]

(* Test action mod 3 (returns 1 or 2) *)
TestGlue[q_,f_]:=Block[{g,m,es},
		g = GlueF[q,f,3];
		m = GlueQ[q,3];
		es = FindGlue[m,g];
		m = Transpose[es].m.es;
		m = Mod[27 m,3];
		Return[m[[2,2]]]
]

(*-------------------Salem Part-------------------*)
qs = QEuler[salem];
fs = Mat[salem];
poly = (2+x)^2      (* This makes glue group Z/9 x Z/27 *)
(* twists = FixSig[qs,fs,(2+x)^2] *)

(* Two possible twists.  Only one will agree with projective BK map *)
twists =  {(2 + x)^2*(-2*x - x^2 + x^3), -(x*(1 + x)*(2 + x)^2)}
polys = twists[[2]]
qs = QTwist[qs,fs,polys];
ms = GlueQ[qs,p]
gs = GlueF[qs,fs,p]
vss = GlueVectors[qs,p]
mod = GlueMod[qs,p]
phis = FindGlue[ms,-gs]

Print["Note that glue form represents 1"]
DiGlue[qs,fs]

(*-------------------Cyclotomic Part-------------------*)
fc = Mat[x^9-1]; 
qc = Id[9];
polyc = (2x^2 - 5);
qc = QTwist[qc,fc,polyc];

(* Look at possible extra twists *)
(* extras = Extras[qc,fc]  *)
extras = 
	{-4 - 2*x + 3*x^2,
	-4*x + x^2 + x^3,
	-2 - 3*x + x^2 + x^3,
	-x^2 - x^3 + x^4, 
	2 - 5*x - 3*x^2 + x^3 + x^4,
	-3*x - 2*x^2 + x^3 + x^4,
	-4*x^2 - 2*x^3 + 3*x^4,
	 2*x - 5*x^2 - 4*x^3 + 4*x^4}

(* Pick one and twist by it*)
(* This fixes signature to (7,2) and fixed sign mod 3 *)
(* but adds more to the glue group *)
(* extra = extras[[1]] *)
(* extra = (x^2-x-1)(x^2+x-1)(x^2-2) *)
extra = extras[[1]]
qc = QTwist[qc,fc,extra]

mc = GlueQ[qc,p]
gc = GlueF[qc,fc,p]
vsc = GlueVectors[qc,p]
phic = FindGlue[mc,gc]

(*-------------------qs and qc -------------------*)
q0 = DirSum[qs,qc]
f0 = DirSum[fs,-fc]

{s1,s2} = Transpose[phis].vss
{c1,c2} = Transpose[phic].vsc

v1 = Join[s1,c1]
v2 = Join[s2,c2]

{qi,fi} = Glue[q0,f0,{v1,v2}]

(*-------------------Retreat from Z/4 to Z/27 as glue group-------------------*)
{qj,fj} = UnGlue[qi,fi,salem cyclo6 cyclo18]

(*-------------------Tail lattice and map -------------------*)
qa1 = DiagonalMatrix[{2,2,4*3,4*27}]
fa1 = DiagonalMatrix[{1,1, -1,  -1}]

qa2 = QA[2] 
fa2 = -Cox[qa2]

qt = DirSum[qa1,qa2]
ft = DirSum[fa1,fa2]

v1 = {1/2,   0, 1/12,1/4,-1/3, 1/3}
v2 = {  0, 1/2, 1/4, 1/4,   0,   0}

{qx,fx} = Glue[qt,ft,{v1,v2}]

(*-------------------Complete K3 -------------------*)
{q,f} = Glue1[qj,fj,qx,fx,3]

DiK3[q,f,polys]

(*-------------------Study glue -------------------*)
(* Explain final gluing pattern *)
{qps,fps} = UnGlue[q,f,salem8 cyclo18 (x-1) (x+1)]
{qst,fst} = UnGlue[q,f,salem8 cyclo6]
{qpt,fpt} = UnGlue[q,f,cyclo18 cyclo6 (x-1) (x+1)]
{qpic,fpic} = {qps,fps}

Print["\n		Gluing pattern\n"]
Print["Glue Salem / periodic: ",GlueGp[qps,fps,salem8]]
Print["Glue Salem / transc  : ",GlueGp[qst,fst,salem8]]
Print["Glue trans / periodic: ",GlueGp[qpt,fpt,cyclo6]]

(*-----------------Ideal theory ----------*)
SetOptions[GroebnerBasis,CoefficientDomain->Integers];
Print[""]
Print["		Ideal Theory"]
Print["Groebner Basis for {salem8, (x+1)^5}  :"]
basis = {salem8, (x+1)^5}
Print[GroebnerBasis[basis]]

Print["Groebner Basis for {u^9-1,2u^4-u^2+2}: "]
basis = {u^9-1,2u^4-u^2+2}
Print[GroebnerBasis[basis]]

Print["Groebner Basis for {t^2-3,t^5}: "]
basis = {t^2-3,t^5}
Print[GroebnerBasis[basis]]
