(* This applies the {{0,1},{-1,0}} to an element of 2-dimensional Siegel 
upperhalf space.  *)
applyS[z_] := 
  applyG[z,{{0,0,1,0},{0,0,0,1},{-1,0,0,0},{0,-1,0,0}}];

(* This applies {{1,b},{0,1}} to an element of 2-dimensional Siegel 
upperhalf space.  *)
applyT[z_,b_] :=
  applyG[z,{{1,0,b[[1,1]],b[[1,2]]},{0,1,b[[2,1]],b[[2,2]]},
                 {0,0,1,0},              {0,0,0,1}}];

applyU[z_,u_] := Module[{ui},
  ui = Transpose[Inverse[u]]; 
  applyG[z,{{u[[1,1]],u[[1,2]],0,0},{u[[2,1]],u[[2,2]],0,0},
          {0,0,ui[[1,1]],ui[[1,2]]},{0,0,ui[[2,1]],ui[[2,2]]}}]];
 
applyG[z_,g_] := Module[{totalChange = z[[2]],nz},
  If[SympQ[g],Null,Print["Ooops...",g," is not a symplectic matrix."]];
  totalChange = g.totalChange;
  nz = ({{g[[1,1]],g[[1,2]]},{g[[2,1]],g[[2,2]]}}.z[[1]] + {{g[[1,3]],g[[1,4]]},{g[[2,3]],g[[2,4]]}}).
  Inverse[{{g[[3,1]],g[[3,2]]},{g[[4,1]],g[[4,2]]}}.z[[1]] + {{g[[3,3]],g[[3,4]]},{g[[4,3]],g[[4,4]]}}];
  {nz,totalChange}];

aG[z_, m_] := applyG[{z, IdentityMatrix[4]}, m][[1]];

(* This takes a symplectic matrix and outputs the vector {a,b,c,d} *)
abcd[g_] := {{{g[[1,1]],g[[1,2]]},{g[[2,1]],g[[2,2]]}},{{g[[1,3]],g[[1,4]]},{g[[2,3]],g[[2,4]]}},{{g[[3,1]],g[[3,2]]},{g[[4,1]],g[[4,2]]}},{{g[[3,3]],g[[3,4]]},{g[[4,3]],g[[4,4]]}}};

(* This find a random 2 by 2 unimodular integer matrix *)
ranunimod[n_] := Module[{a,b,c,d,k},
  a = Random[Integer,{-n,n}];
  b = Random[Integer,{-n,n}];
  While[GCD[a,b] =!= 1,a = Random[Integer,{-n,n}];b = Random[Integer,{-n,n}]];
  If[b == 0,{{a,0},{Random[Integer,{-n,n}],{-1,1}[[Random[Integer,{1,2}]]]}},
    If[a == 0,{{0,b},{{-1,1}[[Random[Integer,{1,2}]]],Random[Integer,{-n,n}]}},
      If[Mod[a,b] == 0,c = b*(a-1);d=1,
        d = PowerMod[a,-1,b];
        c = (a*d-1)/b];
      k = Random[Integer,{Max[Ceiling[-n/Abs[b]-d/b],Ceiling[-n/Abs[a]-c/a]],
                          Min[Floor[n/Abs[b]-d/b],Floor[n/Abs[a]-c/a]]}];
      {{a,b},{c+k*a,d+k*b}}]]]

(* This translates an element of 2-dimensional Siegel upperhalf space into 
the strip Abs[Re[z]] <= 1/2.  *)
TransStrip[z_] := applyT[z,
  -{{Round[Re[z[[1,1,1]]]],Round[Re[z[[1,1,2]]]]},
    {Round[Re[z[[1,1,2]]]],Round[Re[z[[1,2,2]]]]}}];

(* This makes a random 2 by 2 symmetric matrix *)
ransym[n_] := Module[{b},
  {{Random[Integer,{-n,n}],b = Random[Integer,{-n,n}]},
   {b                     ,Random[Integer,{-n,n}]}}]

(* This makes a semi random element of the symplectic group *)
(* n is the number of less random symplectic matrices being multiplied
   and m is the size of the coefficients of these matrices *)
ransymp[n_:10,m_:2] := Module[{out,u,ui,b},
  out = IdentityMatrix[4];
  Do[d = Random[Integer,{1,3}];
    If[d == 1,
      u = ranunimod[m];
      ui = Transpose[Inverse[u]]; 
      out = {{u[[1,1]],u[[1,2]],0,0},{u[[2,1]],u[[2,2]],0,0},
             {0,0,ui[[1,1]],ui[[1,2]]},{0,0,ui[[2,1]],ui[[2,2]]}}.out,
      If[d == 2,
        b = ransym[m];
        out = {{1,0,b[[1,1]],b[[1,2]]},{0,1,b[[2,1]],b[[2,2]]},
                 {0,0,1,0},              {0,0,0,1}}.out,
        out = {{0,0,1,0},{0,0,0,1},{-1,0,0,0},{0,-1,0,0}}.out]];
    ,{i,1,n}];
  out]


