(* ::Package:: *)

(* Canonical basis and diff. eq. from 1307.4083 *)
(* -s/m^2 = (1-x)^2/x, -t/m^2 = (1-y)^2/y *)
(* resolves roots Sqrt[(-s)*(4*m^2-s)] = m^2*Sqrt[(1-x)^2/x*(4+(1-x)^2/x)] = m*(1-x)*(1+x)/x *)
(* x=0 \[TwoWayRule] s=Infinity, x=1 \[TwoWayRule] s=0, x=-1 \[TwoWayRule] s=4*m^2, x=-y \[TwoWayRule] u=0*)
(* y=0 \[TwoWayRule] t=Infinity, y=1 \[TwoWayRule] t=0, y=-1 \[TwoWayRule] t=4*m^2 *)


(* ::Subsection::Closed:: *)
(*Basis integrals and boundary conditions*)


(* f1 = (m^2)^eps*Exp[eps*EulerGamma]*eps*G2000 *)
f1 = eps*Gamma[eps]*Exp[eps*EulerGamma];
(* f2 = (m^2)^eps*Exp[eps*EulerGamma]*eps*t*G0201 *)
f2 = -eps*Gamma[1-eps]*Gamma[-eps]*Gamma[1+eps]/Gamma[1-2*eps]*(y/(1-y)^2)^eps*Exp[eps*EulerGamma];


f1exp = FullSimplify[Series[f1,{eps,0,3}]]
f2exp = FullSimplify[Series[f2,{eps,0,3}]]/.{Log[y/(-1+y)^2]->Log[y]-2Log[1-y]}


(* f3 = (m^2)^eps*Exp[eps*EulerGamma]*eps*Sqrt[(-s)*(4m^2-s)]*G2010 *)
(* For s=0 : G2010 = fct[eps]/(m^2)^(1+eps). Thus we have f3(x=1,y) = 0 *)
f3BC = 0;

(* Full result: *)
f3 = eps*Gamma[1+eps]*Exp[eps*EulerGamma]/2*(1-x)*(1+x)/x*Hypergeometric2F1[1,1+eps,3/2,-(1-x)^2/(4x)];


(* f4 = (m^2)^eps*Exp[eps*EulerGamma]*(-2*eps^2)*Sqrt[(4m^2-t)*(-t)]*G1101 *)
(* Root is missing in 1307.4083 *)

(* Determine leading term of f4 as expansion around y = 1, factor out (1-y)^(-2*eps) for expansion. *)
(* Leading expansion term of G1101 was determined in Exercises/MB_MassiveTriangleExpansion.m *)
G1101 = (m^2)^(-1-eps) (Sqrt[\[Pi]] Gamma[1/2-eps]^2 Gamma[1/2+eps])/(2 Gamma[1-2 eps]) ((1-y)^2/y)^(-(1/2)-eps);
FullSimplify[FullSimplify[(1-y)^(2 eps) (-2*eps^2*Exp[eps*EulerGamma]*Sqrt[(4+(1-y)^2/y)*(1-y)^2/y]*G1101)/.m->1,Assumptions->0<y<1] /. {(-2+1/y+y)^-eps->(1-y)^(-2*eps)*y^(eps)}];
f4BC = (1-y)^(-2 eps)*(%/.{y->1})
f4BCexp = Series[%,{eps,0,3}]


(* f5 = (m^2)^eps*Exp[eps*EulerGamma]*(-2*eps^2)*Sqrt[(-s)*(4m^2-s)]*(-t)*G1111 *)
(* For s=0 : G1111 = fct[y,eps]/(t,m^2)^(2+eps). Thus we have f5(x=1,y) = 0 *)
f5BC = 0;


