<< Pari2inMma4.m;

(* This will integrate (1/y)dx and (x/y)dx along the path lst.
   y is understoof as Product[Sqrt[x-rts[[i]]],{i,1,5}] (in particular
   the branch cuts are horisontal to the left of each branch point.
   WARNING: Don't cross them!
   Depending on the optional argument "IntegrationMethod", the following methods
   are available:
     IntegrationMethod -> NIntegrate - use Mathematica's NIntegrate
     IntegrationMethod -> Gauss      - use Gausian Quadrature implemented in 
                                        Mathematica
     IntegrationMethod -> Pari       - use Gausian Quadrature implemented in 
                            Pari, this needs perint7.c to be compiled.
   You can send options to NIntegrate (if this is the method you want to use
   by giving the optional argument "NIntegrateOptions"
*)

Options[gGIntP] = {IntegrationMethod -> Gauss};

gGIntP[rts_,lst_,pres_,opts___] :=  
  Switch[IntegrationMethod /. {opts} /. Options[gGIntP],
    NIntegrate,
      GInt1[rts,lst,pres,opts],
    Gauss,
      GInt2[rts,lst,pres,opts],
    Pari,
      GInt3[rts,lst,pres,opts],
    _,
      GInt1[rts,lst,pres,opts]]

Options[GInt1] = {NIntegrateOptions -> {Compiled -> True}};

GInt1[rts_,lst_,pres_,opts___] := Module[{thef},
  thef[x_] := Evaluate[Product[Sqrt[x - rts[[i]]],{i,1,5}]];
  NIopts = NIntegrateOptions /. {opts} /. Options[GInt1];
  Sum[{Apply[NIntegrate[1/thef[x],{x,lst[[j]],lst[[j+1]]},WorkingPrecision -> pres + 10,#] &,NIopts],Apply[NIntegrate[x/thef[x],{x,lst[[j]],lst[[j+1]]},WorkingPrecision -> pres+10,#]&,NIopts]},{j,1,Length[lst]-1}]]

Options[GInt2] = {GaussN -> 30};
$FirstGInt2 = True;

GInt2[rts_,lst_,pres2_,opts___] := Block[{k,g,st,sum,pres = pres2,n,GPS},
  If[$FirstGInt2,
    Get["~/m/Genus2/GQData.m"];
    $FirstGInt2 = False];
  n = GaussN /. {opts} /. Options[GInt2];
  If[n===30,
    If[pres <= 600,
      GPS = SetPrecision[GP[30],pres+5],
      Print["For Gauss Quadrature with 30 points, Precision must be < 600"];
      Abort[]],
  If[n===90,
    If[pres <= 200,
      GPS = SetPrecision[GP[90],pres+5],
      Print["For Gauss Quadrature with 90 points, Precision must be < 200"];
      Abort[]],
  If[n===150,
    If[pres <= 200,
      GPS = SetPrecision[GP[150],pres+5],
      Print["For Gauss Quadrature with 150 points, Precision must be < 200"];
      Abort[]],
    Print["Only Gaussian Quadrature with 30, 90 or 150 points allowed"];Abort[]]]];
  k = SetPrecision[kaps[n,pres],Round[pres*1.1]];
  g = (Length[rts]-1)/2;
  Print["Integrating ",N[lst]];
  gGIntP2[rts,lst]]

gGIntP2[rts_,lst_] := 
  Module[{a,b},
    If[Length[lst] < 2,0,
     If[lst[[1]] == lst[[2]],gGIntP2[rts,Drop[lst,1]],
      a = lst[[1]];
      b = lst[[2]];
      st = SetPrecision[
        Min[Table[ft1[a,b,rts[[i]],k],{i,1,Length[rts]}]],Round[pres*1.1]];
      sum = 0;
      While[st < 1,
        newp = a + st*(b-a);
        sum = sum + gGInt[rts,a,newp,n];
        a = newp;
        st = SetPrecision[
          Min[Table[ft1[a,b,rts[[i]],k],{i,1,Length[rts]}]],Round[pres*1.1]]];
      sum + gGInt[rts,a,b,n] + gGIntP2[rts,Drop[lst,1]]]]]

(* This will integrate x^i/Product[Sqrt[x-rts[[i]]],{i,1,g}] from c to d, 
   using 
   Gauss Quadrature with n points. Only n = 30, 90 or 150
   will work... 
   For n = 30  precision up to 600 is allowed 
   For n = 90  precision up to 200 is allowed 
   For n = 150 precision up to 200 is allowed *)

gGInt[rts_,c_,d_,n_] := 
 Module[{g,f0,p},
  g = (Length[rts]-1)/2;
  f0[x_] := Evaluate[Product[Sqrt[x-rts[[i]]],{i,1,Length[rts]}]];
  Sum[p = c+GPS[[i,1]]*(d-c);Table[p^j,{j,0,g-1}]*GPS[[i,2]]/f0[p],{i,1,n}]*(d-c)]

ft1[a_,b_,s_,k_] := 
  Module[{c1,c2,c3,qa,qb,qc,tmin,min,tm},
    c1 = Abs[b-a]^2;
    c3 = Abs[a-s]^2;
    c2 = Re[b-a]*Re[a-s]+Im[b-a]*Im[a-s];
    qa = 1-1/k^2;
    qb = 2*c2/c1;
    qc = c3/c1;
    tmin = N[-c2/c1];
    If[tmin < 0,N[k*Abs[a-s]/Abs[b-a]],
      min = c1*tmin^2+2*tmin*c2+c3;
      If[min < 0,min = 0,min = Sqrt[min]];
      tm = N[k*min/Abs[b-a]];
      If[tm > tmin,
        tm,
        N[1/(2*qa)*(-qb-Sqrt[qb^2-4*qa*qc])]]]]

kaps[n_,pres_:60] := N[(10^(-pres)/(Factorial[n+1]^4/((2*n+3)*Factorial[2*n+2]^2)))^(1/(2*n+3))];

GInt3[rts_,lstin_,pres2_,opts___] := Module[{nstr,str,lst},
  lst = {lstin[[1]]};
  Do[If[lstin[[i-1]] == lstin[[i]],Null,AppendTo[lst,lstin[[i]]]],{i,2,Length[lstin]}];
  lst = {Map[MakePariNumber[#,pres2]&,rts],Map[MakePariNumber[#,pres2]&,lst],pres2};
  str = Mma2Pari[lst];
  nstr = StringReplace[str," " -> ""];
  Print["Integrating ",N[lstin]];
  PerIntIt["int:"<>nstr]]

(* This will do the integral to infinity *)
(* Depending on the optional argument "IntegrationMethod", the following 
   methods are available:
     IntegrationMethod -> Gauss      - use Gausian Quadrature implemented in 
                                        Mathematica
     IntegrationMethod -> Pari       - use Gausian Quadrature implemented in 
                            Pari, this
                            needs PariInMath to be loaded, StartPari run and
                            perint5.c compiled into mypari.
*)
Options[GIntInf] = {IntegrationMethod -> Gauss}

GIntInf[rts_,lst_,pres_,opts___] :=  
  Switch[IntegrationMethod /. {opts} /. Options[GIntInf],
    Gauss,
      GIntInf2[rts,lst,pres,opts],
    Pari,
      GIntInf3[rts,lst,pres],
    _,
      GIntInf2[rts,lst,pres,opts]]


GIntInf3[rts_,lstin_,pres2_] := Module[{nstr,str,lst},
  lst = {Map[MakePariNumber[#,pres2]&,rts],Map[MakePariNumber[#,pres2]&,lstin],pres2};
  str = Mma2Pari[lst];
  nstr = StringReplace[str," " -> ""];
  Print["Integrating ",N[lstin]];
  PerIntIt["inf:"<>nstr]]


(* This should at some point lower precision as we go out... *)
GIntInf2[rts_,lstin_,pres2_,opts___] := GInt2[rts,lstin,pres2,opts]