(****************************************************************************

This computes analytic data for a genus 2 hyperelliptic curve given by
y^2 = pol, where pol must be degree 5 or 6. The answer will be
accurate to precision pres. The computed data is:
   {{rts,PerMat},rt6,a,b,{c{int1,int2}}}
where
   rts is the roots (in a specific order) of the quintic. Degree 6 is first
     converted to quintic.
   PerMat is the Period Matrix corresponding to the usuall symplectic basis 
     and the order of the roots as in rts. See Mumford, "Tata lectures on 
     Theta II"
   rt6 is the root that goes to infinity if pol is degree 6. (0 otherwise)
   a is the leading term of pol.
   b is the leading term of the quintic.
   c is a point on the real axis to the right of all rts.
   int1 is the integral from c to +infinity of 
     1/Prod[Sqrt[x-rts[[i]]],{i,1,5}]
   int2 is the integral from c to +infinity of 
     x/Prod[Sqrt[x-rts[[i]]],{i,1,5}]

Possible options are those for gGInt and the option
Root2Infinity. This only makes sense if pol is of degree 6 and then it
should be an exact root of pol. This root is then mapped to infinity
to get the quintic. If pol is degree 6 and this option is omitted a
root will be chosen automatically.

****************************************************************************)
InitCurve[pol_,pres_,opts___] := Module[{dum1,dum2},
  dum1 = GetPerMat[pol,pres,opts];
  dum2 = GetInfInts[dum1,pres,opts];
  Join[dum1,{dum2}]]

(* This get's the period matrix of any genus 2 curve. It returns
   {{rts,PerMat},rt6,a,b}
   where rts is the roots (in a specific order) PerMat is the Period Matrix
   corresponding to the usuall symplectic basis and the order of the roots 
   as in rts. rt6 is the root that goes to infinity if pol is degree 6.
   a is the leading term of pol.
   b is the leading term of the quintic.
*)
Options[GetPerMat] = {Root2Infinity -> Infinity};

GetPerMat[pol_,pres_,opts___] := Module[{art,rt6,rts,a,b},
  {rts,rt6,a,b} = GetRoots3[pol,pres+10,Root2Infinity /. {opts} /. Options[GetPerMat]];
  {GetPerMat2[PerMat[rts,pres,opts],20,pres],rt6,a,b}]

