<< Pari2inMma4.m;

Options[theta] = {ThetaMethod -> NSum};

theta[{cp_,cpp_},z_,ttau_,pres_,opts___] := 
Module[{},
  Switch[ThetaMethod /. {opts} /. Options[theta],
    Pari,
      If[z === {{0},{0}},
        theta[{cp,cpp},z,ttau,pres,opts] = rtheta[{cp,cpp},z,ttau,pres],
        rtheta[{cp,cpp},z,ttau,pres]],
    _,
      If[z === {{0},{0}},
        theta[{cp,cpp},z,ttau,pres,opts] = ntheta[{cp,cpp},z,ttau,pres],
        ntheta[{cp,cpp},z,ttau,pres]]]]

rtheta[{cp_,cpp_},z_,ttau_,pres_:16] := Module[{str},
  str = "th:[" <> Mma2Pari[Join[cp,cpp]] <> ",";
  str = str <> Mma2Pari[Map[MakePariNumber[#,pres+10] &,Transpose[z][[1]]]]<>"~" <> ",";
  str = str <> ToPariMatrix[Map[MakePariNumber[#,pres+10] &,ttau,{2}]] <> "," <> ToString[pres] <> "]";
  str = StringReplace[str," " -> ""];
  PerIntIt[str]]

ntheta[{cp_,cpp_},z_,ttau_,pres_:16] := 
 Module[{nz,v1,l1,mz,mn,tau = ttau,dum},
  nz = z + ttau.Transpose[{cp}] + Transpose[{cpp}];
  v1 = Floor[Inverse[Im[tau]].Im[nz]];
  nz = nz - tau.v1;
  l1 = N[Pi*Min[Re[Eigenvalues[Im[ttau]]]]];
  mz = 2*Pi*Max[Abs[Im[Transpose[nz][[1]]]]];
  mn = Ceiling[N[Sqrt[pres*Log[10.]/l1 + mz*mz/(4*l1*l1)]+mz/(2*l1)]];
  N[Exp[Pi*I*(cp.ttau.Transpose[{cp}] + 2*cp.(z+Transpose[{cpp}]))],pres][[1]]*
  N[Exp[Pi*I*(-Transpose[v1].tau.v1-2*Transpose[v1].nz)],pres][[1,1]]*
  Sum[N[Exp[Pi*I*({m1,m2}.tau.{{m1},{m2}}+2*{m1,m2}.nz)[[1]]],pres],{m1,-mn,mn},{m2,-mn,mn}]];