(* result for f5 from 1307.4083 for checks *)
f5HS = eps*4*HP[{0},x] + eps^2*(4*GP[{0},y]*HP[{0},x] - 8GP[{1},y]*HP[{0},x]) + eps^3*(
  -8*GP[{0},y]*HP[{-1,0},x] + 4*GP[{0},y]*HP[{0,0},x] - 8*HP[{0},x]*GP[{1,0},y] + 16*HP[{0},x]*GP[{1,1},y]
  +4*HP[{0},x]*GP[{-1/x,0},y] - 8*HP[{0},x]*GP[{-1/x,1},y] + 4*HP[{0},x]*GP[{-x,0},y] - 8*HP[{0},x]*GP[{-x,1},y]
  +8*HP[{-1,0},x]*GP[{-1/x},y] + 8*HP[{-1,0},x]*GP[{-x},y] - 4*HP[{0,0},x]*GP[{-1/x},y] - 4*HP[{0,0},x]*GP[{-x},y]
  +4*GP[{-1/x,0,0},y] - 8*GP[{-1/x,0,1},y] - 4*GP[{-x,0,0},y] + 8*GP[{-x,0,1},y] + 8*HP[{-2,0},x] 
  -16*HP[{-1,-1,0},x] + 8*HP[{-1,0,0},x] - 4*HP[{0,0,0},x] + 10*Pi^2/3*GP[{-1/x},y] - 2*Pi^2*GP[{-x},y]
  -2*Pi^2/3*GP[{0},y] - 4*Pi^2/3*HP[{-1},x] - 7*Pi^2/3*HP[{0},x] + 8*Zeta[3]) /. {HP[{-2,0},x]->GP[{0,-1,0},x]} /. {
  HP[ind_List,x]->GP[ind,x]};


(* BCs at x=y=1 with Lnybar = Log[1-y]*)
fVecBC11exp = {f1exp,f2exp,f3BC,f4BCexp,f5BC} /. {Log[1-y]->Lnybar,y->1};


(* ::Subsection::Closed:: *)
(*Differential equations*)


(* differential equation: df = eps*dAtilde*f *)
alphabet = {x,1+x,y,1+y,1-y,x+y,1+x*y};

(* Atilde = Sum[ai*Log[alphabet[[i]]],{i,i,7}] *)

a1 = Normal[SparseArray[{{3,1}->-1,{3,3}->1,{5,2}->4},{5,5}]];
a2 = Normal[SparseArray[{{3,3}->-2,{5,3}->-8,{5,5}->-2},{5,5}]];
a3 = Normal[SparseArray[{{2,2}->1,{4,1}->2,{4,2}->-2,{5,3}->-4},{5,5}]];
a4 = Normal[SparseArray[{{4,4}->2},{5,5}]];
a5 = Normal[SparseArray[{{2,2}->-2,{4,4}->-2,{5,5}->-2},{5,5}]];
a6 = Normal[SparseArray[{{5,3}->4,{5,4}->2,{5,5}->1},{5,5}]];
a7 = Normal[SparseArray[{{5,3}->4,{5,4}->-2,{5,5}->1},{5,5}]];


(* ::Subsection::Closed:: *)
(*Solution*)


(* Choose straight path from {0,0} to {x,0} and then to {x,y} *)
(* Define Goncharov polylogarithms as GP[{indices},z] *)

fsol[n_,fnMINUS1_List] := Module[{res},
(* path from {0,0} to {x,0} *)
res = (a1+a6).fnMINUS1 * GPNewIndex[0,x] + a2.fnMINUS1 * GPNewIndex[-1,x] //. {y->0,GP[ind_List,0]->0};
(* path from {x,0} to {x,y} *)
res += a3.fnMINUS1 * GPNewIndex[0,y] + a4.fnMINUS1 * GPNewIndex[-1,y] + a5.fnMINUS1 * GPNewIndex[1,y] + 
       a6.fnMINUS1 * GPNewIndex[-x,y] + a7.fnMINUS1 * GPNewIndex[-1/x,y];
(* process GPs *)
res = Expand[res] /. {GPNewIndex[i_,x_]*GP[ind_List,x_] :> GP[Flatten[Prepend[ind,{i}]],x]};
res = res /. {GPNewIndex[i_,x_] :> GP[{i},x]};
(* boundary value *)
res + Table[ToExpression["fVec"<>ToString[n]<>"BC"][i,0,0],{i,1,5}]
];