(* This writes a positive definite symetric matrix as ^tB.B  *)
GetB[q_] := Module[{dum,l1,l2,b},
  dum = Eigensystem[q];
  l1 = Sqrt[dum[[2,1]].dum[[2,1]]];
  l2 = Sqrt[dum[[2,2]].dum[[2,2]]];
  b = {{Sqrt[dum[[1,1]]],0},{0,Sqrt[dum[[1,2]]]}}.{dum[[2,1]]/l1,dum[[2,2]]/l2}];

(* Given a positive definite symmetric matrix q this will find integers a and 
b such that {a,b}.q.{{a},{b}} is non-zero minimal *)
(* It uses GetB and then Algorithm 1.3.14 in Cohen... *)
FindShort[q_] := Module[{dum,bb,a,b,sa,sb,n,r,st,t},
  bb = GetB[q];
  a = {bb[[1,1]],bb[[2,1]]};
  b = {bb[[1,2]],bb[[2,2]]};
  sa = a.a; sb = b.b;
  If[sa < sb,dum = a;a = b;b = dum;dum = sb;sb = sa;sa = dum];
  n = a.b; r = Round[n/sb]; st = sa-2*r*n+r^2*sb;
  While[st < sb,
    t = a - r*b; a = b; b = t; sa = sb; sb = st;
    n = a.b; r = Round[n/sb]; st = sa-2*r*n+r^2*sb];
  Round[LinearSolve[bb,b]]];

(* The point here is to find a unimodular integral matrix {{a,c},{b,d}} such
that Transpose[{{a},{b}}].y.{{a},{b}} is minimal and 
     Transpose[{{c},{d}}].y.{{c},{d}} is minimal (given a and b) and
Transpose[{{c},{d}}].y.{{a},{b}} >= 0. *)
MinkowRed[{{al_,ga_},{_,be_}}] := 
  Module[{a,b,c,d,e,oy,outm},
  {a,b} = FindShort[{{al,ga},{ga,be}}];
  If[a == 0,c = -b;d = 0,
    If[b==0,d = a;c = 0,d = PowerMod[a,-1,b];c = (a*d-1)/b]];
  e = Round[-(a*al*c + b*be*d + b*c*ga + a*d*ga)/
            (a^2*al + b^2*be + 2*a*b*ga)];
  oy = Sign[(d + b*e)*(b*be + a*ga) + (c + a*e)*(a*al + b*ga)];
  If[oy == 0,oy = 1];
  outm = {{a,oy*(a*e+c)},{b,oy*(b*e+d)}};
  {Transpose[outm].{{al,ga},{ga,be}}.outm,outm}]

RanSPmat[n_] := Module[{al=0.5/n,ga=1},
  While[ga^2/al > n,al = Random[Real,{0,n}];ga = Random[Real,{-n,n}]];
  {{al,ga},{ga,Random[Real,{ga^2/al,n}]}}]

RanSiegel[n_] := Module[{dum},
  {{Random[Real,{-n,n}],dum = Random[Real,{-n,n}]},
   {   dum     ,Random[Real,{-n,n}]}} + I*RanSPmat[n]];

tstMats = {{{0,0},{0,0}},{{1,0},{0,0}},{{0,0},{0,1}},{{1,0},{0,1}},{{1,0},{0,-1}},{{0,1},{1,0}},{{1,1},{1,0}},{{0,1},{1,1}},{{0,0},{0,0}}}

tstMts2 = {{{0, 0, 0, 1},  {1, 0, 1, 0}, {0, -1, 0, 0},  {0, 0, 1, 0}}, 
           {{0, 0, -1, 0}, {0, 1, 0, 1}, {1, 0, 0, 0},   {0, 0, 0, 1}}, 
           {{0, 0, -1, 0}, {0, 1, 1, 1}, {1, -1, 1, 0},  {0, 0, 1, 1}}, 
           {{0, 0, -1, 0}, {0, 1, 1, 1}, {1, -1, -1, 0}, {0, 0, 1, 1}}}

