
accu=7;

accu0 := accu+1;


(* Ceuler = 0; *)

Clear[GG,Gred,z,gaexp,GabExp,Gab];

Unprotect[Ex,Coll,Coef,Exn,Zrule];

Ex[x_]  := Expand[x];

Exn[x_] := Expand[Normal[x]];

Coef[x__] := Coefficient[x];

Coll[x_,y_] := Collect[Ex[x],y];

Zrule = {
Ceuler -> EulerGamma,
z2->Zeta[2], 
z3->Zeta[3],
z5->Zeta[5], 
z4->Zeta[4],
z6->Zeta[6], 
z7->Zeta[7],
z8->Zeta[8], 
z9->Zeta[9] ,

B4->16*PolyLog[4,1/2]+2/3*Log[2]^4-2/3*Pi^2*Log[2]^2-13/180*Pi^4

}       ;


Protect[Ex,Coll,Coef,Exn,Zrule];




Module[{z,epp},
z[2] = z2;
z[3] = z3;
z[4] = z4;
z[5] = z5;
z[6] = z6;
z[7] = z7;
z[8] = z8;
z[9] = z9;
z[10] = z10;
z[11] = z11;
z[12] = z12;
z[13] = z13;
z[14] = z14;
z[15] = z15;
z[16] = TooHigAccuracy;

gaexp[x_,epp_] := (-Ceuler*x*epp + Sum[( -x*epp)^i/i*(z[i]),{i,2,16}] );
gaexp[x_     ] := (-Ceuler*x     + Sum[( -x    )^i/i*(z[i]),{i,2,16}] );

];


ga[n_,a_] :=  (Gred[n,a] * Exn[ Exp[ gaexp[a] ]  + O[ep]^accu0]) ;  

Poch[a_,n_Integer ] := Pochhammer[a,n];

Gred[n_Integer,a_ ] := Gred[n,a] = Poch[1+a,n-1]; 


Block[{accu0=8},
nG   =  Exn[(ga[2,-2*ep]/(ga[1,ep]*ga[1,-ep]^2)) + O[ep]^accu0];
nG2l =  Exn[(ga[2,-2*ep]/(ga[1,ep]*ga[1,-ep]^2))^2 + O[ep]^accu0];
nG3l =  Exn[(ga[2,-2*ep]/(ga[1,ep]*ga[1,-ep]^2))^3 + O[ep]^accu0];
nG4l =  Exn[(ga[2,-2*ep]/(ga[1,ep]*ga[1,-ep]^2))^4 + O[ep]^accu0];


G[1,1] := 1/ep*Exn[((ga[1,ep]*ga[1,-ep]^2)/ga[2,-2*ep]) + O[ep]^accu0];

G0 = ep*G[1,1];


G2MSbar1l = Block[{Ceuler=0}, Exn[G0]];
G2MSbar2l = Block[{Ceuler=0}, Exn[G0^2+  O[ep]^accu0]];
G2MSbar3l = Block[{Ceuler=0}, Exn[G0^3+  O[ep]^accu0]];
G2MSbar4l = Block[{Ceuler=0}, Exn[G0^4+  O[ep]^accu0]];
G2MSbar5l = Block[{Ceuler=0}, Exn[G0^5+  O[ep]^accu0]];

];

GG[a_Integer,a1_Integer,b_Integer,b1_Integer] :=  Module[ {pref,res,accu=accu0},


If[ ((a <= 0) && (a1 === 0) )|| ((b <= 0) && (b1 === 0)), Return[0] ];



pref = Exn[
 Ex[
 Poch[ ep*(a1+b1+1),a+b-2]
 *
 Gred[2 -a , -(a1+1)*ep]
 *
 Gred[2 -b , -(b1+1)*ep]
 /(
  Gred[a,a1*ep]*Gred[b,b1*ep]*Poch[2 -ep*(a1+b1+2),2-a-b]
  )
     ]
     + O[ep]^(accu+2)
    ];



Which[
(Coef[pref,ep,0]) =!= 0,  
Print["pref0 =", pref//Ex];                             
res =Exn[Ex[pref*(ep*GGab[a1,b1])] +O[ep]^(accu+1)];
Nrint["res =", res//Ex];
Return[Exn[1/ep*res    + O[ep]^(accu)] ]
                               ,

(Coef[pref,ep,1]) =!= 0  ,
Print["pref1 =", pref//Ex];
pref = pref/ep; 
res =Exn[pref*(ep*GGab[a1,b1]) +O[ep]^(accu)];
Nrint["res =", res//Ex];
Return[res]
                               ,

(Coef[pref,ep,2]) =!= 0  ,
Print["pref 2=", pref//Ex];
pref =  pref/ep^2; 
res  =  Exn[pref*(ep*GGab[a1,b1]) +O[ep]^(accu-1)];
Nrint["res =", res//Ex];
Return[Expand[ep^2/ep*res]],

True,

Nrint["Something is wrong in the function GG!!"];
Abort[]

];


];

GG[al_,be_] :=  GG[al/.ep->0, Coef[al,ep] ,be/.ep->0, Coef[be,ep]];


Clear[GGab];

GGab[a_,b_] :=  1/((a+b+1)*ep)*Exn[(1-2*ep)/(1-(2+a+b)*ep)*Exp[

GabExp[a,b] ] + O[ep]^(accu0 +1)
                                     ];                                  



GabExp[a_,b_]  = Ex[
+gaexp[a+b+1,ep] - gaexp[1,ep]
+gaexp[-a-1,ep]  - gaexp[-1,ep]
+gaexp[-b-1,ep]  - gaexp[-1,ep]
-(
 gaexp[a,ep]   +  gaexp[b,ep] 
+gaexp[-a-b-2,ep] - gaexp[- 2,ep]
  )
                  ];


(***********  F-function  /from       *************)



Prefactor[a1_,a2_,a3_,a4_,a5_] := Prefactor[a1,a2,a3,a4,a5]=Module[{aa,i,res},
aa[1] = a1;
aa[2] = a2;
aa[3] = a3;
aa[4] = a4;
aa[5] = a5;

aa[6]  = -(a1+a2+a3+a4+a5);
aa[7]  = a1+a5+aa[6];
aa[8]  = a2+a4+aa[6];
aa[9]  = a3+a4+a5;
aa[10] = a1+a2+a3;

res = (Normal[
Product[ ga[1,ep*(-1/2-aa[i])]/ga[1,ep*(aa[i]-1/2)],{i,1,10} ]^(1/2)/(1-2*ep)/ga[1,-ep]^2
+O[ep]^6
             ]);


Return[res]

                                        ];


F0[a1_,a2_,a3_,a4_,a5_] := fg6[a1,a2,a5,a3,a4];


U62 =- 79/806400*(z8*9450) + 3/4*z3*z5 - 9/8*z62;


fg6[a1_,a2_,a3_,a4_,a5_]  := Normal[(
  fb6[a1+1/2, a2+1/2,a3+1/2,a4+1/2,a5+1/2]
 *Prefactor[a1+1/2, a2+1/2,a3+1/2,a4+1/2,a5+1/2]*nG2l
                            + O[ep]^6
                                  )]//Expand;



fb6[a1_,a2_,a3_,a4_,a5_] := Module[{aa,i,res,S2,S4},

aa[1] = a1; 
aa[2] = a2;
aa[3] = a3; 
aa[4] = a4;
aa[5] = a5;

aa[6]  = -(a1+a2+a3+a4+a5); 
aa[7]  = a1+a5+aa[6];
aa[8]  = a2+a4+aa[6]; 
aa[9]  = a3+a4+a5;
aa[10] = a1+a2+a3;

S2 = Sum[aa[i]^2,{i,1,10}]; 
S4 = Sum[aa[i]^4,{i,1,10}]; 

Cb[3,0] = 6*z3;

Cb[3,1] = 9*z4;

Cb[3,2] = 3/4*z5;

Cb[5,2] = 5/2*z5;

Cb[3,3] = -105/8*z6-15/4*z3^2;

Cb[5,3] = 25/4*z6-1/2*z3^2;

Cb[3,4] = -549/16*z7 -45/4*z4*z3;

Cb[5,4] = 29/4*z7-3/2*z4*z3;

Cb[7,4]  = 7/2*z7;

Cd[7,4] =-7/128*z7;

Cb[3,5]  =-147/2*z8-87/8*z5*z3;

Cb[5,5] = 35/16*z8-3*z5*z3;

Cb[7,5] = 49/4*z8-z5*z3;

Cd[7,5]= 551/1024*z8-9/16*z5*z3-1/4*U62;

Cb[x_,y_]  :=0;

Cd[x_,y_]  :=0;

Cb[3] = Sum[Cb[3,i]*ep^i,{i,0,5}];

Cb[5] = Sum[Cb[5,i]*ep^i,{i,0,5}];

Cb[7] = Sum[Cb[7,i]*ep^i,{i,0,5}];

Cd[7] = Sum[Cd[7,i]*ep^i,{i,0,5}];

res = Normal[Cb[3] +Cb[5]*S2+Cb[7]*S4 +Cd[7]*(16*S4-4*S2^2+12*S2-63) +O[ep]^6]

                            ]; 










 (* 

m14= Exn[1/ep^3*GG[ep,2*ep] +O[ep]^6];


m11= Coll[Exn[1/ep^3*GG[ep,ep] +O[ep]^6],ep];


m23= Coll[Exn[1/ep^3*GG[1,1+eep] +O[ep]^5],ep];


 *)











