(***********************************************************************)
(*                                                                     *)
(*                PROGRAM:        INSTANTON    (1.0)                   *)
(*                                                                     *)
(*                   send suggestions to A. Klemm                      *)
(*                        klemm@nxth21.cern.ch                         *)
(*                                                                     *)
(***********************************************************************)

(* This a Mathematica program, which calculates the instanton expansion
   for hypersurfaces and complete intersetions in toric varieties, as
   described in the papers S. Hosono, A Klemm, S.Theisen and S.T. Yau:
   [1] "Mirror Symmetry, Mirror Map and Applications to Calabi-Yau
   Hypersurfaces" hep-th/9308122, to be published in Comm. Math. Phys.
   [2] "Mirror Symmetry, Mirror Map and Applications to Complete
   Intersection Calabi-Yau Spaces" CERN-TH.7303/June 1994.

   Strip it from the above TeX file and called it e.g. "inst.m".


   It works in two modes for the cases:

   a.) Hypersurfaces and complete intersections in non singular
       (products of weighted) projective spaces, such as
       the Quintic or the Tian-Yau manifold. Load the program in a
       Mathematica session ( <<inst.m ), invoke the
       function "casea" and follow interactively the instructions.
       If you wish to bypass this input routine specify the
       l^{(i)} (comp ref. [2]) in listform e.g. for the Tian-Yau case.

       $l={{{-3, 0,-1},{1,1,1,1,0,0,0,0}},
           {{ 0,-3,-1},{0,0,0,0,1,1,1,1}}}

       and type "run".


   b.) Singular toric varieties, such as the hypersurface of
       degree eight in IP(2,2,2,1,1). Here you have to specify
       the generators of the Mori cone l^{(i)} in the form, which
       e.g. for the example above reads (comp. ref. [1])

       lm[number of model]={{{-4},{1, 1, 1, 0, 0, 1}},
                           {{ 0},{0, 0, 0, 1, 1,-2}}}

       and the principle parts of the Picard-Fuchs equation see ref.
       [1],[2], in the example

       orth[number of model]={o[1]^2 (o[1]-2*o[2]),
                              o[2]^2}

       and then you invoke "caseb[number of model]". You will be asked
       for the normalisation for J^3, which you can find in ref [1].

       There are several examples of this type included below.

       For both versions there will be several intermediate
       data available after running the program. I.e. the
       pure powerseries part of the periods will be stored:
       S[0] is the pure powerseries solution, S[i],....,S[dim Modulispace]
       is the powerseries part for the single logarithmic solutions,
       S[i,j] i<j are the ones of the double logarithmic solutions.
       The mirror map z_i(q) is stored in qs[i],
       its inversion in the rule rzq. The Yukawa couplings are named
       nyk[i,j,k]. The instanton numbers will
       be written in the file "inst.out".


       Runtime examples: Two moduli cases up to total degree six in the
                         instanton expansion will need few minutes
                         on a (1994) workstation. Five moduli examples up to
                         degree eight can take some hours. You can easily
                         crash your system by demanding too high degrees.



**************************************************************************)


(*                 DATA  FOR   THE  EXAMPLES                             *)


(*Example 1: Degree eight hypersurface in IP(2,2,2,1,1)         run caseb[1]
   *)



lm[1]={{{-4},{1, 1, 1, 0, 0, 1}},
      {{ 0},{0, 0, 0, 1, 1,-2}}}

orth[1]={o[1]^2 (o[1]-2*o[2]),
      o[2]^2}


(*Example 2: Complete intersection X_{4,6}(2,2,2,2,1,1)    run caseb[2]  *)

lm[2]={{{-2,-3},{1, 1, 1, 1, 0, 0, 1}},
       {{ 0, 0},{0, 0, 0, 0, 1, 1,-2}}}

orth[2]={o[1]^2 (o[1]- 2 o[2]),
         o[2]^2}



(*Example 3: Complete intersection X_{4,4,4}(2,2,2,2,2,1,1)  run caseb[3] *)

lm[3]={{{-2,-2,-2},{1, 1, 1, 1, 1, 0, 0, 1}},
       {{ 0, 0, 0},{0, 0, 0, 0, 0, 1, 1,-2}}}

orth[3]={o[1]^2 (o[1] - 2 o[2]),
         o[2]^2}



(*Example 4: Degree fourteen hypersurface in IP(7,2,2,2,1)     run caseb[4]
    *)


lm[4]={{{-7},{0, 1, 1, 1, -3, 7}},
      {{ 0}, {1, 0, 0, 0, 1,-2}}}

orth[4]={o[1]^2 (7*o[1]-2*o[2]),
         o[2]*(o[2]-3*o[1]) }




(*Example 5: Degree twelve hypersurface in IP(4,3,2,2,1)      run caseb[5]
 *)


lm[5]={{{-6},{2, 0, 1, 1,-1, 3}},
      {{ 0}, {0, 1, 0, 0, 1,-2}}}

orth[5]={o[1]^2 (3*o[1]-2*o[2]),
         o[2]*(o[2]-o[1])}






(*Example 6: Degree eighteen hypersurface in IP(9,6,1,1,1)      run caseb[6]
  *)


lm[6]={{{-6},{3, 2, 0, 0, 0, 1}},
      {{ 0}, {0, 0, 1, 1, 1,-3}}}

orth[6]={o[1] (o[1]-3*o[2]),
         o[2]^3 }


(*Example 11-15: Degree eighteen hypersurface in IP(6,6,3,2,1)             *)


(* subdivison A:  run caseb[11] *)
lm[11]={{{-3},{ 1, 1, 0, 0,-1, 1, 0, 1, 0}},
    {{ 0},{ 0, 0, 1, 0, 0,-1, 0, 1,-1}},
    {{ 0},{ 0, 0, 0, 0, 1,-1, 0,-1, 1}},
    {{ 0},{ 0, 0, 0, 0, 0, 1, 1,-1,-1}},
    {{ 0},{ 0, 0, 0, 1, 0, 0,-2, 1, 0}}}


orth[11]={o[1]*(o[1]-o[2]-o[3]+o[4])*(o[1]+o[2]-o[3]-o[4]+o[5]),
      o[1]*(o[2]-o[3]+o[4]),
      o[5]*(o[1]+o[2]-o[3]-o[4]+o[5]),
      (o[1]-o[2]-o[3]+o[4])*(o[4]-2*o[5]),
      o[2]*(o[1]+o[2]-o[3]-o[4]+o[5]),
      (o[1]-o[3])*(o[2]-o[3]+o[4]),
      (o[1]-o[3])*o[2],
      (o[1]-o[3])*o[5],
      o[5]*(o[1]-o[2]-o[3]+o[4]),
      (o[1]-o[3])*(o[4]-2*o[5]),
      o[2]*(o[4]-2*o[5])}

(* subdivision B: run caseb[12] *)
lm[12]={{{-3},{ 1, 1, 0, 0,-1, 1, 0, 1, 0}},
    {{ 0},{ 0, 0,-1, 0, 0, 1, 0,-1, 1}},
    {{ 0},{ 0, 0, 1, 0, 1,-2, 0, 0, 0}},
    {{ 0},{ 0, 0, 1, 0, 0, 0, 1, 0,-2}},
    {{ 0},{ 0, 0, 0, 1, 0, 0,-2, 1, 0}}} ;


orth[12]={o[1]*(o[1]+o[2]-2*o[3])*(o[1]-o[2]+o[5]),
      o[1]*(-o[2]+2*o[4]),
      o[5]*(o[1]-o[2]+o[5]),
      (o[1]+o[2]-2*o[3])*(o[4]-2*o[5]),
      (-o[2]+o[3]+o[4])*(o[1]-o[2]+o[4]),
      (o[1]-o[3])*(-o[2]+2*o[4]),
      (o[1]-o[3])*(-o[2]+o[3]+o[4]),
      (o[1]-o[3])*o[5],
      o[5]*(o[1]+o[2]-2*o[3]),
      (o[1]-o[3])*(o[4]-2*o[5]),
      (-o[2]+o[3]+o[4])*(o[4]-2*o[5]),
      (***** addition ****)
      o[1]*(-o[2] + o[3] + o[4])*(o[1] - o[2] + o[5])}



(* subdivision C: run caseb[13] *)
lm[13]={{{-3},{ 1, 1, 0, 0, 0, 0, 0, 0, 1}},
     {{ 0},{ 0, 0, 1, 0, 1,-2, 0, 0, 0}},
     {{ 0},{ 0, 0, 0, 0,-1, 1, 0, 1,-1}},
     {{ 0},{ 0, 0, 0, 0, 1, 0, 1,-2, 0}},
     {{ 0},{ 0, 0, 0, 1, 0, 0,-2, 1, 0}}} ;

orth[13]={o[1]*(-2*o[2]+o[3])*(o[3]-2*o[4]+o[5]),
      o[1]*(-o[1]+o[3]),
      o[5]*(o[3]-2*o[4]+o[5]),
      (-2*o[2]+o[3])*(o[4]-2*o[5]),
      o[2]*(o[3]-2*o[4]+o[5]),
      (-2*o[2]+o[3])*(o[3]-2*o[4]+o[5]),
      (-o[2]+o[3]-o[4])*o[2],
      (-o[2]+o[3]-o[4])*o[5],
      o[5]*(-2*o[2]+o[3]),
      (-o[2]+o[3]-o[4])*(o[4]-2*o[5]),
      o[2]*(o[4]-2*o[5])}

(* subdivision D:  run caseb[14] *)
lm[14]={{{-3},{ 1, 1, 0, 0,-1, 1, 0, 1, 0}},
    {{ 0},{ 0, 0, 1, 0, 0, 0, 1, 0,-2}},
    {{ 0},{ 0, 0, 0, 0, 1, 0, 1,-2, 0}},
    {{ 0},{ 0, 0, 0, 0, 0,-1,-1, 1, 1}},
    {{ 0},{ 0, 0, 0, 1, 0, 1,-1, 0,-1}}} ;

orth[14]={o[1]*(o[1]-o[4]+o[5])*(o[1]-2*o[3]+o[4]),
      o[1]*(2*o[2]-o[4]+o[5]),
      o[5]*(o[1]-2*o[3]+o[4]),
      (o[1]-2*o[3]+o[4])*(2*o[2]-o[4]+o[5]),
      o[2]*(o[1]-2*o[3]+o[4]),
      (o[1]-o[3])*(2*o[2]-o[4]+o[5]),
      (o[1]-o[3])*o[2],
      (o[1]-o[3])*o[5],
      o[5]*(o[1]-o[4]+o[5]),
      (o[1]-o[3])*(o[2]+o[3]-o[4]-o[5]),
      o[2]*(o[2]+o[3]-o[4]-o[5]),
      (**** addition ****)
      o[1]*(o[2] + o[3] - o[4] - o[5])*(o[1] - o[4] + o[5])}