Fundom[z_] := Module[{nz,fl,k},
  If[Dimensions[z] === {2},
    nz = z,
    nz = {z,IdentityMatrix[4]}];
  fl = True;
  While[fl,
    nz = TransStrip[applyU[nz,Transpose[MinkowRed[Im[nz[[1]]]][[2]]]]];
    k=1;
    While[Abs[Det[nz[[1]]+tstMats[[k]]]] >= 1 && k<9,k++];
    If[k == 9,
      k=1;
      While[Abs[Det[nz[[1]]-tstMats[[k]]]] >= 1 && k<9,k++];
      If[k==9,
        If[Abs[nz[[1,1,1]]] < 1,
          nz = applyG[nz,tstMts2[[2]]],
          If[Abs[nz[[1,2,2]]] < 1,
            nz = applyG[nz,tstMts2[[3]]],
            If[Abs[nz[[1,1,1]]+nz[[1,2,2]]-2*nz[[1,1,2]]+1] < 1,
              nz = applyG[nz,tstMts2[[3]]],
              If[Abs[nz[[1,1,1]]+nz[[1,2,2]]-2*nz[[1,1,2]]-1] < 1,
                nz = applyG[nz,tstMts2[[4]]],
                fl = False]]]],
        nz = applyG[nz,{{0,0,-1,0},{0,0,0,-1},{1,0,-tstMats[[k,1,1]],-tstMats[[k,1,2]]},{0,1,-tstMats[[k,2,1]],-tstMats[[k,2,2]]}}]],
(*        Print["This is where the mistake was"]], *)
      nz = applyG[nz,{{0,0,-1,0},{0,0,0,-1},{1,0,tstMats[[k,1,1]],tstMats[[k,1,2]]},{0,1,tstMats[[k,2,1]],tstMats[[k,2,2]]}}]]];
  If[Precision[nz] == 16,Print["Precision is only 16..."]];
  nz]

SympQ[{{a11_,a12_,b11_,b12_},
       {a21_,a22_,b21_,b22_},
       {c11_,c12_,d11_,d12_},
       {c21_,c22_,d21_,d22_}}] := Apply[And,
  {a12*c11 - a11*c12 + a22*c21 - a21*c22 == 0, 
   -(a12*c11) + a11*c12 - a22*c21 + a21*c22 == 0, 
   -1 - b11*c11 - b21*c21 + a11*d11 + a21*d21 == 0, 
   -(b11*c12) - b21*c22 + a12*d11 + a22*d21 == 0, 
   -(b12*c11) - b22*c21 + a11*d12 + a21*d22 == 0, 
   -1 - b12*c12 - b22*c22 + a12*d12 + a22*d22 == 0, 
   b12*d11 - b11*d12 + b22*d21 - b21*d22 == 0}]


(* This takes two taus and if they are equivalent under Symp(4,Z) it returns
a symplectic matrix such that aG[tau1,smat] = tau2. *)
(* Well it should, but: this only checks whether an integer translation and 
one of the tstMts2 takes Fundom[tau1] to Fundom[tau2]. Of course the ultimate
solution would be to make Fundom precision proof and not include all 
boundaries in the fundamental domain... *)
EquivQold[tau1_,tau2_] := Module[{ftau1,ftau2},
  ftau1 = Fundom[tau1];
  ftau2 = Fundom[tau2];
  EquivQ1[ftau1,ftau2]]

EquivQ1[ftau1_,ftau2_] := Module[{dum,k},
  dum = equivqsubroutine[ftau1,ftau2];
  k=1;
  While[Not[dum[[1]]] && (k < 5),
    dum = equivqsubroutine[applyG[ftau1,tstMts2[[k]]],ftau2];
    k = k+1];
  If[k == 5,False,Inverse[ftau2[[2]]].dum[[2,2]]]]

equivqsubroutine[tau1_,tau2_] := Module[{dum1,dum2},
  dum1 = tau1[[1]]-tau2[[1]];
  dum2 = Map[Round,dum1,{2}];
  If[Max[Map[Abs,Flatten[dum2-dum1]]] > 10^(-10),{False,Null},{True,applyT[tau1,-dum2]}]]

(* The above does not work:
It claims the two matrices:
  {{-2*x - x^3/2, 3 + x^2}, { 3 + x^2, x^3/4}} 
and
  {{-2*x - x^3/2,-3 - x^2}, {-3 - x^2, x^3/4}} 
where x is -Sqrt[-5 + Sqrt[5]] are inequivalent, but 
  {{1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, -1}}
takes the one to the other...
Both these matrices seem to be in Gottschling's fundamental domain...

So as an ugly fix, try *)
EquivQ[tau1_,tau2_] := Module[{ftau1,ftau2,dum,k},
  ftau1 = Fundom[tau1];
  ftau2 = Fundom[tau2];
  If[(dum = EquivQ1[ftau1,ftau2]) =!= False,dum,
    k = 1;
    While[((dum = EquivQ1[Fundom[applyG[ftau1,ransymp[]]],ftau2]) === False) && (k<10),
      k = k + 1];
    If[k==10,False,dum]]]