(***********************************************************************)
(*                                                                     *)
(*                PROGRAM:        INSTANTON    (1.?)                   *)
(*                                                                     *)
(*                   send suggestions to A. Klemm                      *)
(*                     aklemm@physics.wisc.edu                         *)
(*                                                                     *)
(***********************************************************************)



(*   Some improvements to INSTANTON: it now checks integrability:        

     This is a version, which does not requieres any more the ring of
     topological triple couplings nor the principal part of the differential
     operators. 

     If these data are not given it uses the fact that the Instanton 
     corrected couplings  must integrate to a prepotential to fix the 
     relative topological constants. The overall constant may be fixed 
     using the Eulernumber.       

     The usage is :
       
     cases[No. of the model, Expansion precision],

     where for the model with No. (No. of model) the moricone 
     has to be predefined as in the examples below. 
     
     If the program does find in addition ring[No. of model] or 
     orth[No. of model] it uses these data (ring is prefered if both are
     given) and checks integrability.

     If these data are not found, the ring is generated from the 
     integrability constraints.   
        
    
**************************************************************************)



(*  	   with "newinvert" by M. Kreuzer / Tue Nov 26 2002                 *)

ITime:=FIRSTTIME=LASTTIME=AbsoluteTime[];	(* Initialize realtime info *) 
RTime:=Block[{rtime=Round[AbsoluteTime[]-LASTTIME]},(* incremental Realtime *)
             LASTTIME=AbsoluteTime[];Return[rtime]];
TTime:=Block[{rtime=Round[AbsoluteTime[]-FIRSTTIME]},	  (* Total realtime *)
             FIRSTTIME=AbsoluteTime[];Return[rtime]];
MBmem:=Quotient[MemoryInUse[],1000000];			    (* Memory in MB *)


(*                 DATA  FOR  THE EXAMPLES                             *)

(* LG examples: Notation lm[weights], ring[weights] is moricone and
   intersection ring for the examples *) 


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

ring[11169]=9*J[1]^3 + 3*J[1]^2*J[2] + J[1]*J[2]^2 

lm[quintic]={{{-5},{1,1,1,1,1}}}

ring[quintic]=5*J[1]^3  

lm[sextic]={{{-6},{1,1,1,1,2}}}

ring[sextic]=3*J[1]^3  




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


(**************************INPUT ROUTINES******************************)
Clear[nrc]



cases[n_,in2_,$str_,lflag_] := Block[{hi},
   If[Head[lflag]==0,$lflag=0;cases[n,in2,$str],$lflag=1;lcases[n,in2,$str]]];
  

cases[n_,in2_,$str_] := Block[{hi},       			ITime;
       in1=2;
       $lflag=0;
       $strf:=$str;
       $strflong:=StringJoin["Long",$str]; 
       Clear[nrc]; 
       $l=lm[n];
       If[Head[lb[n]]===lb,
       $limit={Table[lx0[i],{i,1,Length[$l[[1,1]] ]}],
           Table[lx1[i],{i,1,Length[$l[[1,2]]]}]},   
           $limit=lb[n]];
       Print["Limit ",$limit];
       $orth=orth[n];
       $m=Length[$l];
       $Js=Table[J[i],{i,1,$m}];
       $qs=Table[q[i],{i,1,$m}];
       Print["findring=0, series=1 and full=2 "];
       flag=in1;  
       If[flag==0 || flag==2,
       If[Head[ring[n]]===ring,
          If[Head[orth[n]]===orth,
            Print["No data for the ring or the Ideal of 
                   differential operators"];
            $ring=setring[$m];             
            $findringflag=1;
            tr[$ring,3] 
          ,
          findring[$orth,3];
          $findringflag=0]
        ,
        $ring=ring[n];
        tr[$ring,3];
        $findringflag=0
        ],00];
       If[flag==1 || flag==2,
       Print["Expansion Precision, the highest (multi)degree evaluated"];
       $no=in2;
       Print["Evaluating the Period ... after ",RTime," sec"];
       If[flag==1,pals,pals],00];
       If[flag==2,
       yukawa;
       eqlout=Union[$eql]; 
       Write[$strf,"(*********************************************)"];
       Save[$strf,eqlout];
       Write[$strf,"(*********************************************)"];
       Save[$strf,$treq];
       Write[$strf,"(*********************************************)"];
       Save[$strf,$l];
       Write[$strf,"(*********************************************)"];
       Save[$strf,$ring];
       Write[$strf,"(*********************************************)"];
       Save[$strf,nrc];
       Write[$strflong,"(*********************************************)"];
       Save[$strflong,rzq];
       Write[$strflong,"(*********************************************)"];
       Save[$strflong,qs];
       Write[$strflong,"(*********************************************)"];
       Save[$strflong,Sr]; 
       Write[$strf,"(*******************************************)"];
       Save[$strflong,df];
       Write[$strflong,"(*********************************************)"];
       Save[$strflong,nyk],00];	Print["Done                            ", 
		RTime," sec for 3pt fun -> total real time: ",TTime," sec"]];

(* this is the local version *)
lcases[n_,in2_,$str_] := Block[{hi},    
       Clear[nyk];     
       in1=2;
       $lflag=1;
       $strf:=StringJoin["local",$str];
       $strflong:=StringJoin["locallong",$str]; 
       Clear[nrc]; 
       $l=lm[n];
       $lc=lc[n]; 
       If[Head[lb[n]]===lb,
       $limit={Table[lx0[i],{i,1,Length[$l[[1,1]] ]}],
           Table[lx1[i],{i,1,Length[$l[[1,2]]]}]},   
           $limit=lb[n]];
       $orth=orth[n];
       $m=Length[$l];
       $Js=Table[J[i],{i,1,$m}];
       $qs=Table[q[i],{i,1,$m}];

       Print["findring=0, series=1 and full=2 "];
       flag=in1;  
       If[flag==0 || flag==2,
       If[Head[ring[n]]===ring,
          If[Head[orth[n]]===orth,
            Print["No data for the ring or the Ideal of 
                   differential operators"];
            $ring=setring[$m,2];             
            $findringflag=1;
            tr[$ring,3-1] 
          ,
          findring[$orth,3-1];
          $findringflag=0]
        ,
        $mc=1; 
        $ring=ring[n]; 
        tr[$ring,3-1];
        $findringflag=0
        ],00];
       If[flag==1 || flag==2,
       Print["Expansion Precision, the highest (multi)degree evaluated"];
       $no=in2;
       Print["Evaluating the Period"];
       If[flag==1,pals,pals],00];
       If[flag==2,
       lyukawa;
       eqlout=Union[$eql]; 
       Write[$strf,"(*********************************************)"];
       Save[$strf,eqlout];
       Write[$strf,"(*********************************************)"];
       Save[$strf,$treq];
       Write[$strf,"(*********************************************)"];
       Save[$strf,$l];
       Write[$strf,"(*********************************************)"];
       Save[$strf,$ring];
       Write[$strf,"(*********************************************)"];
       Save[$strf,nrc];
       Write[$strf,"(*********************************************)"];
       Save[$strflong,rzq];
       Write[$strf,"(*********************************************)"];
       Save[$strflong,df];
       Write[$strflong,"(*********************************************)"];
       Save[$strflong,nyk],00]]



setring[m_]:=Block[{hi},
$Js=Table[J[i],{i,1,m}];
lpw=Map[#1-1&,Flatten[Array[a,Table[4,{i,1,m}]]]/.{a->List}];
lpw=Select[lpw,Plus @@ #1 == 3&];
$unbek=Table[aa[i],{i,1,Length[lpw]}];
$ring=Sum[aa[i]*xlinamo[lpw[[i]]],{i,1,Length[lpw]}]]

setring[m_,d_]:=Block[{hi},
$Js=Table[J[i],{i,1,m}];
lpw=Map[#1-1&,Flatten[Array[a,Table[d+1,{i,1,m}]]]/.{a->List}];
lpw=Select[lpw,Plus @@ #1 == d&];
$unbek=Table[aa[i],{i,1,Length[lpw]}];
$ring=Sum[aa[i]*xlinamo[lpw[[i]]],{i,1,Length[lpw]}]]


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

 

pals/: pals := Block[{hi,mon,li0,li1,np,lh,pg0,pg1,dpg1,dpg0,
                     mk1,mk2,ip,isum},
   Clear[Sr];
   lim0=-$limit[[1]];
   lim1= $limit[[2]];
   (*Print["lim0 ",lim0," lim1 ",lim1];*)
   li0=-Map[Drop[#1,-1][[1]]&,$l]; (* {{-l_0^{(1)}},...{-l_0^{(k)}} *)
   li1= Map[Drop[#1,1][[1]]&,$l];
   (*Print["li0 ",li0," li1 ",li1];*)
   var=Table[z[i],{i,1,$m}];
   Do[Sr[i]=0;Do[Sr[i,j]=0,{j,1,$m}],{i,0,$m}];
   Do[ip[i]=0,{i,1,$m}];
   isum:=-1;
   ip[1]=-1;
   is1=1;
   While[is1<=$m,
     If[isum<$no,
     ip[is1]=ip[is1]+1;
     lh=Table[ip[j],{j,1,$m}]; 
     mon=Inner[Power,var,lh,Times]/.z[u_]:>h*z[u];
     (*Print[mon];*)
     nl0=Inner[Times,lh,li0];
     nl1=Inner[Times,lh,li1];
     np0=Flatten[Position[Negative[nl0],True]]; 
     np1=Flatten[Position[Negative[nl1],True]];
     nn0=Select[nl0,Negative[#1]&];
     pn0=Select[nl0,#>=0&];
     nn1=Select[nl1,Negative[#1]&];
     pn1=Select[nl1,#>=0&];  
     (*Print["nl0/nl1 ",nl0," / ",nl1," pn0/pn1 ",pn0," / ",pn1];*)
     (*Print["nn0/nn1 ",nn0," / ",nn1," np0/np1 ",np0," / ",np1];*)
     case=Length[nn1]-Length[nn0];
     If[case < 0,Print["Inconsistent mori generators"],00  ];
       facl=(-1)^(Apply[Plus,nn1]+Apply[Plus,nn0])*
            Apply[Times,Map[lim1[[#]]&,np1]]/
            Apply[Times,Map[lim0[[#]]&,np0]];
       w0=facl*
          Apply[Times,Map[#!&,pn0]]/
          Apply[Times,Map[#!&,pn1]]*
          Apply[Times,Map[(-1-#)!&,nn1]]/
          Apply[Times,Map[(-1-#)!&,nn0]];
     If[case == 0,
          Sr[0]=Sr[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];
            Sr[b]=Sr[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];
            Sr[c,b]=Sr[c,b] + w0*(Ac*Ab+Bbc)*mon,
         {c,1,b}],
        {b,1,$m}]
       ,
     (* Excess of zero's is one *)
     If[case == 1,
           If[Length[nn0]>0,Print["consistency check required"], 00];
           (*Print["case 1"];*)
           pg0=Map[pg[#]&,pn0];
           pg1=Map[pg[#]&,pn1];
           pgmk1=Map[pg[-nl1[[#]]-1]&,np1];
           pgmk0=Map[pg[-nl0[[#]]-1]&,np0];
           Do[
             Sr[b]=Sr[b]+w0*mon*
                   (Sum[(*  (-1)^(-nl1[[np1[[i1]]]]+1)*  *)
                      -li1[[b]][[np1[[i1]]]]/
                      lim1[[np1[[i1]]]],
                   {i1,1,Length[np1]}]+
                    Sum[(*   (-1)^(-nl0[[np0[[i1]]]])*  *)
                        li0[[b]][[np0[[i1]]]]/
                       lim0[[np0[[i1]]]],
                   {i1,1,Length[np0]}]);                      
                 Ab=Inner[Times,Drop[li0[[b]],np0],pg0]+
                 Inner[Times,Map[li0[[b]][[#]]&,np0],pgmk0]-
                 Inner[Times,Drop[li1[[b]],np1],pg1]-
                 Inner[Times,Map[li1[[b]][[#]]&,np1],pgmk1];
              Do[
              Ac=Inner[Times,Drop[li0[[c]],np0],pg0]+
                 Inner[Times,Map[li0[[c]][[#]]&,np0],pgmk0]-
                 Inner[Times,Drop[li1[[c]],np1],pg1]-
                 Inner[Times,Map[li1[[c]][[#]]&,np1],pgmk1];
                 (*Print["Ab/Ac", Ab,"  / ", Ac];*)
                 Sr[c,b]=Sr[c,b]+w0*mon*
                    (Ac*(
                    Sum[
                  (*(-1)^(-nl1[[np1[[i1]]]]+1)*)
                     -li1[[b]][[np1[[i1]]]]/
                      lim1[[np1[[i1]]]],
                    {i1,1,Length[np1]}]+
                    Sum[
                     (*(-1)^(-nl0[[np0[[i1]]]])*) 
                     li0[[b]][[np0[[i1]]]]/
                     lim0[[np0[[i1]]]],
                    {i1,1,Length[np0]}])+
                    Ab*(
                    Sum[
                     (*(-1)^(-nl1[[np1[[i1]]]]+1)*)
                     -li1[[c]][[np1[[i1]]]]/
                     lim1[[np1[[i1]]]],
                    {i1,1,Length[np1]}]+
                    Sum[
                     (*(-1)^(-nl0[[np0[[i1]]]])*)
                     li0[[c]][[np0[[i1]]]]/
                     lim0[[np0[[i1]]]],
                    {i1,1,Length[np0]}])),
           {c,1,b}],
        {b,1,$m}]
      ,
       If[case == 2,
        (*Print["case 2"];*)
        araw=Apply[Times,Map[#!&,pn0]]/
             Apply[Times,Map[#!&,pn1]]*
             Apply[Times,Map[(-1-#)!&,nn1]]/
             Apply[Times,Map[(-1-#)!&,nn0]];
        Do[
           sum11=Sum[(*(-1)^(nl1[[np1[[i1]]]]+nl1[[np1[[i2]]]])*)
                  (li1[[c]][[np1[[i1]]]]*li1[[b]][[np1[[i2]]]]+
                   li1[[b]][[np1[[i1]]]]*li1[[c]][[np1[[i2]]]])/
                   (lim1[[np1[[i1]]]]*lim1[[np1[[i2]]]]),
                   {i2,1,Length[np1]},{i1,1,i2-1}];
           sum00=Sum[(*(-1)^(nl0[[np0[[i1]]]]+nl0[[np0[[i2]]]])*)
                  (li0[[c]][[np0[[i1]]]]*li0[[b]][[np0[[i2]]]]+
                   li0[[b]][[np0[[i1]]]]*li0[[c]][[np0[[i2]]]])/
                   (lim0[[np0[[i1]]]]*lim0[[np0[[i2]]]]),
                   {i2,1,Length[np0]},{i1,1,i2}];
           sum01=Sum[(*(-1)^(nl0[[np0[[i1]]]]+nl1[[np1[[i2]]]])**)
                  -(li0[[c]][[np0[[i1]]]]*li1[[b]][[np1[[i2]]]]+
                   li0[[b]][[np0[[i1]]]]*li1[[c]][[np1[[i2]]]])/
                   (lim0[[np0[[i1]]]]*lim1[[np1[[i2]]]]),
                   {i1,1,Length[np0]},{i2,1,Length[np1]}];
     (*Print[c,"/",b," sums ",sum11,"  ",sum01,"  ",sum00," araw ",araw];*)
           Sr[c,b]=Sr[c,b]+w0*(sum11+sum00+sum01)*mon,
         {b,1,$m},{c,1,b}]
     ,00],00],00];
     isum=isum+1;
     is1=1
    ,
     isum=isum-ip[is1];
     ip[is1]=0;
     is1=is1+1;
     ]];
    Write[$strflong,"(*********************************************)"];
    Save[$strflong,Sr]; 
    Print["evaluate qs                     ", RTime," sec for Periods"];
    Do[
      qs[i]=Normal[sexpand[Series[h z[i] Exp[Sr[i]/Sr[0]],{h,0,$no}]]],
    {i,1,$m}];
    Save[$strflong,qs]]

flipsign:=Block[{hi},Do[qs[i]=qs[i]/.z[1]->-z[1];Sr[i]=Sr[i]/.z[1]->-z[1];Do[Sr[i,j]=Sr[i,j]/.z[1]->-z[1],{j,1,$m}],{i,1,$m}]]


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


pspals[m_,no_]:=Block[{hi},
      $m=m;
      $no=no;
      Do[Sr[i]=((Sr[i]/.h->1)/.z[a__]:>z[a]*h),{i,0,$m}];
      Do[
      qs[i]=Normal[sexpand[Series[h z[i] Exp[Sr[i]/Sr[0]],{h,0,$no}]]],
      {i,1,$m}];
      Do[
          Sr[i]=Series[Sr[i],{h,0,$no-1}]
        ,{i,1,$m}];
        Sr[0]=Series[Sr[0],{h,0,$no-1}];
        Do[
           Sr[i,j]=Sr[i,j]-sexpand[Sr[i]*Sr[j]/Sr[0]];
           Sr[i,j]=Normal[sexpand[Sr[i,j]/Sr[0]]]/.h->1,
        {j,1,$m},{i,1,j}];
        Print["Evaluating the inverse series"];
        invert]

   

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

degz/: degz[p_]:=Sum[Exponent[p,z[ix]],{ix,1,$m}]





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

yukawa := Block[{hi,logs,ws,cc},
        Do[
          Sr[i]=Series[Sr[i],{h,0,$no-1}]
        ,{i,1,$m}];
        Sr[0]=Series[Sr[0],{h,0,$no-1}];
        Do[
           Srr[i,j]=Sr[i,j]; 
           Sr[i,j]=Sr[i,j]-sexpand[Sr[i]*Sr[j]/Sr[0]];
           Sr[i,j]=Normal[sexpand[Sr[i,j]/Sr[0]]]/.h->1,
        {j,1,$m},{i,1,j}];
        Print["Evaluating the inverse series   ",RTime," sec for qs"];
        invert;
	Print["Invert done: ",MBmem,"MB                ",RTime,
		" sec for inversion to order O(",$no,")"]; 
        Save["inter",rzq]; 
        Print["Evaluating the threepoint functions"]; 
        Do[
          Print["calculating period entry :",i+$m+1," in q coord."];
          If[$m == 1,
          df[i]=Sum[Coefficient[pv[i],deg2[[j]]]*
          Sr[Sequence @@ convlist[xmonali[ deg2[[j]] ] ]],{j,1,Length[deg2]}];
          df[i]=Expand[nsexpand[Normal[df[i]/.rzq]]],
          df[i]=Expand[Sum[Coefficient[pv[i],deg2[[j]]]*
          Sr[Sequence @@ convlist[xmonali[deg2[[j]]]]]/.z[a__]:>h*z[a],
          {j,1,Length[deg2]}]];
          (* --------------------------
	Albrecht1: An dieser Stelle ist df[i]= h a_k z[k] 
		+ h^2 a_{jk} z[j] z[k]+... O(z^($no+1)) 
        In df[i] muss man z[k]->q[k] ersetzen. Die neue serie soll 
		bis O(q^{$no+1}) richtig sein.
        Unten eine Version, die Dein psexpandinsert benutzt, allerdings treten
	mehr z und q potenzen in den Schritten im Loop auf, als gut ist. 
	Die Instantonen kommen richtig raus. 
             ------------------------- *) 
          df[i]=Expand[df[i]/h];
df[i]=insertZ2Q[df[i]];
(*          Do[df[i]=psexpandinsert[z[k],zn[k]/h,df[i],h,$no],{k,1,$m}]; *)
          df[i]=df[i]/.h->1/.q[i__]:>q[i] h],  
        {i,1,$m}]; 
       Do[
         nyk[i,j,k]=Coefficient[Expand[$ring+dum],J[i] J[k] J[j]]+
                    Expand[dq[dq[df[k],i],j]],
        {k,1,$m},{j,1,k},{i,1,j}];       
        alcc;
        If[$findringflag==1,
           $loes=Flatten[Solve[$eql,$unbek]];
           $ring=$ring/.$loes;
           Print["Ring as determined by integrability : "];
           Print[$ring];  
           checktop;
          Print["Fix remaining constant by specifying Eulernumber or type No"];
           hieul=Input[];
           If[hieul===No, 00,
           $loes2=Flatten[Solve[hieul==ceul,$unbek]];
           $ring=$ring/.$loes2;
           Print["$ring   :",$ring];
           checktop;
           Print["You may set now ring[No. of model]=$ring"] ]   
            ,Print["Check integrability : ",Union[$eql]];
             Print["Instantons stored in in file ",$strf," and variable nrc, type: ?nrc"]]]    
         
lyukawa := Block[{hi,logs,ws,cc},
        Do[
          Sr[i]=Series[Sr[i],{h,0,$no-1}]
        ,{i,1,$m}];
        Sr[0]=Series[Sr[0],{h,0,$no-1}];
        Do[
           Srr[i,j]=Sr[i,j]; 
           Sr[i,j]=Sr[i,j]-sexpand[Sr[i]*Sr[j]/Sr[0]];
           Sr[i,j]=Normal[sexpand[Sr[i,j]/Sr[0]]]/.h->1,
        {j,1,$m},{i,1,j}];
        Print["Evaluating the inverse series"];
        invert;
        Save["inter",rzq]; 
        Print["Evaluating the threepoint functions"]; 
        Do[
          Print["calculating period entry :",i+$m+1," in q coord."];
          If[$mc == 1,
           df[i]=Sum[Coefficient[pv[i],deg2[[j]]]*
           Sr[Sequence @@ convlist[xmonali[ deg2[[j]] ] ]],{j,1,Length[deg2]}];
           df[i]=Expand[nsexpand[Normal[df[i]/.rzq]]]
          ,
           df[i]=Expand[Sum[Coefficient[pv[i],deg2[[j]]]*
           Sr[Sequence @@ convlist[xmonali[deg2[[j]]]]]/.z[a__]:>h*z[a],
          {j,1,Length[deg2]}]];
          Do[df[i]=nsub[df[i],k1],{k1,1,$m}]],  
        {i,1,$mc}]; 
       
       Do[
         nyk[i,j,k]=Coefficient[Expand[$ring+dum],J[i] J[k] J[j]]+
                    Expand[dq[dq[df[i],k],j]],
        {k,1,$m},{j,1,k},{i,1,$mc}];       
        alcc;
        If[$findringflag==1,
           $loes=Flatten[Solve[$eql,$unbek]];
           $ring=$ring/.$loes;
           Print["Ring as determined by integrability : "];
           Print[$ring];  
           checktop;
          Print["Fix remaining constant by specifying Eulernumber or type No"];
           hieul=Input[];
           If[hieul===No, 00,
           $loes2=Flatten[Solve[hieul==ceul,$unbek]];
           $ring=$ring/.$loes2;
           Print["$ring   :",$ring];
           checktop;
           Print["You may set now ring[No. of model]=$ring"] ]   
            ,Print["Check integrability : ",Union[$eql]];
             Print["Instantons stored in in file ",$strf," and variable nrc, type: ?nrc"]]]    
         

checkint:=Block[{hi},
                Do[Ri[i,j]=schwarz[i,j],{i,1,$m},{j,i+1,$m}]]

schwarz[i_,j_]:=Expand[dq[df[i],j]-dq[df[j],i]]

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 *******************)



oldinvert := Block[{st},
     If[$m == 1,
        rzq= {z[1] ->
              InverseSeries[Series[(qs[1]/.{h->1,z[1]->x}),{x,0,$no}],q]}/.
              q->q[1];
       zn[1]=Normal[z[1]/.rzq]/.q[1]->h*q[1],
       getfpow[$no];
       Print["got formal powerseries"];
       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],
            {k,1,$m}],
          {i,1,$m}],
       {st,2,$no}];
       rzq=Table[z[i]->sexpand[Series[zn[i],{h,0,$no}]],{i,1,$m}];
(*       Do[$rzq[k]=Table[z[i]->Series[rzq[[i,2]],{h,0,k}],{i,1,$m}],
         {k,1,$no}]*)
]]




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]]]


nsub[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,1,$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[{sl},
                    Print["applying series expand"]; 
                    sl=Normal[Series[ss,{h,0,$no}]];
                    Print["step1"]; 
                    sl=Series[Expand[sl],{h,0,$no}]]



sexpand[ss_] := Block[{sl},                    
                    sl=Series[ss,{h,0,$no}];
                    Do[sl[[3,i]]=Expand[sl[[3,i]]],{i,1,Length[sl[[3]]]}];
                    Return[sl]]

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

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



(************New extraction of the Instanons checks integrability************) 

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


alcc:=Block[{hi},
    init;
    Clear[nrc];
    $flagint=0;
    $eql={}; 
    $treq={}; 
    Do[
    allc[i1,i2,i3],
   {i1,1,$m},{i2,i1,$m},{i3,i2,$m}]]


init/: init := Block[{hi},
               $fser=Series[xx/(1-xx),{xx,0,$no}];
               Do[
               nyk[k1,k2,k3]= Expand[Normal[Series[nyk[k1,k2,k3],
               {h,0,$no}]]/.{h->1}],
               {k1,1,$m},{k2,1,$m},{k3,1,$m}]]





allc[i1_,i2_,i3_] := 
             Block[{hi},
             ref=Table[0,{i,1,$m}];
             saywhichcoup={i1,i2,i3};
             $nkk=nyk[i1,i2,i3]-(nyk[i1,i2,i3]/.q[i__]:>0);
             lnyk=If[FreeQ[FullForm[$nkk],Plus] && $nkk=!=0,1,Length[$nkk]];
             deg=If[lnyk==1,{Xmonali[$nkk]},
                         Table[Xmonali[$nkk[[i]]],{i,1,lnyk}]];
                     While[Length[deg] =!= 0, 
                          If[deg[[1]] =!= ref,
                           nofc=irc[{i1,i2,i3},deg[[1]],saywhichcoup];
                         ,
                           Drop[deg,1]]]]



irc[li_,il_,swc_] := Block[{hi},
             gcd=GCD[Sequence @@ il];
             ilr=il/gcd;
             lil=Length[il];
             nml=Floor[$no/Apply[Plus,ilr]];
             degm=Table[ilr*i,{i,1,nml}]; 
             degm=Select[degm,MemberQ[deg,#1]&];
             (*Print["deg    ",deg];
             Print["degm   ",degm];*)
             deg=Complement[deg,degm];
             Do[curli=degm[[i]];
               curm=Xlinamo[curli];
                (* The next two lines are a funny work-
                   arround to what I think 
                   is a mathematica bug *)
                exp=Expand[$nkk/.curm->X];
                cc=(Coefficient[exp,X])/.{q[i___]:>0};
                (*Print["li,curli,ilr,cc,curm,nml",li,curli,ilr,cc,curm,nml];*)
                If[$lflag==0,
                nrc[Sequence @@ curli]=cc/
                Product[weightf[curli[[li[[k]]]]],{k,1,3}];
                checkintegrability[curli,nrc[Sequence @@ curli],swc];
                ,
                nrc[Sequence @@ curli]=Simplify[cc/(Sum[$lc[[k]]*curli[[k]],{k,1,Length[curli]}]*
                Product[curli[[li[[k]]]],{k,2,3}])];
                (*Print[cc,"  ",swc,"  ",curli,"  ",nrc[Sequence @@ curli]]*)];                 
                (*Print[swc,"  ",curli,"  ",nrc[Sequence @@ curli]];
                Print[swc,"  ",curli,"  ",nrc[Sequence @@ curli]];*)
                $nkk=Expand[$nkk-cc*(Normal[Series[$fser,{xx,0,
                Ceiling[nml/i]}]]
                /.xx->curm)] 
             ,{i,1,Length[degm]}]]

weightf[as_]:= (If[as==0,1,as])

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


checkintegrability[cl_,value_,swc_]:=Block[{hi},
              zpos=Flatten[Position[cl,0]];                
              If[Length[Complement[Union[swc],zpos]]=!=Length[Union[swc]],
                 $eql=Join[$eql,{value==0}];
                 If[value==0,00,$treq=Join[$treq,{{cl,swc,value,0}}]],   
                 If[$flagint==0,
                  $flagint=1;
                  $curinst={value};
                  $curdeg={cl},
                  pos=Position[$curdeg,cl];
                  If[pos=={},
                    $curinst=Join[$curinst,{value}];
                    $curdeg=Join[$curdeg,{cl}],
                    $eql=Join[$eql,{value==$curinst[[pos[[1]][[1]]]]}];
                    If[value==$curinst[[pos[[1]][[1]]]],00,   $treq=
                     Join[$treq,{{cl,swc,value,$curinst[[pos[[1]][[1]]]]}}]]
  ]]]]
              

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




findring[orto_,dim_] := Block[{ch,cjjj,ns,ms1,mon,degrees,degi},
       $Js=Table[J[i],{i,1,$m}];
       $qs=Table[q[i],{i,1,$m}];
       ort=orto/.{o[i_]:>J[i]};
       degrees=Map[#1-1&,Flatten[Array[a,Table[dim+1,{i,1,$m}]]]/.a->List];
       Do[
         free$ring[i]=Map[xlinamo[#1]&,Select[degrees,(Plus @@ #1)==i&]],
       {i,0,dim}];
       orth3={};
       Do[
         degi=getdeg[ort[[i]]];
         Print[degi];
         fr=free$ring[dim-degi];
         Do[
           orth3=Join[orth3,{convpolyvec[fr[[j]]*ort[[i]],dim]}],
         {j,1,Length[fr]}],
       {i,1,Length[ort]}];
      ns1=NullSpace[orth3];
      ns=Sum[ns1[[i]]*y[i],{i,1,Length[ns1]}];     
      $ring=convvecpoly[Flatten[ns],dim];
      ch=Expand[$ring+xxdum][[2]];
      cjjj=ch/.J[i__]:>1;
      mon=ch/cjjj;
      Print[$ring];
      Print["Normalisation of the ",mon," intersection required"];
      norm=Input[];
      $ring=Expand[(norm)*$ring/cjjj];
      tr[$ring,dim]]

getdeg[poly_]:=Blxock[{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[free$ring[deg]];
           npoly=Expand[poly];
           vec=Table[Coefficient[npoly,free$ring[deg][[ii]]],{ii,1,dimv}]]

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



tr[ringg_,dim_]:= Block[{hi,hlist},
          hlist=Map[#1-1&,Flatten[Array[a,Table[4,{i,1,$m}]]]/.a->List];
          If[$lflag==0,
          deg2=Map[xlinamo[#1]&,Select[hlist,(Plus @@ #1)==dim-1&]],
          deg2=Map[xlinamo[#1]&,Select[hlist,(Plus @@ #1)==dim&]]];
          kopl=Map[xlinamo[#1]&,Select[hlist,(Plus @@ #1)==dim&]];
          (*Print[Head[ringg]]*); 
          If[Head[ringg] === Plus,
          topring=Sum[ringg[[i]]/Apply[Times,Map[#!&,xmonali[ringg[[i]]]]],
                   {i,1,Length[ringg]}],
          topring=ringg/Apply[Times,Map[#!&,xmonali[ringg]]]];
          If[$lflag==1,
              pv[1]=topring;$mc=1
             ,   
              Do[pv[i]=D[topring,J[i]],{i,1,$m}];$mc=$m]]







checktop := Block[{hli},
            li0=Map[Drop[#1,-1][[1]]&,$l];
            li1=Map[Drop[#1,1][[1]]&,$l];
            Clear[clk];
            clk[i_, j_, k_] := clk[Apply[Sequence, Sort[{j, i, k}]]]; 
            Do[
              clk[Sequence @@ convlist[xmonali[kopl[[j]]]]]=
              Coefficient[Expand[$ring/.{kopl[[j]]->X}],X],
            {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}];
       Write[$strf,"*********************************************"];
       Save[$strf,ceul];
       Write[$strf,"*********************************************"];
       Save[$strf,c2J];
       Write[$strf,"*********************************************"]]



(* Change of basis: Let J  be the old and JJ the new basis for the 
   K"ahlercone and $ma be the matrix  such that JJ_i=ma_{i,k} J_k. *)


ma1={{1,0,0},{0,1,1/2},{0,0,1}}

ma2={{1,0,0},{0,1,0},{0,1,1}}

chbask[ma_]:=Block[{hi,qq},
            $m=Length[ma];
            mat=Transpose[ma]; 
            chrule=Table[q[i]->Product[qq[j]^mat[[i,j]],{j,1,$m}],{i,1,$m}];
            Clear[nykold];
            nyk[i_, j_, k_] := nyk[Apply[Sequence, Sort[{j, i, k}]]]; 
            Do[
               nykold[i,j,k]=nyk[i,j,k],
            {k,1,$m},{j,1,$m},{i,1,$m}];
            Clear[nyk];
            Do[nyk[i,j,k]=
               Expand[(Sum[mat[[i1,i]]*mat[[i2,j]]*mat[[i3,k]]*nykold[i1,i2,i3],
                   {i1,1,$m},{i2,1,$m},{i3,1,$m}]/.chrule)/.{qq[i__]:>h q[i]}],
             {k,1,$m},{j,1,k},{i,1,j}]]


(* Change of basis: Let t  be the old and tt the new basis for the 
   Moricone and mb be the matrix  such that tt_i=mb_{i,k} t_k. *)


ma1={{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{1,0,2,1,1},{0,0,0,0,1}}
ma2={{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{3,0,3,1,0},{0,0,0,0,1}}
ma3={{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,1,0,1,0},{0,0,0,0,1}}

chbasm[ma_,nofnewbas_]:=Block[{hi,qq},
            $m=Length[ma];
            mat=Inverse[ma]; 
            chrule=Table[q[i]->Product[qq[j]^mat[[i,j]],{j,1,$m}],{i,1,$m}];
            lm[nofnewbas]=Table[Sum[ma[[i,j]]*$l[[j]],{j,1,$m}],{i,1,$m}];
            Clear[nykold];
            nyk[i_, j_, k_] := nyk[Apply[Sequence, Sort[{j, i, k}]]]; 
            Do[
               nykold[i,j,k]=nyk[i,j,k],
            {k,1,$m},{j,1,$m},{i,1,$m}];
            Clear[nyk];
            ring[nofnewbas]=0; 
            Do[nyk[i,j,k]=
               Expand[(Sum[mat[[i1,i]]*mat[[i2,j]]*mat[[i3,k]]*nykold[i1,i2,i3],
                   {i1,1,$m},{i2,1,$m},{i3,1,$m}]/.chrule)/.{qq[i__]:>h q[i]}];
               ring[nofnewbas]=ring[nofnewbas]+
                J[i]*J[j]*J[k]*(nyk[i,j,k]/.{q[a__]:>0}),
             {k,1,$m},{j,1,k},{i,1,j}]]



   count[slots_,prec_]:=Block[{hi},
                       Do[ip[i]=0,{i,1,slots}];
                       isum:=0;
                       i=1;
                       While[i<=slots,
                        If[isum<prec,
                           ip[i]=ip[i]+1;
                           Print[Table[ip[j],{j,1,slots}]]; 
                           isum=isum+1;
                           i=1,
                           isum=isum-ip[i];
                           ip[i]=0;
                           i=i+1]]]

                           


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

topdata[nm_,ndim_] := Block[{h,vol,red,amb,nor,hip,ll},
   Clear[J];
   $l=lm[nm];   
   $m=Length[$l];
   $Js=Table[J[i],{i,1,$m}];
   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)==ndim&]];
   vol=Product[Sum[dd[[i,j]] h J[j],{j,1,Length[dd[[1]]]}],
         {i,1,Length[dd]}];
   If[vol==0,vol=1,00]; 
   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[j],{j,1,Length[dd[[1]]]}]),
                  {i,1,Length[dd]}];
   Print["amb: ", amb,"  nor:  ", nor,"  vol  ",vol,"  red:  ",red];
   hip=xlinamo[Table[Length[ww[[i]]]-1,{i,1,Length[ww]}]];
   Print["coefficientof ",hip];
   ll=Length[Flatten[ww]]-$m;
   c3=Coefficient[Expand[Normal[Series[amb/nor,{h,0,ndim}]]],h^ndim];
   Print[c3];
   eul=Coefficient[Expand[(Normal[Series[ h^ndim*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];
   Print[c2];
   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=Transpose[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];

degz/: degz[p_]:=Sum[Exponent[p,z[ix]],{ix,1,$m}]



(* Defines the ring of a a local subsystem of a global CY 
   if t[li[[i]]] go to infinity *)
limit[ring_,li_]:=Block[{hi},
                tr[ring,3];
                hi={};
        Do[hi=Union[hi,Select[deg2,MemberQ[#1,J[li[[i]]]]&]],{i,1,Length[li]}];
    rep=Table[hi[[i]]->H[i],{i,1,Length[hi]}];
                ss=Expand[Sum[x[i]D[topring,J[i]]/.rep/.J[a__]:>0,{i,1,$m}]];
   eqs=Table[Coefficient[ss,H[i]]==0,{i,1,Length[rep]}];               
                l3=Flatten[Solve[eqs,Table[x[i],{i,1,$m}]]];
            res=Expand[Sum[x[i] D[topring,J[i]]/.rep/.l3/.H[i__]:>0,{i,1,$m}]]]



(* =====================        N E W   I N V E R T     ==================== *)

IinvTime:=LASTinvTIME=AbsoluteTime[];				(* Init *)
RinvTime:=Block[{rtime=Round[AbsoluteTime[]-LASTinvTIME]},	(* Real *)
             LASTinvTIME=AbsoluteTime[];Return[rtime]];
IinsTime:=LASTinsTIME=AbsoluteTime[];				(* Init *)
RinsTime:=Block[{rtime=Round[AbsoluteTime[]-LASTinsTIME]},	(* Real *)
             LASTinsTIME=AbsoluteTime[];Return[rtime]];

Do[z2qZXlist[a]=List[], {a,1,6}];	(* global variables for z2q-info *)
z2qPlist=List[];			(* as computed in psexpandinvert *)

psexpandmakeinsertlist[VS_,H_,N_]:=		(* VSli[[x]]:=VS^x to O(N-x) *)
  Block[{x,a,b,VSli},VSli=List[Collect[VS,H,Expand]]; Do[AppendTo[VSli,Sum[H^b
    *Expand[Sum[ Coefficient[VSli[[1]],H,a]*Coefficient[VSli[[x-1]],H,b-a]
		  ,{a,0,b}]] ,{b,0,N-x}]] ,{x,2,N}]; Return[VSli]]

psexpandsolve[Q_,QS_,Z_,H_,N_]:=    (* solve Q=QS[Z] for Z to order H^N *)
  Block[{a,b,Zli,Zz,Tz,czli,cle}, Zz=Series[Q-QS+Z,{H,0,N-1}];
    (* *) Tz=Series[Coefficient[Normal[Zz],Z],{H,0,N-1}];    (* improve *)
    Zz=Map[Expand,Map[Expand,Zz-Z*Tz] / Map[Expand,(1-Tz)]];
    (* Zz=Map[Expand,Map[Expand,Zz-Z*Tz]/(1-Tz)]; *)  (* end of improve *)  
    Zli=psexpandmakeinsertlist[Normal[Zz],H,N];
    Zz=SeriesCoefficient[Zli[[1]],0];
    If[Zz!=Q,Print["QS != Z H +O(H^2) ... bad equation"];Abort[]]; Zz=Zli[[1]];
    Do[czli=Coefficient[Zz,H,a]; Zz=Collect[Zz-czli*H^a,H]; 
      czli=CoefficientList[czli,Z]; cle=Length[czli]; 
      Zz=Collect[Zz+czli[[1]]*H^a,H];  			Do[Zz=Collect[Zz+
	Collect[H^a*czli[[b+1]]*Normal[Series[Zli[[b]],{H,0,N-a-1}]],H,Expand]
	,H,Expand]	,{b,1,cle-1}] ,{a,1,N-1}];	Return[Zz]]; 

psexpandinsertlist[V_,Vli_,X_,H_,N_]:= 	(* return X /. V^a->Vli[[a]] *)
  Block[{a,b,XX,coli,cole},		       (* X, VS are order 0 thru N-1 *)
    XX=Coefficient[X,V,0];
    Do[coli=CoefficientList[ Coefficient[X,H,a], V];
       cole=Length[coli];		(* add  coli[[b+1]] = coef(V^b H^a) *)
       Do[ XX=Collect[XX+ Collect[H^a*Normal[Series[Vli[[b]],{H,0,N-a-1}]]
	      *coli[[b+1]] ,H,Expand],H] ,{b,1,cole-1}]	(* 1 <= b <= a+1 *)
    			,{a,N-1,0,-1}]; Return[XX]]     (* 0 <= a <= N-1 *)

psexpandinsert[V_,VS_,X_,H_,N_]:=   	      (*	return X /. V->VS    *)
  Return[psexpandinsertlist[V,psexpandmakeinsertlist[VS,H,N],X,H,N]];

psexpandinvert[Q_,QS_,Z_,ZS_,M_,H_,N_]:=	(* solve Q[i]=QS[Z] for Z[j] *)
  Block[{j,k,J,K,pi=List[],leli=List[],Zli}, IinvTime;		   (* i,j<=M *)
    Do[If[0!=QS[j]/.{H->0}, Print["QS not of order H !!!"]; (* test *)
      Abort[]], {j,1,M}];  (* next: permutation for solving short eqs. first *)
    Do[AppendTo[leli,Length[Expand[QS[k]]]]; AppendTo[pi,k], {k,1,M}];	      
    Do[If[leli[[pi[[k]]]]>leli[[pi[[j]]]], K=pi[[k]];pi[[k]]=pi[[j]];
      pi[[j]]=K],{k,1,M-1},{j,k+1,M}]; Print["|qs|=",leli," pi=",pi];
    z2qPlist=pi;
    Do[ZS[k]=QS[k]/H,{k,1,M}]; Do[ K=pi[[k]];                     (* QS /= H *)
      ZS[K]=psexpandsolve[Q[K],ZS[K],Z[K],H,N]; IinsTime;
      z2qZXlist[K]=Zli=psexpandmakeinsertlist[ZS[K],H,N];
      Do[J=pi[[j]]; If[J!=K, ZS[J]=psexpandinsertlist[Z[K],Zli,ZS[J],H,N]]
      ,{j,1,M}];Print["z[",K,"] ",RinvTime,"sec (",RinsTime,"ins) ",MBmem,"MB"]
    ,{k,1,M} ]; Do[ZS[j]=Expand[ZS[j]*H],{j,1,M}];		  (* ZS *= H *)
    leli=List[];Do[AppendTo[leli,Length[ZS[j]]],{j,1,M}];Print["|zn|=",leli]]

psexpandinsertZ2Q[Z_,X_,M_,H_,N_]:=	(* X /. Z[i]^a -> z2qZXlist[i][[a]] *)
  Block[{a,A,XX},XX=X; Do[A=z2qPlist[[a]]; 
    XX=psexpandinsertlist[Z[A],z2qZXlist[A],XX,H,N],{a,1,M}];Return[XX]];


insertZ2Q[X_]:=psexpandinsertZ2Q[z,X,$m,h,$no];(* use z2qZXlist for z -> zn *)

newinvert:=(psexpandinvert[q,qs,z,zn,$m,h,$no];
	rzq=Table[z[i]->sexpand[Series[zn[i],{h,0,$no}]],{i,1,$m}])

invert:=newinvert



(*	=================	D5 data and testing	  =============      *)

lm[d5]={{{-2,0},{0,0,0,0,0,2,0,0,0,-1,1}},{{0,-2},{0,0,0,0,2,0,0,0,-1,1,
        0}},{{0,0},{1,1,-1,0,0,0,0,0,0,0,-1}},{{0,0},{0,0,1,1,0,0,0,0,0,
        0,-2}},{{0,0},{0,0,0,0,-2,0,0,1,1,0,0}},{{0,0},{0,0,0,0,1,-1,1,-1,0,0,
        0}}}
ring[d5]=8J[1]^3+16J[1]^2J[2]+2J[1]^2J[3]+3J[1]^2J[4]+24J[1]^2J[5]+16J[1]^2J[
        6]+16J[1]J[2]^2+4J[1]J[2]J[3]+6J[1]J[2]J[4]+24J[1]J[2]J[5]+16J[1]J[
        2]J[6]+J[1]J[3]J[4]+6J[1]J[3]J[5]+4J[1]J[3]J[6]+
    J[1]J[4]^2+9J[1]J[4]J[5]+6J[1]J[4]J[6]+24J[1]J[5]^2+16J[1]J[5]J[
        6]+16J[2]^3+4J[2]^2J[3]+6J[2]^2J[4]+24J[2]^2J[5]+16J[2]^2J[6]+2J[2]J[
        3]J[4]+6J[2]J[3]J[5]+4J[2]J[3]J[6]+2J[2]J[4]^2+9J[2]J[4]J[5]+6J[2]J[
        4]J[6]+24J[2]J[5]^2+16J[2]J[5]J[6]+3J[3]J[4]J[5]+2J[3]J[4]J[6]+6J[
        3]J[5]^2+4J[3]J[5]J[6]+3J[4]^2J[5]+2J[4]^2J[6]+9J[4]J[5]^2+6J[4]J[5]J[
        6]+24J[5]^3+16J[5]^2J[6]

OLDpsexpandinversiontest[M_,H_,N_]:=
  Block[{a,b,x}, RTime; Do[QZQ[a]=Normal[Series[qs[a]/H,{H,0,N-1}]]; 
  Do[QZQ[a]=psexpandinsert[z[b],zn[b]/H,QZQ[a],H,N],{b,M,1,-1}];
  x=Series[QZQ[a]-q[a]+H^N,{H,0,N-1}][[4]];
  If[x<N, Print["q[",a,"]   ERROR  at  order   ",x+1," ... ",RTime,"sec"],
     Print["q[",a,"] correct ... ",RTime,"sec"]]       ,{a,1,M}]]

psexpandinversiontest[M_,H_,N_]:=
  Block[{a,b,x}, RTime; Do[QZQ[a]=Normal[Series[qs[a]/H,{H,0,N-1}]]; 
  QZQ[a]=insertZ2Q[QZQ[a]];
  x=Series[QZQ[a]-q[a]+H^N,{H,0,N-1}][[4]];
  If[x<N, Print["q[",a,"]   ERROR  at  order   ",x+1," ... ",RTime,"sec"],
     Print["q[",a,"] correct ... ",RTime,"sec"]]       ,{a,1,M}]]

inversiontest:=psexpandinversiontest[$m,h,$no]

comparetest:=(invert:=oldinvert; Print[Timing[cases[d5,NN,"D5.inst.old"]]];
	Do[auxZN[a]=zn[a],{a,1,$m}]; newinvert; 
	Do[Print[If[auxZN[a]==zn[a],True,False,False]], {a,1,$m}];)

invtest:=(Timing[newinvert;]; inversiontest)

d5old:=(invert:=oldinvert; Print[Timing[cases[d5,NN,"D5.inst.old"]]]);

d5new:=(invert:=newinvert; Print[Timing[cases[d5,NN,"D5.inst.new"]]]);










