AllLoopSunset/dimension-raising.mpl at main · pierrevanhove/AllLoopSunset · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
# This code computes the dimensions raising operator for the all equal mass case sunset integral
with(DEtools): with(PDETools):

# converting epsilon to d
epsilon:=(2-d)/2;

# debug
printdebug:=false;

# this file is from [arXiv:2401.09908] and https://github.com/pierrevanhove/TwistedGriffithsDwork#readme
read("PFSunsetEqualMass-1to20-epsilon.txt"):

# Put the differential operators in a list
PFSunsetEqualMass:=[PFEqualMass1,
PFEqualMass2,
PFEqualMass3,
PFEqualMass4,
PFEqualMass5,
PFEqualMass6,
PFEqualMass7,
PFEqualMass8,
PFEqualMass9,
PFEqualMass10,
PFEqualMass11,
PFEqualMass12,
PFEqualMass13,
PFEqualMass14,
PFEqualMass15,
PFEqualMass16,
PFEqualMass17,
PFEqualMass18,
PFEqualMass19,
PFEqualMass20];

Operator_O:=proc(loop::integer)
   global printdebug;
   if printdebug then print("Generate the coefficients of the operator O") end if;
   return sum(a[r](t)*Dt^r,r=0..loop-1)
end:

Operator_LO:=proc(loop::integer)
   local Ldp2;
   global printdebug;
   if printdebug then print("Generate the operator L(d+2)*O") end if;
   Ldp2:=subs(d=d+2,PFSunsetEqualMass[loop]);
   return mult(Ldp2,Operator_O(loop),[Dt,t]);
end:

Pol_gen:=proc(k::integer,c):
   sum(c[r]*t^r,r=0..k)
end:

Source_dimraising:=proc(loop::integer):
   Pol_gen(loop-1,SC);
end:

Reduce_Order:=proc(loop::integer,r::integer):
   local Ltmp,fLtmp,listtmp,itmp,fstart,fredtmp;
   global printdebug;
   if printdebug then print("Reduce derivative order", r," modulo PF equation") end if;
   rightdivision(Dt^r,PFSunsetEqualMass[loop],[Dt,t]);
end:

Reduce_LO:=proc(loop::integer,Source):
   local LOtmp,ordertmp,redtmp,subtmp,itmp,sourcetmp1,sourcetmp2,Ldp2;
   global printdebug;
   if printdebug then print("Reduce the operator L(d+2)O") end if;
   LOtmp:=Operator_LO(loop);
   redtmp:=rightdivision(LOtmp,PFSunsetEqualMass[loop],[Dt,t]);
   sourcetmp1:=expand(subs(f(t)=Source,diffop2de(redtmp[1],f(t),[Dt,t])));
   Ldp2:=diffop2de(subs(d=d+2,PFSunsetEqualMass[loop]),f(t),[Dt,t]);
   sourcetmp2:=expand(subs(f(t)=Source_dimraising(loop),Ldp2));
   [redtmp[2],sourcetmp1+sourcetmp2];
end:

System_A:=proc(loop::integer,Source):
   local Ldtmp,LOtmp;
   global printdebug;
   LOtmp:=Reduce_LO(loop,Source);
   if printdebug then print("Extract the coefficients") end if;
   [coeffs(collect(LOtmp[1],Dt),Dt),Source-LOtmp[2]];
end:

System_C:=proc(loop::integer,Source):
   local systmp,substmp,itmp,syssubtmp,ptmp,dentmp,syssubtmpnum;
   global printdebug;
   systmp:=System_A(loop,Source);
   substmp:=[seq(a[itmp](t)=Pol_gen(loop+itmp,c[itmp]),itmp=0..loop-1)];
   syssubtmpnum:=numer(expand(subs(substmp,systmp)));
   if printdebug then print("build the system of equations") end if;
   seq(coeffs(collect(ptmp,t),t),ptmp in syssubtmpnum);
end:

Solve_DimRaising:=proc(loop::integer,S):
   local sysCtmp,vartmp,soltmp,listatmp,itmp,listRtmp;
   sysCtmp:=System_C(loop,S);
   vartmp:=indets([sysCtmp])minus{d,S};
   soltmp:=solve([sysCtmp],vartmp);
   listatmp:=collect(subs(soltmp,[seq(a[itmp](t)=Pol_gen(loop+itmp,c[itmp]),itmp=0..loop-1)]),t,factor);
   listRtmp:=collect(subs(soltmp,Source_dimraising(loop)),t,factor);
   [op(listatmp),R[loop](t)=listRtmp];
end:

integral_4d:=proc(loop::integer,Source)
   local dimension_raising_tmp, operator_finite_4d, pole_part,optmp,optmpser,sourcetmp,r;
   dimension_raising_tmp:=Solve_DimRaising(loop,Source);
   optmp:=[seq(subs(dimension_raising_tmp,a[r](t)*(2/(2-d))^loop),r=0..loop-1)];
   optmpser:=[seq(convert(series(optmp[r],d=2,loop+1),polynom),r=1..loop)];
   sourcetmp:=-1*(loop+1)!*subs(dimension_raising_tmp,R[loop](t))*(GAMMA((2-d)/2))^loop*exp(-1/2*(d-2)*loop*gamma);
   operator_finite_4d:=collect(sum(optmpser[r]*Dt^(r-1),r=1..loop),Dt,factor);
   pole_part:=collect(subs(d=2-2*eps,convert(series(sourcetmp,d=2,loop+1),polynom)),eps,simplify);
   return [operator_finite_4d,pole_part]
end: