(* ::Package:: *)

(* ::Subtitle:: *)
(*Restriction z1 \[RightArrow] 0 with logarithmic corrections*)


SetDirectory @ NotebookDirectory[];
Import["../utils.m"]


(* ::Section:: *)
(*Pfaffian matrices*)


(* ::Text:: *)
(*Pfaffian matrices for the massive Bhabha box*)


{P1, P2} = Import["massive_pfaffians.m"];


(* ::Text:: *)
(*P1 not in normal form: 2nd order pole (want 1st order)*)


P1 // mseries1[z1, 0]


(* ::Text:: *)
(*P2 not in normal form: 1st order pole (want finite)*)


P2 // mseries1[z1, 0]


(* ::Section:: *)
(*Normal form by Moser reduction*)


(* ::Text:: *)
(*Shearing gauge transformation*)


G1 = {z1, 1, 1, 1, 1} // DiagonalMatrix


(* ::Text:: *)
(*New P1 has 1st order pole at z1=0*)


P1G1 = P1 // gauge[G1, z1] // Factor;
P1G1 // mseries1[z1, 0]


(* ::Text:: *)
(*New P2 is finite at z1=0*)


P2G1 = P2 // gauge[G1, z2] // Factor;
P2G1 // mseries1[z1, 0]


(* ::Text:: *)
(*Residue matrix for new P1*)


P1G1r = P1G1 // res[z1]


(* ::Text:: *)
(*Spectrum is non-resonant*)


P1G1r // Eigenvalues


(* ::Text:: *)
(*Check that integrability still holds*)


integ[P1G1, P2G1][z1, z2]


(* ::Section:: *)
(*Restriction: eigenvalue 0*)


(* ::Text:: *)
(*We work with the rank-3 Jordan block corresponding to eigenvalue 0*)


P1G1r // jord


(* ::Text:: *)
(*Three master integrals for the massless system*)


masters0 = Import["massless_master_integrals.m"] /. j[_,\[Nu]__] :> j[\[Nu]]


(* ::Text:: *)
(*Five master integrals for the massive system*)


mastersm = Import["massive_master_integrals.m"] /. j[_,\[Nu]__] :> j[\[Nu]]


(* ::Text:: *)
(*The following constraints hold.*)
(*They tell us that the massless tadpole j[0,0,1,0] vanishes, and there is an IBP between  the massless triangle j[0,1,1,1] and the massless buble j[0,1,0,1]*)


constraint = P1G1r . mastersm == 0 // Thread // DeleteCases[#, True]&;
constraint // Column[#, Spacings -> 2]&
constraintSol = Solve[constraint, {j[0,0,1,0], j[0,1,1,1]}] // Flatten // ReplaceAll


(* ::Text:: *)
(*masters0 = B . mastersm*)


B = CoefficientArrays[masters0, mastersm] [[2]] // Normal;
B // MatrixForm


(* ::Text:: *)
(*R-matrix*)


R = P1G1r // rowReduce


(* ::Text:: *)
(*M-matrix*)


M = Join[B, R]


(* ::Text:: *)
(*It has the following property*)


M . mastersm // constraintSol // Factor


(* ::Text:: *)
(*Restriction*)


Q2 = (D[M,z2] + M . (P2G1 /. z1 -> 0)) . Inverse[M] // RightComposition[
	#[[;;3 , ;;3]]&,
	Factor,
	Identity
];


(* ::Text:: *)
(*This is the Pfaffian matrix for the basis masters0 in the variable z2 = t/s*)


Q2 // MatrixForm


(* ::Text:: *)
(*Export the holomorphic solution vector (this is eq. (5.51) in the paper)*)


Isol = G1 . Inverse[M] . Join[masters0, {0, 0}]
Export["holomorphic_solution_vector.m", Isol]


(* ::Section:: *)
(*Restriction: eigenvalue -\[Epsilon]*)


(* ::Text:: *)
(*We need to shift the residue matrix P1G1r by \[Lambda]2 a few times, so we define this diagonal matrix*)


\[Lambda]2 = -\[Epsilon] * IdentityMatrix[5]


(* ::Subsection:: *)
(*Resolve constraint (P1G2r - \[Lambda]2)^2 . I = 0*)


(* ::Text:: *)
(*Solution vector called I^(\[Lambda]2,0) in the text*)


I0 = {i1[z2], i2[z2], i3[z2], i4[z2], i5[z2]}


(* ::Text:: *)
(*R-matrix from the square of the residue matrix*)


R = With[{P = P1G1r - \[Lambda]2}, P . P // rowReduce]
R // MatrixForm


(* ::Text:: *)
(*Constraints imposed by R*)


constraint = R . I0 // RightComposition[
	Factor,
	Map @ Numerator,
	Equal[#, 0]&,
	Thread
];

constraint // Column


(* ::Text:: *)
(*We solve the constraint for, say, {i2, i3, i5}*)


sol = Solve[constraint, {i2[z2], i3[z2], i5[z2]}] // Flatten


(* ::Text:: *)
(*New solution vector which solves R . I0 = 0*)


I0 = I0 //. sol // FullSimplify
R . I0 // Factor


(* ::Text:: *)
(*The solution vector called I^(\[Lambda]2,1)*)


I1 = (P1G1r - \[Lambda]2) . I0 // FullSimplify


(* ::Subsection:: *)
(*Pfaffian system*)


(* ::Text:: *)
(*There are two indepedent functions i1[z2] and i4[z2] in I0.*)
(*Let's call them J1[z1] and J2[z2] as in the paper.*)


toJ = {i1 -> J1, i4 -> J2} // ReplaceAll;
J = {i1[z2], i4[z2]} // toJ


(* ::Text:: *)
(*A nice choice of B.*)
(*We could choose anything we want, as long as M becomes invertible.*)


B = {
	{0, 0, 0, \[Epsilon]*z2 / (\[Epsilon]-1), \[Epsilon]*z2 / (2(\[Epsilon]-1))}, 
	{0, 0, 0, 1, 0}
}


(* ::Text:: *)
(*M-matrix*)


M = Join[B, R]


(* ::Text:: *)
(*M has this property*)


M . I0 // toJ // Factor


(* ::Text:: *)
(*Q2 is the Pfaffian matrix for the basis J*)


Q2 = (D[M,z2] + M . (P2G1 /. z1 -> 0)) . Inverse[M] // RightComposition[
	#[[;;2 , ;;2]]&,
	Factor,
	Identity
];

Q2 // MatrixForm // Magnify


(* ::Subsection:: *)
(*Solve the Pfaffian system*)


(* ::Text:: *)
(*Pfaffian system*)


system = D[J, z2] - Q2 . J == 0 // Thread


(* ::Text:: *)
(*It is simple enough that it can be solved exactly by Mathematica*)


systemsol = DSolve[system, J, z2] // Flatten // ReplaceAll


(* ::Text:: *)
(*Now we can write the full solution vectors*)


{I0sol, I1sol} = {I0, I1} // toJ // systemsol


(* ::Text:: *)
(*The final result includes our gauge transformation G1 and a factor of z1^\[Lambda]2*)


Isol =  z1^(-\[Epsilon]) * G1 . (I0sol + I1sol * Log[z1]) // FullSimplify // ReplaceAll[Log[z1]-Log[z2] -> Log[z1/z2]]


(* ::Text:: *)
(*We still need to solve for c1[\[Epsilon]] and c2[\[Epsilon]]. *)
(*This is done in the next file, so we export Isol*)


Export["log_solution_vector.m", Isol]