GPtoLogLi = {GP[{0},x_]:>Log[x],GP[{1},x_]:>Log[1-x],GP[{0,0},x_]:>Log[x]^2/2,GP[{1,0},x_]:>PolyLog[2,x]+Log[x]*Log[1-x],GP[{0,1},x_]:>-PolyLog[2,x],
  GP[{1,1},x_]:>1/2 Log[1-x]^2,GP[{-1,0},x_]:>PolyLog[2,-x]+Log[x] Log[1+x]};

GPSpecialValues = {GP[{0},1]->0,GP[{1},1]->Lnybar,GP[{-1},1]->Log[2],GP[{0,0},1]->0,GP[{1,0},1]->\[Pi]^2/6,GP[{0,1},1]->-(\[Pi]^2/6),GP[{1,1},1]->Lnybar^2/2,GP[{-1,0},1]->-(\[Pi]^2/12),
  GP[{0,0,0},1]->0,GP[{1,0,0},1]->-Zeta[3],GP[{-1,0,0},1]->(3 Zeta[3])/4,GP[{0,1,0},1]->2*Zeta[3],GP[{0,-1,0},1]->-((3 Zeta[3])/2),GP[{0,0,1},1]->-Zeta[3],
  GP[{0,1,1},1]->Zeta[3],GP[{-1,-1,0},1]->-(1/12) \[Pi]^2 Log[2]+Zeta[3]/8,GP[{-1,0,1},1]->-(1/6) \[Pi]^2 Log[2]+(5 Zeta[3])/8,GP[{1,1,0},1]->Lnybar \[Pi]^2/6+Zeta[3],GP[{1,0,1},1]->-Lnybar \[Pi]^2/6-2 Zeta[3],GP[{1,1,1},1]->Lnybar^3/6};
(* with Lnybar = Log[1-y] *)


(* From boundary conditions we obtain at O(eps^0) *)
fVec0 = {1,1,0,0,0};


(* Result at O(eps) *)
fVec1 = fsol[1,fVec0] /. {fVec1BC[i_,0,0]->0}
(* Check BCs *)
Coefficient[fVecBC11exp,eps] - fVec1 /. {x->1,y->1} /. GPSpecialValues
(* Result in terms of Log and PolyLog *)
fVec1Alt = fVec1 /. GPtoLogLi


(* Checks *)
Simplify[fVec1Alt[[1]] - Coefficient[f1exp,eps]]
Simplify[fVec1Alt[[2]] - Coefficient[f2exp,eps]]
FullSimplify[fVec1Alt[[3]] - Coefficient[Series[f3,{eps,0,1}],eps],Assumptions->0<x<1]
Simplify[Coefficient[f5HS,eps,1] - fVec1[[5]]]


(* Result at O(eps^2) *)
fVec2 = fsol[2,fVec1] /. {fVec2BC[1,0,0]->\[Pi]^2/12,fVec2BC[2,0,0]->-\[Pi]^2/12,fVec2BC[3,0,0]->\[Pi]^2/6,fVec2BC[4,0,0]->-4*\[Pi]^2/3,fVec2BC[5,0,0]->0}
(* Check BCs *)
Simplify[Coefficient[fVecBC11exp,eps,2] - fVec2 /. {x->1,y->1} /. GPSpecialValues]
(* Result in terms of Log and PolyLog *)
fVec2Alt = Simplify[fVec2 /. GPtoLogLi]


(* Checks *)
Simplify[fVec2Alt[[1]] - Coefficient[f1exp,eps,2]]
Simplify[fVec2Alt[[2]] - Coefficient[f2exp,eps,2]]
Simplify[Coefficient[f5HS,eps,2] - fVec2[[5]]]


(* Result at O(eps^3) *)
fVec3 = fsol[3,fVec2] /. {fVec3BC[1,0,0]->-Zeta[3]/3,fVec3BC[2,0,0]->-7Zeta[3]/3,fVec3BC[3,0,0]->2Zeta[3],fVec3BC[4,0,0]->-10 Zeta[3],fVec3BC[5,0,0]->8 Zeta[3]}
(* Check BCs *)
Simplify[Coefficient[fVecBC11exp,eps,3] - fVec3 /. {x->1,y->1} /. GPSpecialValues]


(* Checks *)
Simplify[Coefficient[f5HS,eps,3] - fVec3[[5]]]