(* This takes the righthandside of a hyperelliptic equation and tries to put 
it in the form y^2 = quintic. It returns the roots to precision pres. It
also returns the root that was sent to infinity to make the sextic quintic *)
GetRoots3[pol_,pres_,art_:Infinity] := Module[{dum,npol,rt,a,b},
  If[Exponent[pol,x] == 5,
    a = 0;
    npol = pol; 
    rt = Infinity,
    a = Coefficient[pol,x,6];
    npol = pol/a;
    If[art === Infinity,
      dum = Factor[npol];
      If[dum[[0]] === Times,
        If[NumberQ[dum[[1]]],dum = dum[[2]],dum = dum[[1]]]];
      rt = (x /. Solve[dum==0,x])[[1]],
      If[FullSimplify[pol /. x -> art] === 0,rt = art,Print["The third argument of GetRoots3 should be an exact root of the sextic."];Abort[]]];
    npol = Expand[Fold[x*#1+#2&,0,Reverse[Map[RootReduce,CoefficientList[(npol /. x -> (1/x + rt))*x^6,x]]]]]];
  b = Coefficient[npol,x,5];
  npol = npol/b;
  Print[npol];
  {HypRoots[npol,pres],rt,a,b}]

HypRoots[f_,pres_] := Module[{},
  deg = Exponent[f,x];
  If[(deg < 1) || (Mod[deg,2] != 1),
    Print["f does not have an appropriate degree"],
    If[Coefficient[f,x,deg] != 1,
      Print["f is not monic"],
      Map[#[[2]] &,Apply[List,NRoots[f == 0,x,pres+10]]]]]]

(* This computes the Period Matrix of an hyperelliptic curve y^2 = f(x) *)
PerMat[irts_,pres_,opts___] := Module[{rts = SetPrecision[irts,pres],dum,Imrts,sas,basepnt,cutlst,aints,pths,htn,vcuts,vc1,vc2,cnt,hc1,hc2,g,nlps,iml,lmat,loops,emat,eform,lstps,lstms},
  dum = Union[Im[rts]];
  Imrts = {dum[[1]]};
  Do[If[dum[[i-1]] == dum[[i]],Null,AppendTo[Imrts,dum[[i]]]],{i,2,Length[dum]}];
  sas = Map[Function[x,Sort[Re[Select[rts,Im[#] == x &]]]],Imrts];
  basepnt = Max[Re[rts]] + 5;
  cutlst = Join[{First[Imrts]-1/2},Table[(Imrts[[i]]+Imrts[[i+1]])/2,{i,1,Length[Imrts]-1}],{Last[Imrts]+1/2}];
  aints = Table[0,{i,1,Length[rts]}];
  pths = Table[0,{i,1,Length[rts]}];
  Do[htn = Position[Imrts,Im[rts[[i]]]][[1,1]];
    vcuts = Join[{First[sas[[htn]]]-1/2},Table[(sas[[htn,j]]+sas[[htn,j+1]])/2,{j,1,Length[sas[[htn]]]-1}],{Last[sas[[htn]]]+1/2}];
    vc1 = Max[Select[vcuts,(# < Re[rts[[i]]]) &]];
    vc2 = Min[Select[vcuts,(# > Re[rts[[i]]]) &]];
    cnt = Count[sas[[htn]],x_/;(x > vc1)];
    hc1 = I*Min[Select[cutlst,(# > Im[rts[[i]]]) &]];
    hc2 = I*Max[Select[cutlst,(# < Im[rts[[i]]]) &]];
    Print["Doing root ",i];
    pths[[i]] =  {basepnt,basepnt+hc1,vc1+hc1,vc1+hc2,vc2+hc2,vc2+hc1,basepnt+hc1,basepnt};
    aints[[i]] = 
      2*gGIntP[rts,{basepnt,basepnt+hc1,vc2+hc1},pres,opts]+
      gGIntP[rts,{vc2+hc1,vc1+hc1,vc1+I*Im[rts[[i]]]},pres,opts]+
      (-1)^cnt*gGIntP[rts,{vc1+I*Im[rts[[i]]],vc1+hc2,vc2+hc2,vc2+I*Im[rts[[i]]]},pres,opts] -
      gGIntP[rts,{vc2+I*Im[rts[[i]]],vc2+hc1},pres,opts],
    {i,1,Length[rts]}];
  g = (Length[rts]-1)/2;
  loops = Table[Join[pths[[i]],Drop[pths[[i+1]],1]],{i,1,2*g}];
  lmat = Transpose[Map[RWNs[#,rts] &,loops]];
  iml = Makeiml[g];
  nlps = Map[LinearSolve[lmat,#] &,iml];
  lstps = Map[Apply[Join,Table[Table[i+(1-Sign[#[[i]]])/2,{j,1,Abs[#[[i]]]}],{i,1,2*g}]] &,nlps];
  lstms = Map[Apply[Join,Table[Table[i+(1+Sign[#[[i]]])/2,{j,1,Abs[#[[i]]]}],{i,1,2*g}]] &,nlps];
  dum = Map[Apply[Plus,Map[Function[x,aints[[x]]],#]] &,lstps] -
        Map[Apply[Plus,Map[Function[x,aints[[x]]],#]] &,lstms];
  emat = Table[Irts[rts,lstps[[i]],lstms[[i]],lstps[[j]],lstms[[j]]],{i,1,2*g},{j,1,2*g}];
  eform[lst1_,lst2_] := (lst1.emat.Transpose[{lst2}])[[1]];
  dum = Transpose[HypSympBase[IdentityMatrix[2*g],eform].dum];
  {irts,aints,pths,dum,tauFromPerM[dum]}]

(* This computes winding numbers around branch points on a
hyperelliptic Riemann surface *)
(* The idea is to count the number and direction in which the path
crosses a branch cuts going to the left of the branch points. Note
that everytime you cross a branch cut your on a different sheet so the
orientation reverses. *)
RWNs[pth_,rts_] := Module[{wn,sgn,mpt},
  Table[
    wn = 0;
    sgn = 1;
    Do[
      If[SCut[pth[[j]],pth[[j+1]],rts[[i]]] != 0,
        mpt = Re[pth[[j]]] + (Im[rts[[i]]]-Im[pth[[j]]])/(Im[pth[[j+1]]]-Im[pth[[j]]])*(Re[pth[[j+1]]]-Re[pth[[j]]]) + Im[rts[[i]]]*I;
        sgn = sgn*(-1)^Sum[SCut[pth[[j]],mpt,rts[[k]],rts[[i]]],{k,1,Length[rts]}];
        wn = wn + sgn;
        sgn = -sgn*(-1)^Sum[SCut[pth[[j+1]],mpt,rts[[k]],rts[[i]]],{k,1,Length[rts]}],
        sgn = sgn*(-1)^Sum[SCut[pth[[j]],pth[[j+1]],rts[[k]]],{k,1,Length[rts]}]],
      {j,1,Length[pth]-1}];
    If[sgn == -1,Print["This is not a closed path"]];
    wn,
    {i,1,Length[rts]}]]

(* This decides whether the line from p1 to p2 crosses the branch cut going horizontally from pt to the left. If Im part of p1 or p2 equalls that of pt then it's considered to be bigger *)
SCut[p1_,p2_,pt_] := 
  If[Im[p1] >= Im[pt],
    If[Im[pt] > Im[p2],
      If[Re[p1] + (Im[pt]-Im[p1])/(Im[p2]-Im[p1])*(Re[p2]-Re[p1]) < Re[pt],
        1,
        0],
      0],
    If[Im[pt] <= Im[p2],
      If[Re[p1] + (Im[pt]-Im[p1])/(Im[p2]-Im[p1])*(Re[p2]-Re[p1]) < Re[pt],
        1,
        0],
      0]]

(* This does the same as above except that now p2's imaginary part is considered to be the same as that of pr *)
SCut[p1_,p2_,pt_,pr_] := 
  If[Im[p1] >= Im[pt],
    If[Im[pt] > Im[p2],
      If[Re[p1] + (Im[pt]-Im[p1])/(Im[p2]-Im[p1])*(Re[p2]-Re[p1]) < Re[pt],
        1,
        0],
      If[Im[pt] == Im[p2],
        If[(Re[p2] < Re[pt]) && (Re[pt] < Re[pr]),
          1,
          0],
        0]],
    If[Im[pt] < Im[p2],
      If[Re[p1] + (Im[pt]-Im[p1])/(Im[p2]-Im[p1])*(Re[p2]-Re[p1]) < Re[pt],
        1,
        0],
      If[Im[pt] == Im[p2],
        If[(Re[p2] < Re[pt]) && (Re[pt] > Re[pr]),
          1,
          0],
        0]]]

(* This makes B's then A's *)
Makeiml[g_] := Module[{iml},
  iml = Table[Table[0,{i,1,2*g+1}],{j,1,2*g}];
  Do[iml[[i+g,2*i-1]] = 1,{i,1,g}];
  Do[iml[[i+g,2*i]] = -1,{i,1,g}];
  Do[Do[iml[[i,j]] = (-1)^j,{j,2*i,2*g+1}],{i,1,g}];
  iml]

(* This assumes that (eform[rbas[[i]],rbas[[j]]])_(i,j) mod 2 is J. It then 
finds a change of variable with a matrix of the form I + 2*A such that
(eform[rbas[[i]],rbas[[j]]])_(i,j) = J  *)
HypSympBase[inbas_,eform_] := Module[{rbas = inbas,n,fl,dum,mnm,k,mink,minj},
  n = Length[rbas];
  If[n == 0,{},
    fl = 1;
    mnm = 0;
    Do[
      Do[
        If[0 < (dum = eform[rbas[[j]],rbas[[k]]]),
          If[(dum < mnm) || (fl == 1),
            mnm = dum;mink = k;minj = j;fl = 0]],
        {k,1,n}],
      {j,1,n}];
    fl = 1;
    Do[
      If[Mod[eform[rbas[[minj]],rbas[[j]]],mnm] != 0,
        rbas[[j]] = rbas[[j]] - 2*Round[eform[rbas[[minj]],rbas[[j]]]/(2*mnm)]*rbas[[mink]];
        fl = 0],
      {j,1,n}];
    If[fl == 0,HypSympBase[rbas,eform],
      Do[
        If[Mod[eform[rbas[[j]],rbas[[mink]]],mnm] != 0,
          rbas[[j]] = rbas[[j]] - 2*Round[eform[rbas[[j]],rbas[[mink]]]/(2*mnm)]*rbas[[minj]];
          fl = 0],
        {j,1,n}];
      If[fl == 0,HypSympBase[rbas,eform],
        If[mnm != 1,Print["This bilinear form does not have det 1"]];
        If[minj < mink,rbas[[minj]] = -rbas[[minj]]];
        Do[
          k = skip[j,minj,mink];
          rbas[[k]] = perp[rbas[[minj]],rbas[[mink]],rbas[[k]],eform],
          {j,1,n-2}];
        dum = HypSympBase[Table[rbas[[skip[j,minj,mink]]],{j,1,n-2}],eform];
        If[minj < mink,
          Join[Take[dum,{1,minj-1}],{rbas[[minj]]},Take[dum,{minj,mink-2}],{rbas[[mink]]},Take[dum,{mink-1,n-2}]],
          Join[Take[dum,{1,mink-1}],{rbas[[mink]]},Take[dum,{mink,minj-2}],{rbas[[minj]]},Take[dum,{minj-1,n-2}]]]]]]]

perp[a0_,a1_,b_,eform_] := (b - (eform[b,a1]/eform[a0,a1])*a0
                  + (eform[b,a0]/eform[a0,a1])*a1);

skip[j_,a0_,a1_] := Module[{a0p,a1p},
  a0p = Min[a0,a1];
  a1p = Max[a0,a1];
  If[j < a0p,
    j,
    If[j+1 < a1p,
      j+1,
      j+2]]]

(* This finds the intersection of two curves built up of pairs of half loops all starting at the same point and going around a single branch point *)
Irts[rts_,lst1p_List,lst1m_List,lst2p_List,lst2m_List] := 
  If[lst1p == {},0,Irts[rts,lst1p[[1]],lst1m[[1]],lst2p,lst2m]+Irts[rts,Drop[lst1p,1],Drop[lst1m,1],lst2p,lst2m]];

Irts[rts_,a_Integer,b_Integer,lst2p_List,lst2m_List] :=
  If[lst2p == {},0,Irts[rts,a,b,lst2p[[1]],lst2m[[1]]]+Irts[rts,a,b,Drop[lst2p,1],Drop[lst2m,1]]];

(* Find the intersection number at the intersection of (a,b) and (c,d) where a,b,c and d are half loops and b follows a and d follows c. *)
Irts[rts_,a_Integer,b_Integer,c_Integer,d_Integer] := Irts[rts,a,c,d] - Irts[rts,b,c,d];

(* This gives 1 if (i,j) splits to the right of (i,k) *)
Irts[rts_,i_,j_,k_] := Irts[rts,i,j]*Irts[rts,i,k]*Irts[rts,j,k];

(* For two unequal half loops, this gives 1 if i is to the right of j, -1 otherwise. If the loops are equall it gives 0. *)
Irts[rts_,i_,j_] := 
  If[i == j,0,
  If[Im[rts[[i]]] == Im[rts[[j]]],
    If[Re[rts[[i]]] < Re[rts[[j]]],1,-1],
    If[Im[rts[[i]]] < Im[rts[[j]]],-1,1]]];


tauFromPerM[PerM_] := Module[{g},
  g = Length[PerM];
  Inverse[Table[Take[PerM[[i]],-g],{i,1,g}]].Table[Take[PerM[[i]],g],{i,1,g}]]

GetPerMat2[pm_,num_,pres_] := Module[{llst,len,dum,sz,bper},
  llst = Permutations[Table[i,{i,1,5}]];
  len = Length[llst];
  dum = Table[{dum = llst[[Random[Integer,{1,len}]]],tsttau[tauFromPerM[PerMatPer[pm,dum][[2]]]]},{k,1,num}];
  sz = Max[Transpose[dum][[2]]];
  bper = dum[[Position[Transpose[dum][[2]],sz][[1,1]]]];
  dum = pm;
  dum[[2]] = SetPrecision[pm[[2]],3*pres];
  dum = PerMatPer[dum,bper[[1]]];
  dum[[2]] = SetPrecision[dum[[2]],pres];
  dum]

(* This takes the output of PerMat and a permutation of the roots and computes tau *)
PerMatPer[{orts_,oaints_,opths_,_,_},perm_] := Module[{dum,aints,pths,g,nlps,iml,lmat,loops,emat,eform,lstps,lstms,rts},
  g = (Length[orts]-1)/2;
  rts = Table[orts[[perm[[i]]]],{i,1,2*g+1}];
  aints = Table[oaints[[perm[[i]]]],{i,1,2*g+1}];
  pths = Table[opths[[perm[[i]]]],{i,1,2*g+1}];
  loops = Table[Join[pths[[i]],Drop[pths[[i+1]],1]],{i,1,2*g}];
  lmat = Transpose[Map[RWNs[#,rts] &,loops]];
  iml = Makeiml[g];
  nlps = Map[LinearSolve[lmat,#] &,iml];
  lstps = Map[Apply[Join,Table[Table[i+(1-Sign[#[[i]]])/2,{j,1,Abs[#[[i]]]}],{i,1,2*g}]] &,nlps];
  lstms = Map[Apply[Join,Table[Table[i+(1+Sign[#[[i]]])/2,{j,1,Abs[#[[i]]]}],{i,1,2*g}]] &,nlps];
  dum = Map[Apply[Plus,Map[Function[x,aints[[x]]],#]] &,lstps] -
        Map[Apply[Plus,Map[Function[x,aints[[x]]],#]] &,lstms];
  emat = Table[Irts[rts,lstps[[i]],lstms[[i]],lstps[[j]],lstms[[j]]],{i,1,2*g},{j,1,2*g}];
  eform[lst1_,lst2_] := (lst1.emat.Transpose[{lst2}])[[1]];
  dum = Transpose[HypSympBase[IdentityMatrix[2*g],eform].dum];
  {rts,dum}]

(* This test tau. If the answer is small theta will take long. 0.25 and g=3 
will take about 40 seconds for a theta evaluation *)
tsttau[tau_] := Min[Re[Eigenvalues[Im[tau]]]];

(* This computes the InfInts for the data given by GetPerMat *)
GetInfInts[gpm_,pres_,opts___] := Module[{rts,basepnt},
  rts = gpm[[1,1]];
  basepnt = Ceiling[Max[Re[rts]] + 5];
  {basepnt,GIntInf[rts,{basepnt,10^(2*(pres+20))},pres+20,opts]}]