(* subdivision E: run caseb[15] *)
lm[15]={{{-3},{ 1, 1, 0, 0,-1, 1, 0, 1, 0}},
    {{ 0},{ 0, 0, 1, 1, 0, 1, 0, 0,-3}},
    {{ 0},{ 0, 0, 0, 0, 1, 0, 1,-2, 0}},
    {{ 0},{ 0, 0, 0, 1, 0, 0,-2, 1, 0}},
    {{ 0},{ 0, 0, 0,-1, 0,-1, 1, 0, 1}}} ;

orth[15]={o[1]*(o[1]+o[2]-o[5])*(o[1]-2*o[3]+o[4]),
      o[1]*(3*o[2]-o[5]),
      (o[2]+o[4]-o[5])*(o[1]-2*o[3]+o[4]),
      (o[1]-2*o[3]+o[4])*(3*o[2]-o[5]),
      o[2]*(o[1]-2*o[3]+o[4]),
      (o[1]-o[3])*(3*o[2]-o[5]),
      (o[1]-o[3])*o[2],
      (o[1]-o[3])*(o[2]+o[4]-o[5]),
      (o[3]-2*o[4]+o[5])*(3*o[2]-o[5]),
      (o[1]-o[3])*(o[3]-2*o[4]+o[5]),
      o[2]*(o[3]-2*o[4]+o[5]),
      (*** addition ***)
      o[2]*(o[2]+o[4]-o[5])*(o[1]+o[2]-o[5]),
      o[1]*(o[1]+o[2]-o[5])*(o[3]-2*o[4]+o[5]),
      o[1]*(o[2]+o[4]-o[5])*(o[1]+o[2]-o[5])}






(***************************** PROGRAM ********************************)


(**************************INPUT ROUTINES******************************)

casea := Block[{npol,nfac,wi,nweights},
      Print["Dimension of the manifold"];
      ndim=Input[];
      Print["Specify the number of factors of projective spaces"];
      nfac=Input[];
      nweights=0;
      offset[1]=0;
      Do[
         Print["Specify the dimension of the ",i," th  projective space"];
         dimension[i]=Input[];
         nweights=nweights+dimension[i]+1;
         offset[i+1]=nweights;
         j=1;
         wi=1;
         Print["Specify the weights"];
         Print["(If you input zero ordinary projective is assumed)"];
         While[j <= dimension[i]+1 && wi =!= 0,
           Print[j," th weight of the ",i," th space"];
           wi=Input[];
           If[wi == 0,Do[weight[i,k]=1,{k,1,dimension[i]+1}],weight[i,j]=wi];
           j++],{i,1,nfac}];
         npol=nweights-nfac-ndim;
         Print["Specify the degrees of the ",npol," polynomials"];
         Do[
         Print["Degree of the ",i," th polynomial in the coordinates of the ",
              j," th projective space"];
         degrees[i,j]=Input[],{i,1,npol},{j,1,nfac}];
         $l={};
         Do[
           hl1=Table[0,{ii,1,nweights}];
           Do[hl1[[kk+offset[i]]]=weight[i,kk],{kk,1,dimension[i]+1}];
           hl0=Table[-degrees[kk,i],{kk,1,npol}];
           hl=Join[{hl0},{hl1}];
           $l=Join[$l,{hl}],{i,1,nfac}];run]


caseb[n_] := Block[{hi},
       $l=lm[n];
       $orth=orth[n];
       $m=Length[$l];
       $Js=Table[J[i],{i,1,$m}];
       $qs=Table[q[i],{i,1,$m}];
       findring[$orth];
       Print["Expansion Precision, the highest (multi)degree evaluated"];
       $no=Input[];
       Print["Evaluating the Period"];
       pals;
       yukawa;
       Write["inst.out","*********************************************"];
       Save["inst.out",$l];
       Write["inst.out","*********************************************"];
       Save["inst.out",ring];
       Write["inst.out","*********************************************"];
       Save["inst.out",nrc]]


run := Block[{hi},
       $m=Length[$l];
       $Js=Table[J[i],{i,1,$m}];
       $qs=Table[q[i],{i,1,$m}];
       Print["Expansion Precision, the highest (multi)degree  evaluated"];
       $no=Input[];
       Print["Evaluating the Period"];
       pals;
       topdata;
       yukawa;
       Write["inst.out","*********************************************"];
       Save["inst.out",$l];
       Write["inst.out","*********************************************"];
       Save["inst.out",ring];
       Write["inst.out","*********************************************"];
       Save["inst.out",nrc]]




(****************** CONSTRUCTION OF THE SOLUTIONS**************)

pals/: pals := Block[{hi},
   Clear[S];
   li0=-Map[Drop[#1,-1][[1]]&,$l];
   li1=Map[Drop[#1,1][[1]]&,$l];
   var=Table[z[i],{i,1,$m}];
   lpw=Map[#1-1&,Flatten[Array[a,Table[$no+1,{i,1,$m}]]]/.{a->List}];
   lpw=Select[lpw,Plus @@ #1 <= $no&];
   Do[S[i]=0;Do[S[i,j]=0,{j,1,$m}],{i,0,$m}];
   Do[
     lh=lpw[[k]];
     mon=Inner[Power,var,lh,Times]/.z[u_]:>h*z[u];
     (*Print[mon];*)
     nl0=Inner[Times,lh,li0];
     nl1=Inner[Times,lh,li1];
     (*Print[nl1];*)
     np=Flatten[Position[Negative[nl1],True]];
     (*Print[np];*)
     case=Length[np];
     If[case == 0,
       w0=Apply[Times,Map[#!&,nl0]]/Apply[Times,Map[#!&,nl1]];
       S[0]=S[0]+w0*mon;
       pg0=Map[pg[#]&,nl0];
       pg1=Map[pg[#]&,nl1];
       dpg0=Map[dpg[#]&,nl0];
       dpg1=Map[dpg[#]&,nl1];
       Do[
         Ab=Inner[Times,li0[[b]],pg0]-Inner[Times,li1[[b]],pg1];
         S[b]=S[b]+w0*Ab*mon;
         Do[
           Ac=Inner[Times,li0[[c]],pg0]-Inner[Times,li1[[c]],pg1];
           Bbc=Inner[Times,Inner[Times,li0[[b]],li0[[c]],List],dpg0]-
               Inner[Times,Inner[Times,li1[[b]],li1[[c]],List],dpg1];
            S[c,b]=S[c,b] + w0*(Ac*Ab+Bbc)*mon,
         {c,1,b}],
        {b,1,$m}],
     If[case == 1,
        kk=np[[1]];
        mk=-nl1[[kk]];
        nl1=Select[nl1,#1 >= 0&];
        pg0=Map[pg[#]&,nl0];
        pg1=Map[pg[#]&,nl1];
        pgmk=pg[mk-1];
        Do[
           A2=(-1)^(mk+1)*Apply[Times,Map[#!&,nl0]]*(mk-1)!/
              Apply[Times,Map[#!&,nl1]];
           S[b]=S[b]+li1[[b]][[kk]]*A2*mon;
           Ab=Inner[Times,li0[[b]],pg0]-Inner[Times,Drop[li1[[b]],{kk}],pg1]-
              li1[[b]][[kk]]*pgmk;
           Do[
             Ac=Inner[Times,li0[[c]],pg0]-Inner[Times,Drop[li1[[c]],{kk}],pg1]-
             li1[[c]][[kk]]*pgmk;
             S[c,b]=S[c,b]+A2*(Ab*li1[[c]][[kk]]+Ac*li1[[b]][[kk]])*mon,
           {c,1,b}],
        {b,1,$m}]
      ,
      If[case == 2,
        k1=np[[1]];
        k2=np[[2]];
        mk1=-nl1[[k1]];
        mk2=-nl1[[k2]];
        nl1=Select[nl1,#1 >= 0&];
        A3=(-1)^(mk1+mk2)*Apply[Times,Map[#!&,nl0]]*(mk1-1)!*(mk2-1)!/
           Apply[Times,Map[#!&,nl1]];
        Do[
           S[c,b]=S[c,b]+A3*(li1[[c]][[k1]]*li1[[b]][[k2]]+
                             li1[[b]][[k1]]*li1[[c]][[k2]])*mon,
         {b,1,$m},{c,1,b}]
     ,],],],
    {k,1,Length[lpw]}];
    Do[
      qs[i]=Normal[sexpand[Series[h z[i] Exp[S[i]/S[0]],{h,0,$no}]]],
    {i,1,$m}]]

pg[a_] := PolyGamma[a+1]/.EulerGamma->0
dpg[a_]:= PolyGamma[1,a+1]/.Pi->0;




(*********************TOPOLOGICAL DATA*********************)

topdata := Block[{h,vol,red,amb,nor,hip,dd,ww,ll},
   conv;
   kopl=Map[#1-1&,Flatten[Array[a,Table[4,{i,1,$m}]]]/.a->List];
   deg2=Map[h^3 xlinamo[#1]&,Select[kopl,(Plus @@ #1)==2&]]/.h->1;
   kopl=Map[h^3 xlinamo[#1]&,Select[kopl,(Plus @@ #1)==3&]];
   vol=Product[Sum[dd[[i,j]] h J[i],{i,1,$m}],
         {j,1,Length[dd[[1]]]}];
   red=Product[Product[ww[[i,j]],{j,1,Length[ww[[i]]]}],
         {i,1,$m}];
   amb=Product[Product[(1+ww[[i,j]] h J[i]),{j,1,Length[ww[[i]]]}],
         {i,1,$m}];
   nor=Product[(1+Sum[dd[[i,j]] h J[i],{i,1,$m}]),
                  {j,1,Length[dd[[1]]]}];
   hip=xlinamo[Table[Length[ww[[i]]]-1,{i,1,Length[ww]}]];
   ll=Length[Flatten[ww]]-$m;
   c3=Coefficient[Expand[Normal[Series[amb/nor,{h,0,3}]]],h^3];
   eul=Coefficient[Expand[(Normal[Series[ h^3*c3*
                        vol/red,{h,0,ll}]]/.{ h->1})],hip];
   Print["Topological Data :  "];
   Print["Eulernumber == ",eul];
   Print["Evaluation of second Chern class : "];
   c2=Coefficient[Expand[Normal[Series[amb/nor,{h,0,2}]]],h^2];
   Do[
      c2j[i]=Coefficient[Expand[(Normal[Series[ h^3 c2 J[i]*
                        vol/red,{h,0,ll}]]/.{ h->1})],hip];
      Print["c2.J[",i,"]     ==  ",c2j[i]],
   {i,1,$m}];
   Print["Topological couplings: "];
   Do[
      kk[i]=Coefficient[Expand[(Normal[Series[kopl[[i]]*
                        vol/red,{h,0,ll}]]/.{h->1})],hip];
      Print[(kopl[[i]])/.h->1,"      == ",kk[i]],
   {i,1,Length[kopl]}];
   kopl=kopl/.h->1;
   ring=Sum[kk[i]*kopl[[i]],{i,1,Length[kopl]}];
   topring=Sum[kk[i]*kopl[[i]]/Apply[Times,Map[#!&,xmonali[kopl[[i]]]]],
           {i,1,Length[kopl]}];
   Do[pv[i]=D[topring,J[i]],{i,1,$m}]]



conv := Block[{i},
        dd=Table[-$l[[i,1]],{i,1,$m}];
        ww=Table[Select[$l[[i,2]],#=!=0&],{i,1,$m}]]


xmonali/: xmonali[p_]:=Table[Exponent[p,J[ix]],{ix,1,$m}]
xlinamo/: xlinamo[l_List]:= Inner[Power,$Js,l,Times];


(************CALCULATION OF THE YUKAWA COUPLING EXPANSION*******)

yukawa := Block[{hi,logs,ws,cc,df},
        Do[
          S[i]=Series[S[i],{h,0,$no-1}]
        ,{i,1,$m}];
        S[0]=Series[S[0],{h,0,$no-1}];
        Do[
           S[i,j]=S[i,j]-sexpand[S[i]*S[j]/S[0]];
           S[i,j]=Normal[sexpand[S[i,j]/S[0]]]/.h->1,
        {j,1,$m},{i,1,j}];
        Print["Evaluating the inverse series"];
        invert;
        Print["Evaluating the threepoint functions"];
        Do[
          df[i]=Sum[Coefficient[pv[i],deg2[[j]]]*
          S[Sequence @@ convlist[xmonali[ deg2[[j]] ] ]],{j,1,Length[deg2]}];
          df[i]=Expand[nsexpand[Normal[df[i]/.rzq]]],
        {i,1,$m}];
       Do[
         nyk[i,j,k]=Expand[dq[dq[df[k],i],j]],
        {k,1,$m},{j,1,k},{i,1,j}];alcc]

dq[ex_,i_]:= q[i]*D[ex,q[i]]

convlist[l_]:=Block[{deg},
              deg=Plus @@ l;
              cl={};
              Do[
                hi=Flatten[Map[Table[#1,{k,1,i}]&,Flatten[Position[l,i]]]];
                cl=Flatten[Join[cl,hi]],
              {i,1,deg}];Return[Sort[cl]]]
(*****************************************************************)



(**************************INVERSION OF SERIES *******************)



invert := Block[{st},
     If[$m == 1,
        rzq= {z[1] ->
              InverseSeries[Series[(qs[1]/.{h->1,z[1]->x}),{x,0,$no}],q]}/.
              q->q[1],
       getfpow[$no];
       Do[
         zn[i]=Expand[q[i]*h - qs[i] + h*z[i]],
      {i,1,$m}];
       Do[
          Print[st];
          Do[
            Do[
              zn[i]=sub[zn[i],k]
              (*Print[Length[zn[i]]]*),
            {k,1,$m}],
          {i,1,$m}],
       {st,2,$no}];
       rzq=Table[z[i]->sexpand[Series[zn[i],{h,0,$no}]],{i,1,$m}]]]


sub[ex_,k_] := Block[{sx,er,ccs},
                sx=zn[k];
                (*Print["k ",k];
                Print[ex];*)
                Do[cc[i]=Coefficient[sx,h^i],{i,1,$no}];
                res=ex/.z[k]->0;
                er=ex-res;
                Do[
                  (*Print[h^l];*)
                  ccs=Coefficient[er,h^l];
                  (*Print["ccs  ",ccs];*)
                  ccs=ccs/.{z[k]^e_:>pow[e,$no-l+e]/h^e,z[k]->
                                  pow[1,$no-l+1]/h};
                  res=res+ccs*h^l,
                {l,2,$no}];
               Return[Expand[res]]]

pow[n_,ord_] := Block[{hi},
                   (*Print["pow: ",n,"  ",ord];*)
                   If[n>ord,
                     pres=0,
                     pres=fp[n,ord]/.Table[cs[i]->cc[i],{i,1,ord-n+1}]];
                   (*Print["pres ",pres];*)
                   Return[pres]]

getfpow[n_]:=Block[{sh,cs},
           sh=Sum[cs[i] h^i,{i,1,n}];
           Do[
             fph=sexpand[Series[sh^l,{h,0,n}]];
             Do[
               fp[l,k]=Normal[Series[fph,{h,0,k}]],
             {k,1,n}]
           ,{l,1,n}]]

sexpand[ss_] := Block[{s},
                    s=ss;
                    Do[s[[3,i]]=Expand[s[[3,i]]],{i,1,Length[s[[3]]]}];
                    Return[s]]
nsexpand[ss_] :=Block[{s},
                    s=ss;
                    If[Head[s] === SeriesData,Do[s[[3,i]]=Expand[s[[3,i]]],
                    {i,1,Length[s[[3]]]}],];
                    Return[s]]

expa/: expa[st_] := Series[(st/.{z[i_] :> h z[i]}),{h,0,$no}];



(*********************EXTRACTION OF THE INSTANTON NUMBERS**************)



Xmonali/: Xmonali[p_]:=Table[Exponent[p,q[ix]],{ix,1,$m}]
Xlinamo/: Xlinamo[l_List]:= Inner[Power,$qs,l,Times];




init/: init := Block[{hi},
               Do[
               nkk[k1,k2,k3]= Expand[Normal[nyk[k1,k2,k3]]/.{h->1}],
               {k1,1,$m},{k2,k1,$m},{k3,k2,$m}]]





irc[il_] := Block[{hi},
             gcd=GCD[Sequence @@ il];
             ilr=il/gcd;
             lil=Length[il];
             mon=Xlinamo[il];
             Do[If[il[[i]] =!= 0, lyk={i,i,i}; npos=i, ],{i,lil,1,-1}];
                  pf=il[[npos]]^3/(gcd)^3;
                  nuc=seln[Coefficient[nkk[Sequence @@ lyk],mon]]/pf;
                  (*Print[il,"  ",nuc,"  ",lyk];*)
                  If[gcd =!= 1,
                  Do[
                   If[Mod[gcd,pp] == 0,
                   (*Print[pp,"  ", qsnc[ilr*pp]];*)
                   nuc=nuc-pp^3*qsnc[ilr*pp],],
                 {pp,1,gcd-1}]];
                 Return[nuc=nuc/(gcd)^3]]

qsnc/: qsnc[sil_] := Block[{sgcd,silr,srt,smon,spf,snc},
             sgcd=GCD[Sequence @@ sil];
             silr=sil/sgcd;
             smon=Xlinamo[sil];
             (*Print["smon  ",smon];*)
             spf=sil[[npos]]^3/(sgcd)^3;
             sncr=seln[Coefficient[nkk[Sequence @@ lyk],smon]]/spf;
             If[sgcd =!= 1,
               Do[
                   If[Mod[sgcd,spp] == 0,sncr=sncr-spp^3*qsnc[silr*spp],],
                 {spp,1,sgcd-1}]];
                 Return[sncr=sncr/(sgcd)^3]]

seln[ps_]:= Select[ps+dh1+dh2,NumberQ[#1]&]

snc/: snc[sil_] := Block[{sgcd,silr,srt,smon,spf,snc},
             sgcd=GCD[Sequence @@ sil];
             silr=sil/sgcd;
             smon=Product[q[i]^sil[[i]],{i,1,lil}];
             (*Print["smon  ",smon];*)
             spf=Product[If[sil[[i]] =!= 0, sil[[i]]^np[i],1],{i,1,lil}]/
                 (sgcd)^3; (*Print["spf  ",spf];*)
             sncr=seln[Coefficient[nkk[Sequence @@ lyk],smon]]/spf;
             (*Print["sncr ",sncr];*)
             If[sgcd =!= 1,
               Do[
                   If[Mod[sgcd,spp] == 0,sncr=sncr-spp^3*snc[silr*spp],],
                 {spp,1,sgcd-1}]];
                 Return[sncr=sncr/(sgcd)^3]]


allc[i1_,i2_,i3_] := Block[{hi},
                     ref=Table[0,{i,1,$m}];
                     lnkk=Length[nkk[i1,i2,i3]];
                     deg=Table[Xmonali[nkk[i1,i2,i3][[i]]],{i,1,lnkk}];
                     Do[If[deg[[i]] =!= ref,
                          nofc=irc[deg[[i]]];
                          Print[deg[[i]],"  ",nofc];
                          $samcu=Join[{{deg[[i]],nofc}},$samcu],],
                      {i,1,lnkk}]]


alcc := Block[{hi},
        Print["Extracting the instantons"];
        init;
        $samcu={};
        Do[Print["  "];Print["From Y",i1,i2,i3," we obtain"];
            allc[i1,i2,i3],
        {i1,1,$m},{i2,i1,$m},{i3,i2,$m}];
        $samcu=Union[$samcu];
        $samcu=Sort[$samcu,! OrderedQ[{Abs[#1[[1]]],Abs[#2[[1]]]}]];
        $samcu=
        Sort[$samcu,(GCD[Sequence @@ #1[[1]]] <=
             GCD[Sequence @@ #2[[1]]])&];
        Do[nrc[Sequence @@ $samcu[[i,1]]]=$samcu[[i,2]],
        {i,1,Length[$samcu]}]]

(**********************************************************************)




(*********************EXTRACTION OF THE INSTANTON NUMBERS**************)


Xmonali/: Xmonali[p_]:=If[p=!=0,Table[Exponent[p,q[ix]],{ix,1,$m}],0]
Xlinamo/: Xlinamo[l_List]:= Inner[Power,$qs,l,Times];




init/: init := Block[{hi},
               Do[
               nkk[k1,k2,k3]= Expand[Normal[nyk[k1,k2,k3]]/.{h->1}],
               {k1,1,$m},{k2,k1,$m},{k3,k2,$m}]]





irc[il_] := Block[{hi},
             gcd=GCD[Sequence @@ il];
             ilr=il/gcd;
             lil=Length[il];
             mon=Xlinamo[il];
             Do[If[il[[i]] =!= 0, lyk={i,i,i}; npos=i, ],{i,lil,1,-1}];
                  pf=il[[npos]]^3/(gcd)^3;
                  nuc=seln[Coefficient[nkk[Sequence @@ lyk],mon]]/pf;
                  (*Print[il,"  ",nuc,"  ",lyk];*)
                  If[gcd =!= 1,
                  Do[
                   If[Mod[gcd,pp] == 0,
                   (*Print[pp,"  ", qsnc[ilr*pp]];*)
                   nuc=nuc-pp^3*qsnc[ilr*pp],],
                 {pp,1,gcd-1}]];
                 Return[nuc=nuc/(gcd)^3]]

qsnc/: qsnc[sil_] := Block[{sgcd,silr,srt,smon,spf,snc},
             sgcd=GCD[Sequence @@ sil];
             silr=sil/sgcd;
             smon=Xlinamo[sil];
             (*Print["smon  ",smon];*)
             spf=sil[[npos]]^3/(sgcd)^3;
             sncr=seln[Coefficient[nkk[Sequence @@ lyk],smon]]/spf;
             If[sgcd =!= 1,
               Do[
                   If[Mod[sgcd,spp] == 0,sncr=sncr-spp^3*qsnc[silr*spp],],
                 {spp,1,sgcd-1}]];
                 Return[sncr=sncr/(sgcd)^3]]

seln[ps_]:= Select[ps+dh1+dh2,NumberQ[#1]&]

snc/: snc[sil_] := Block[{sgcd,silr,srt,smon,spf,snc},
             sgcd=GCD[Sequence @@ sil];
             silr=sil/sgcd;
             smon=Product[q[i]^sil[[i]],{i,1,lil}];
             (*Print["smon  ",smon];*)
             spf=Product[If[sil[[i]] =!= 0, sil[[i]]^np[i],1],{i,1,lil}]/
                 (sgcd)^3; (*Print["spf  ",spf];*)
             sncr=seln[Coefficient[nkk[Sequence @@ lyk],smon]]/spf;
             (*Print["sncr ",sncr];*)
             If[sgcd =!= 1,
               Do[
                   If[Mod[sgcd,spp] == 0,sncr=sncr-spp^3*snc[silr*spp],],
                 {spp,1,sgcd-1}]];
                 Return[sncr=sncr/(sgcd)^3]]


allc[i1_,i2_,i3_] := Block[{hi},
                     ref=Table[0,{i,1,$m}];
                     lnkk=If[FreeQ[FullForm[nkk[i1,i2,i3]],Plus] &&
                             nkk[i1,i2,i3]=!=0, 1,
                             Length[nkk[i1,i2,i3]] ];
                     deg=If[lnkk==1,{Xmonali[nkk[i1,i2,i3]]},
                         Table[Xmonali[nkk[i1,i2,i3][[i]]],{i,1,lnkk}]];
                     Do[If[deg[[i]] =!= ref,
                          nofc=irc[deg[[i]]];
                          Print[deg[[i]],"  ",nofc];
                          $samcu=Join[{{deg[[i]],nofc}},$samcu],],
                      {i,1,lnkk}]]

kttt[i1_,i2_,i3_]:= Module[{prod},
                     prod[s1_,s2_,s3_]:=s1*s2*s3;
                     ref=Table[0,{i,1,$m}];
                     lnkk=If[FreeQ[FullForm[nkk[i1,i2,i3]],Plus], 1,
                             Length[nkk[i1,i2,i3]] ];
                     deg=If[lnkk==1,{Xmonali[nkk[i1,i2,i3]]},
                         Table[Xmonali[nkk[i1,i2,i3][[i]]],{i,1,lnkk}]];
               If[lnkk=!=1 && deg[[1]]=!=ref && deg[[1]]=!=0,
                (nkk[i1,i2,i3][[1]]/.{q[aa_]->0}) +
                  Sum[
                   irc[deg[[i]]]*
                    prod[deg[[i]][[i1]],deg[[i]][[i2]],deg[[i]][[i3]] ]*
                      Xlinamo[deg[[i]]]/(1-Xlinamo[deg[[i]]]),
                       {i,1,lnkk}],
               If[lnkk==1 && deg[[1]]=!=ref && deg[[1]]=!=0,
                   irc[deg[[1]]]*
                    prod[deg[[1]][[i1]],deg[[1]][[i2]],deg[[1]][[i3]] ]*
                     Xlinamo[deg[[1]]]/(1-Xlinamo[deg[[1]]]),
               If[deg[[1]]=!=0,
                  (nkk[i1,i2,i3]/.{q[a_]->0}) +
                   Sum[
                    irc[deg[[i]]]*
                     prod[deg[[i]][[i1]],deg[[i]][[i2]],deg[[i]][[i3]] ]*
                       Xlinamo[deg[[i]]]/(1-Xlinamo[deg[[i]]]),
                       {i,2,lnkk}],
                  nkk[i1,i2,i3]
                    ]]]        ]

ktttexp[i1_,i2_,i3_]:=ReleaseHold[Normal[Series[
                      kttt[i1,i2,i3]/.{q[s_]->h q[s]},{h,0,$no}]]/.h->1]

(******************************************************)

alcc := Block[{hi},
        Print["Extracting the instantons"];
        init;
        $samcu={};
        Do[Print["  "];Print["From Y",i1,i2,i3," we obtain"];
            allc[i1,i2,i3],
        {i1,1,$m},{i2,i1,$m},{i3,i2,$m}];
        $samcu=Union[$samcu];
        $samcu=Sort[$samcu,! OrderedQ[{Abs[#1[[1]]],Abs[#2[[1]]]}]];
        $samcu=
        Sort[$samcu,(GCD[Sequence @@ #1[[1]]] <=
             GCD[Sequence @@ #2[[1]]])&];
        Do[nrc[Sequence @@ $samcu[[i,1]]]=$samcu[[i,2]],
        {i,1,Length[$samcu]}]]

(**********************************************************************)


(*******************OTHER UTILITIES******************************)

 prepot := Block[{hi},
                prep={};
                Do[
                   part=Expand[Integrate[Expand[Integrate[
                        Expand[Integrate[nyk[i,j,k]/q[i],q[i]]]/
                        q[j],q[j]]]/q[k],q[k]]];
                Print[part];
                If[Head[part]===Plus,part=part/.Plus->List,part={part}];
                prep=Union[prep,part/.Plus->List],
                {i,1,$m},{j,i,$m},{k,j,$m}];
                prepotential=prep/.List->Plus;
                prepotential=prepotential+(topring/.J[a__]:>Log[q[a]])]






findring[ort_] := Block[{degrees,degi,orth},
       orth=ort/.{o[i_]:>J[i]};
       degrees=Map[#1-1&,Flatten[Array[a,Table[4,{i,1,$m}]]]/.a->List];
       Do[
         freering[i]=Map[xlinamo[#1]&,Select[degrees,(Plus @@ #1)==i&]],
       {i,0,3}];
       orth3={};
       Do[
         degi=getdeg[orth[[i]]];
         fr=freering[3-degi];
         Do[
           orth3=Join[orth3,{convpolyvec[fr[[j]]*orth[[i]],3]}],
         {j,1,Length[fr]}],
       {i,1,Length[orth]}];
      ring=convvecpoly[Flatten[NullSpace[orth3]],3];
      cjjj=Coefficient[Expand[ring],J[1]^3];
      Print["Normalisation of the J^3 intersection required"];
      norm=Input[];
      ring=Expand[(norm/cjjj)*ring];
      tr[ring]]

getdeg[poly_]:=Block[{hi},
               npoly=Expand[poly];
               If[Head[npoly]===Plus,
                 deg=Max[Map[Apply[Plus,xmonali[#1]]&,Apply[List,npoly]]],
                 deg=Apply[Plus,xmonali[npoly]]];
               Return[deg]]


convpolyvec[poly_,deg_] := Block[{hi},
           dimv=Length[freering[deg]];
           npoly=Expand[poly];
           vec=Table[Coefficient[npoly,freering[deg][[i]]],{i,1,dimv}]]

convvecpoly[vec_,deg_] := Block[{hi},
           Inner[Times,vec,freering[deg],Plus]]



tr[ring_]:= Block[{hi,hlist},
          hlist=Map[#1-1&,Flatten[Array[a,Table[4,{i,1,$m}]]]/.a->List];
          deg2=Map[xlinamo[#1]&,Select[hlist,(Plus @@ #1)==2&]];
          kopl=Map[xlinamo[#1]&,Select[hlist,(Plus @@ #1)==3&]];
          topring=Sum[ring[[i]]/Apply[Times,Map[#!&,xmonali[ring[[i]]]]],
                   {i,1,Length[ring]}];
          Do[pv[i]=D[topring,J[i]],{i,1,$m}]]



checktop := Block[{hli},
            li0=Map[Drop[#1,-1][[1]]&,$l];
            li1=Map[Drop[#1,1][[1]]&,$l];
            Do[
              clk[Sequence @@ convlist[xmonali[kopl[[j]]]]]=
              Coefficient[ring,kopl[[j]]],
            {j,1,Length[kopl]}];
            ceul=(1/3)*Sum[clk[i1,i2,i3]*(
                          Sum[li0[[i1,ii]]*li0[[i2,ii]]*li0[[i3,ii]],
                             {ii,1,Length[li0[[i1]]]}]+
                          Sum[li1[[i1,ii]]*li1[[i2,ii]]*li1[[i3,ii]],
                                {ii,1,Length[li1[[i1]]]}]),
                         {i1,1,$m},{i2,1,$m},{i3,1,$m}];
            Print["Eulernumber = ",ceul];
            Do[
                c2J[i1]=(1/2)*Sum[clk[i1,i2,i3]*(
                          Sum[li0[[i2,ii]]*li0[[i3,ii]],
                             {ii,1,Length[li0[[i2]]]}]-
                          Sum[li1[[i2,ii]]*li1[[i3,ii]],
                                {ii,1,Length[li1[[i1]]]}]),
                 {i2,1,$m},{i3,1,$m}];
             Print["c2j[",i1,"] = ",c2J[i1]],
             {i1,1,$m}]]







clk[i_, j_, k_] := clk[Apply[Sequence, Sort[{j, i, k}]]]