Bug Summary

File:generators/phokhara/phokhara/eemmg-lib/src/pointer.h
Warning:line 51, column 17
Use of memory after it is freed

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -O3 -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name minor.cpp -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/data/b2soft/buildbot/development/build -fcoverage-compilation-dir=/data/b2soft/buildbot/development/build -resource-dir /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++ -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/x86_64-redhat-linux -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/c++/backward -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/python3.12 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/CLHEP -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/Geant4 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/root -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/include/belle_legacy -I include/ -D _PACKAGE_="generators" -D G4UI_USE_TCSH -D RaveDllExport= -D HAS_SQLITE -D HAS_CALLGRIND -D HAVE_CONFIG_H -I generators/phokhara/phokhara/eemmg-lib/include -I /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/include/libxml2 -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++ -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/x86_64-redhat-linux -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/backward -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/lib/clang/21/include -internal-isystem /usr/local/include -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-00/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../x86_64-redhat-linux/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -Wno-missing-braces -Wno-unused-command-line-argument -std=c++20 -fdeprecated-macro -ferror-limit 19 -fgnuc-version=4.2.1 -fno-implicit-modules -fskip-odr-check-in-gmf -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /scan_build/2026-05-31-004316-385593-1 -x c++ generators/phokhara/phokhara/eemmg-lib/src/minor.cpp

generators/phokhara/phokhara/eemmg-lib/src/minor.cpp

→
1/*
2 * minor.cpp - constructors, signed minors and integrals
3 *
4 * this file is part of PJFry library
5 * Copyright 2011 Valery Yundin
6 */
7
8#include "minor.h"
9#include "cache.h"
10
11const unsigned char MinorBase::ti2[8]={0, 1, 3, 6, 10, 15, 21, 28};
12const unsigned char MinorBase::ti3[8]={0, 1, 4, 10, 20, 35, 56, 84};
13const unsigned char MinorBase::ti4[8]={0, 1, 5, 15, 35, 70, 126, 210};
14const unsigned char MinorBase::ti5[8]={0, 1, 6, 21, 56, 126, 252, 0};//462};
15
16const double MinorBase::teps=1e-14;
17const double MinorBase::heps=1e-15;
18
19const double MinorBase::ceps=5e-2;
20
21const double MinorBase::deps1=5e-2;
22const double MinorBase::deps2=5e-2;
23const double MinorBase::deps3=5e-2;
24
25const double MinorBase::seps1=1e-8;
26const double MinorBase::seps2=1e-5;
27
28const double MinorBase::epsir1=5e-6;
29const double MinorBase::epsir2=5e-6; // m.b. 5e-5 is better
30
31double MinorBase::deps=1e-14;
32double MinorBase::meps=1e-10;
33double MinorBase::m3eps=1;
34
35void MinorBase::Rescale(double factor)
36{
37 meps*=factor;
38 m3eps*=factor*factor;
39}
40
41const unsigned char MinorBase::idxtbl[64]={
42 0, // 0 0b0
43 0, //0 1 0b1
44 1, //1 2 0b10
45 0, //01 3 0b11
46 2, //2 4 0b100
47 1, //02 5 0b101
48 5, //12 6 0b110
49 0, //012 7 0b111
50 3, //3 8 0b1000
51 2, //03 9 0b1001
52 6, //13 10 0b1010
53 1, //013 11 0b1011
54 9, //23 12 0b1100
55 4, //023 13 0b1101
56 10, //123 14 0b1110
57 104, // 15 0b1111
58 4, //4 16 0b10000
59 3, //04 17 0b10001
60 7, //14 18 0b10010
61 2, //014 19 0b10011
62 10, //24 20 0b10100
63 5, //024 21 0b10101
64 11, //124 22 0b10110
65 104, // 23 0b10111
66 12, //34 24 0b11000
67 7, //034 25 0b11001
68 13, //134 26 0b11010
69 104, // 27 0b11011
70 16, //234 28 0b11100
71 104, // 29 0b11101
72 104, // 30 0b11110
73 105, // 31 0b11111
74 5, //5 32 0b100000
75 4, //05 33 0b100001
76 8, //15 34 0b100010
77 3, //015 35 0b100011
78 11, //25 36 0b100100
79 6, //025 37 0b100101
80 12, //125 38 0b100110
81 104, // 39 0b100111
82 13, //35 40 0b101000
83 8, //035 41 0b101001
84 14, //135 42 0b101010
85 104, // 43 0b101011
86 17, //235 44 0b101100
87 104, // 45 0b101101
88 104, // 46 0b101110
89 105, // 47 0b101111
90 14, //45 48 0b110000
91 9, //045 49 0b110001
92 15, //145 50 0b110010
93 104, // 51 0b110011
94 18, //245 52 0b110100
95 104, // 53 0b110101
96 104, // 54 0b110110
97 105, // 55 0b110111
98 19, //345 56 0b111000
99 104, // 57 0b111001
100 104, // 58 0b111010
101 105, // 59 0b111011
102 104, // 60 0b111100
103 105, // 61 0b111101
104 105, // 62 0b111110
105 255, // 63 0b111111
106};
107
108// Find first three indices which are not occupied by set[] and put them in free[]
109void MinorBase::freeidxM3(int set[], int free[])
110{
111 free[0]=0;
112 free[1]=1;
113 free[2]=2;
114
115 int ic=0;
116 int nc=0;
117 while (ic < 3 && nc < 3) {
118 if (free[nc]==set[ic]) {
119 for (int cc=nc; cc<3; cc++) {
120 free[cc]++;
121 }
122 ic++;
123 } else {
124 nc++;
125 }
126 }
127}
128
129/* ------------------------------------------------------------
130* ------------------------------------------------------------
131* Minor2 section
132* ------------------------------------------------------------
133* ------------------------------------------------------------
134*/
135
136// Constructor from higher rank minor
137Minor2::Minor2(const Kinem2& k, Minor5::Ptr mptr5, int s, int t, int u, int is)
138 : kinem(k), pm5(mptr5), ps(s), pt(t), pu(u), offs(is)
139{
140// printf("Minor2 from Minor5\n");
141}
142
143/* ========================================================
144 * ========================================================
145 *
146 * Minor3 section
147 *
148 * ========================================================
149 * ========================================================
150 */
151
152// Constructor from higher rank minor
153Minor3::Minor3(const Kinem3& k, Minor5::Ptr mptr5, int s, int t, int is)
154 : kinem(k), pm5(mptr5), ps(s), pt(t), offs(is)
155{
156// printf("Minor3 from Minor5\n");
157}
158
159/* ========================================================
160 * ========================================================
161 *
162 * Minor4 section
163 *
164 * ========================================================
165 * ========================================================
166 */
167
168// Direct construction is proxied through dummy Minor5
169// see MCache::getMinor4 in cache.cpp
170
171// Constructor from higher rank minor
172Minor4::Minor4(const Kinem4 &k, Minor5::Ptr mptr5, int s, int is)
173 : kinem(k), pm5(mptr5), ps(s), offs(is)
174{
175// printf("Minor4 from Minor5\n");
176}
177
178/* ========================================================
179 * ========================================================
180 *
181 * Minor5 section
182 *
183 * ========================================================
184 * ========================================================
185 */
186
187/* --------------------------------------------------------
188 * Minors with 4 scratched lines
189 * --------------------------------------------------------
190 */
191double Minor5::M4ii(int u, int v, int i)
192{
193 return (Cay[nss(u,u)]*Cay[nss(v,v)]-Cay[nss(u,v)]*Cay[nss(u,v)]);
194}
195
196double Minor5::M4ui(int u, int v, int i)
197{
198 return (Cay[nss(u,v)]*Cay[ns (i,v)]-Cay[ns (i,u)]*Cay[nss(v,v)]);
199}
200
201double Minor5::M4vi(int u, int v, int i)
202{
203 return (Cay[nss(u,v)]*Cay[ns (i,u)]-Cay[ns (i,v)]*Cay[nss(u,u)]);
204}
205
206double Minor5::M4uu(int u, int v, int i)
207{
208 return (Cay[nss(i,i)]*Cay[nss(v,v)]-Cay[ns (i,v)]*Cay[ns (i,v)]);
209}
210
211double Minor5::M4vu(int u, int v, int i)
212{
213 return (Cay[ns (i,v)]*Cay[ns (i,u)]-Cay[ns (u,v)]*Cay[nss(i,i)]);
214}
215
216double Minor5::M4vv(int u, int v, int i)
217{
218 return (Cay[nss(i,i)]*Cay[nss(u,u)]-Cay[ns (i,u)]*Cay[ns (i,u)]);
219}
220
221
222/* --------------------------------------------------------
223 * Preprocessor definitions
224 * --------------------------------------------------------
225 */
226#define k5s1s12,p3,p4,p5,s45,s34,m2,m3,m4,m5 s12,p3,p4,p5,s45,s34,m2,m3,m4,m5
227#define k5s2p1,s23,p4,p5,s45,s15,m1,m3,m4,m5 p1,s23,p4,p5,s45,s15,m1,m3,m4,m5
228#define k5s3p1,p2,s34,p5,s12,s15,m1,m2,m4,m5 p1,p2,s34,p5,s12,s15,m1,m2,m4,m5
229#define k5s4p1,p2,p3,s45,s12,s23,m1,m2,m3,m5 p1,p2,p3,s45,s12,s23,m1,m2,m3,m5
230#define k5s5p2,p3,p4,s15,s23,s34,m2,m3,m4,m1 p2,p3,p4,s15,s23,s34,m2,m3,m4,m1
231
232#define k5st12s45,p4, p5, m3,m4,m5 s45,p4, p5, m3,m4,m5
233#define k5st13s12,s34,p5, m2,m4,m5 s12,s34,p5, m2,m4,m5
234#define k5st14s12,p3, s45,m2,m3,m5 s12,p3, s45,m2,m3,m5
235#define k5st15p3, p4, s34,m3,m4,m2 p3, p4, s34,m3,m4,m2
236#define k5st23p1, s15,p5, m1,m4,m5 p1, s15,p5, m1,m4,m5
237#define k5st24p1, s23,s45,m1,m3,m5 p1, s23,s45,m1,m3,m5
238#define k5st25s23,p4, s15,m3,m4,m1 s23,p4, s15,m3,m4,m1
239#define k5st34p1, p2, s12,m1,m2,m5 p1, p2, s12,m1,m2,m5
240#define k5st35p2, s34,s15,m2,m4,m1 p2, s34,s15,m2,m4,m1
241#define k5st45p2, p3, s23,m2,m3,m1 p2, p3, s23,m2,m3,m1
242
243#define k5stu123p5, m4, m5 p5, m4, m5
244#define k5stu124s45, m3, m5 s45, m3, m5
245#define k5stu125p4, m4, m3 p4, m4, m3
246#define k5stu134s12, m2, m5 s12, m2, m5
247#define k5stu135s34, m4, m2 s34, m4, m2
248#define k5stu145p3, m3, m2 p3, m3, m2
249#define k5stu234p1, m1, m5 p1, m1, m5
250#define k5stu235s15, m4, m1 s15, m4, m1
251#define k5stu245s23, m3, m1 s23, m3, m1
252#define k5stu345p2, m2, m1 p2, m2, m1
253
254#define m5create4(s) \
255{ \
256 Kinem4 k4=Kinem4(k5s##s); \
257 Minor4::Ptr minor=Minor4::create(k4,self,s, offs); \
258 MCache::insertMinor4(k4,minor); \
259}
260
261#define m5create3(s,t) \
262{ \
263 Kinem3 k3=Kinem3(k5st##s##t); \
264 Minor3::Ptr minor=Minor3::create(k3,self,s,t, offs); \
265 MCache::INSERTMINOR3smartinsertMinor3(k3,minor); \
266}
267
268#define m5create2(s,t,u) \
269{ \
270 Kinem2 k2=Kinem2(k5stu##s##t##u); \
271 Minor2::Ptr minor=Minor2::create(k2,self,s,t,u, offs); \
272 MCache::INSERTMINOR2smartinsertMinor2(k2,minor); \
273}
274
275/* --------------------------------------------------------
276 * Real 5-point kinematics
277 * --------------------------------------------------------
278 */
279Minor5::Minor5(const Kinem5& k) : kinem(k), smax(5), pmaxS4(), pmaxS3()
280{
281#ifdef USE_GOLEM_MODE_6
282 psix=6;
283#endif
284 const double p1=kinem.p1();
285 const double p2=kinem.p2();
286 const double p3=kinem.p3();
287 const double p4=kinem.p4();
288 const double p5=kinem.p5();
289 const double s12=kinem.s12();
290 const double s23=kinem.s23();
291 const double s34=kinem.s34();
292 const double s45=kinem.s45();
293 const double s15=kinem.s15();
294 const double m1=kinem.m1();
295 const double m2=kinem.m2();
296 const double m3=kinem.m3();
297 const double m4=kinem.m4();
298 const double m5=kinem.m5();
299
300 Cay[ 0]=2*m1;
301 Cay[ 1]=m1+m2-p2; Cay[ 2]=2*m2;
302 Cay[ 3]=m1+m3-s23; Cay[ 4]=m2+m3-p3; Cay[ 5]=2*m3;
303 Cay[ 6]=m1+m4-s15; Cay[ 7]=m2+m4-s34; Cay[ 8]=m3+m4-p4; Cay[ 9]=2*m4;
304 Cay[10]=m1+m5-p1; Cay[11]=m2+m5-s12; Cay[12]=m3+m5-s45; Cay[13]=m4+m5-p5; Cay[14]=2*m5;
305
306 // create subkinematics minors
307 Ptr self=Ptr(this);
308 const int offs=0;
309
310 m5create4(1);
311 m5create4(2);
312 m5create4(3);
313 m5create4(4);
314 m5create4(5);
315
316 m5create3(1,2);
317 m5create3(1,3);
318 m5create3(1,4);
319 m5create3(1,5);
320 m5create3(2,3);
321 m5create3(2,4);
322 m5create3(2,5);
323 m5create3(3,4);
324 m5create3(3,5);
325 m5create3(4,5);
326
327 m5create2(1,2,3);
328 m5create2(1,2,4);
329 m5create2(1,2,5);
330 m5create2(1,3,4);
331 m5create2(1,3,5);
332 m5create2(1,4,5);
333 m5create2(2,3,4);
334 m5create2(2,3,5);
335 m5create2(2,4,5);
336 m5create2(3,4,5);
337
338 maxCay(); // triggers chain of evalM1, evalM2 and evalM3
339}
340
341/* --------------------------------------------------------
342 * Dummy 5-from-4 kinematics
343 * --------------------------------------------------------
344 */
345Minor5::Minor5(const Kinem4& k) : smax(1), pmaxS4(), pmaxS3()
346{
347#ifdef USE_GOLEM_MODE_6
348 psix=6;
349#endif
350// 12 pinched dummy 5-point kinematics
351 const double p3=k.p2();
352 const double p4=k.p3();
353 const double p5=k.p4();
354 const double s12=k.p1();
355 const double s34=k.s23();
356 const double s45=k.s12();
357 const double m2=k.m1();
358 const double m3=k.m2();
359 const double m4=k.m3();
360 const double m5=k.m4();
361 kinem=Kinem5(0.,0.,p3,p4,p5,s12,0.,s34,s45,0.,0.,m2,m3,m4,m5);
362
363 Cay[ 0]=0;
364 Cay[ 1]=m2; Cay[ 2]=2*m2;
365 Cay[ 3]=m3; Cay[ 4]=m2+m3-p3; Cay[ 5]=2*m3;
366 Cay[ 6]=m4; Cay[ 7]=m2+m4-s34; Cay[ 8]=m3+m4-p4; Cay[ 9]=2*m4;
367 Cay[10]=m5; Cay[11]=m2+m5-s12; Cay[12]=m3+m5-s45; Cay[13]=m4+m5-p5; Cay[14]=2*m5;
368
369 // create subkinematics minors
370 Ptr self=Ptr(this);
371 const int offs=1;
372
373 m5create4(1);
1
Calling 'Minor4::create'
→
7
←
Returning; memory was released
→
8
←
Calling '~SPtr'
→
374
375 m5create3(1,2);
376 m5create3(1,3);
377 m5create3(1,4);
378 m5create3(1,5);
379
380 m5create2(1,2,3);
381 m5create2(1,2,4);
382 m5create2(1,2,5);
383 m5create2(1,3,4);
384 m5create2(1,3,5);
385 m5create2(1,4,5);
386
387 maxCay(); // triggers chain of evalM1, evalM2 and evalM3
388}
389
390#undef m5create4
391#undef m5create3
392#undef m5create2
393
394/* --------------------------------------------------------
395 *
396 * 5-point signed minors
397 *
398 * --------------------------------------------------------
399 */
400
401void Minor5::maxCay()
402{
403 for (int i=1; i<=DCay-1; i++) {
404 for (int ip=i+1; ip<=DCay-1; ip++) {
405 const double m1=kinem.mass(i);
406 const double m2=kinem.mass(ip);
407 const double maxM = m1>m2 ? m1 : m2;
408 pmaxM2[im2(i,ip)-5] = maxM>meps ? maxM : meps; // NOTE meps depends on mu2 scale
409 }
410 }
411
412 for (int i=1; i<=DCay-1; i++) {
413 for (int j=i; j<=DCay-1; j++) {
414 const double cay=fabs(Cay[nss(i,j)]);
415 for (int s=1; s<=smax; s++) {
416 if (s==i || s==j) continue;
417 if (pmaxS4[s-1] < cay) pmaxS4[s-1]=cay;
418 for (int t=s+1; t<=DCay-1; t++) {
419 if (t==i || t==j) continue;
420 const int idx = im2(s,t)-5;
421 if (pmaxS3[idx] < cay ) pmaxS3[idx]=cay;
422 }
423 }
424 }
425 }
426 if (not fEval[E_M1] ) {
427 evalM1();
428 }
429 // Normalize with |G|/|S|
430 for (int s=1; s<=smax; s++) {
431 pmaxS4[s-1]=fabs(pmaxS4[s-1]*M1(s,s)/M2(0,s,0,s));
432 for (int t=s+1; t<=DCay-1; t++) {
433 const int idx=im2(s,t)-5;
434
435 int i=0;
436 double dsits0t=0;
437 for (int ii=1; ii<=DCay-1; ii++) {
438 if (i==s || i==t) continue;
439 const double val=fabs(M3(0,s,t,ii,s,t));
440 if (val > dsits0t) {
441 dsits0t=val;
442 i=ii;
443 }
444 }
445 imax3[idx]=i;
446
447 const double maxcay3=pmaxS3[idx];
448 const double dstst=M2(s,t,s,t);
449 const double ds0ts0t=M3(0,s,t,0,s,t);
450
451 pmaxS3[idx]=fabs((maxcay3*dstst)/ds0ts0t);
452 pmaxT3[idx]=fabs(ds0ts0t/(maxcay3*M3(0,s,t,i,s,t)));
453 pmaxU3[idx]=fabs(dstst/M3(0,s,t,i,s,t));
454 }
455 }
456}
457
458/* --------------------------------------------------------
459 Return one-index minor with proper sign
460 * --------------------------------------------------------
461 */
462double Minor5::M1(int i, int l)
463{
464 return pM1[is(i,l)];
465}
466
467/* --------------------------------------------------------
468 Return two-index minor with proper sign
469 * --------------------------------------------------------
470 */
471double Minor5::M2(int i, int j, int l, int m)
472{
473 int sign=signM2ud(i,j,l,m);
474 if (sign==0) return 0;
475
476 int uidx=im2(i,j);
477 int lidx=im2(l,m);
478
479 return pM2[is(uidx,lidx)]*sign;
480}
481
482/* --------------------------------------------------------
483 Return three-index minor with proper sign
484* --------------------------------------------------------
485*/
486double Minor5::M3(int i, int j, int k, int l, int m, int n)
487{
488 int sign=signM3ud(i,j,k,l,m,n);
489 if (sign==0) return 0;
490
491 int uidx=im3(i,j,k);
492 int lidx=im3(l,m,n);
493
494 return pM3[is(uidx,lidx)]*sign;
495}
496
497/* --------------------------------------------------------
498 Evaluate all 15 one-index minors (need 2-idx minors)
499* --------------------------------------------------------
500*/
501void Minor5::evalM1()
502{
503 if (not fEval[E_M2] ) {
504 evalM2();
505 }
506#ifndef NDEBUG
507 for (int i=0; i<=DCay-1; i++) {
508 for (int l=0; l<=i; l++) {
509 pM1[iss(l,i)]=std::numeric_limits<double>::quiet_NaN();
510 }
511 }
512#endif
513// for (int i=0; i<=0; i++)
514 {
515 const int i=0;
516// for (int l=0; l<=i; l++) {
517 {
518 const int l=0;
519// int j = i==0 ? 1 : 0;
520 const int j=1;
521
522 double m1ele=0;
523 for (int m=1; m<=DCay-1; m++) {
524 m1ele+=M2(i,j,l,m)*Cay[nss(j,m)];
525 }
526 pM1[is(i,l)]=m1ele;
527 }
528 }
529 const int j=0;
530 for (int i=1; i<=smax; i++) {
531 for (int l=0; l<=i; l++) {
532 double m1ele=0;
533 for (int m=1; m<=DCay-1; m++) {
534 m1ele+=M2(i,j,l,m);
535 }
536 pM1[iss(l,i)]=m1ele;
537 }
538 }
539 fEval[E_M1]=true;
540}
541
542/* --------------------------------------------------------
543 Evaluate Gram3 with the least precision loss
544* --------------------------------------------------------
545*/
546double Minor5::gram3(double p1, double p2, double p3)
547{
548 double g3;
549 if (fabs(p1) > fabs(p2)) {
550 if (fabs(p1) > fabs(p3)) {
551 const double diff=(p1 - p2 - p3);
552 const double subs=(-4.)*p2*p3;
553 g3=diff*diff+subs;
554 }
555 else {
556 const double diff=(p3 - p2 - p1);
557 const double subs=(-4.)*p2*p1;
558 g3=diff*diff+subs;
559 }
560 }
561 else {
562 if (fabs(p2) > fabs(p3)) {
563 const double diff=(p2 - p1 - p3);
564 const double subs=(-4.)*p1*p3;
565 g3=diff*diff+subs;
566 }
567 else {
568 const double diff=(p3 - p2 - p1);
569 const double subs=(-4.)*p2*p1;
570 g3=diff*diff+subs;
571 }
572 }
573 return g3;
574}
575
576/* --------------------------------------------------------
577 Evaluate all 120 two-index minors (need 3-idx minors)
578 * --------------------------------------------------------
579 */
580void Minor5::evalM2()
581{
582 if (not fEval[E_M3] ) {
583 evalM3();
584 }
585#ifndef NDEBUG
586 for (int i=0; i<=DCay-1; i++) {
587 for (int j=i+1; j<=DCay-1; j++) {
588 const int uidx=im2(i,j);
589 for (int l=0; l<=DCay-1; l++) {
590 for (int m=l+1; m<=DCay-1; m++) {
591 int lidx=im2(l,m);
592 if (lidx > uidx) continue;
593 pM2[is(uidx,lidx)]=std::numeric_limits<double>::quiet_NaN();
594 }
595 }
596 }
597 }
598#endif
599// for (int i=0; i<=0; i++)
600 {
601 const int i=0;
602 for (int j=i+1; j<=DCay-1; j++) {
603 const int uidx=im2(i,j);
604 const int k = (j==1 ? 2 : 1);
605 for (int l=0; l<=smax; l++) {
606 for (int m=l+1; m<=DCay-1; m++) {
607 int lidx=im2(l,m);
608 if (lidx > uidx) continue;
609
610 double m2ele=M3(i,j,k,l,m,0);
611 for (int n=1; n<DCay; n++) {
612 m2ele+=M3(i,j,k,l,m,n)*Cay[ns(k,n)];
613 }
614 pM2[is(uidx,lidx)]=m2ele;
615 }
616 }
617 }
618 }
619 const int k=0;
620 for (int i=1; i<=smax; i++) {
621 for (int j=i+1; j<=DCay-1; j++) {
622 const int uidx=im2(i,j);
623 for (int l=0; l<=smax; l++) {
624 for (int m=l+1; m<=DCay-1; m++) {
625 int lidx=im2(l,m);
626 if (lidx > uidx) continue;
627
628 double m2ele=0;
629 for (int n=1; n<DCay; n++) {
630 m2ele+=M3(i,j,k,l,m,n);
631 }
632 pM2[is(uidx,lidx)]=m2ele;
633 }
634 }
635 }
636 }
637 fEval[E_M2]=true;
638}
639
640/* --------------------------------------------------------
641 Evaluate all 210 three-index minors
642 * --------------------------------------------------------
643 */
644void Minor5::evalM3()
645{
646#ifndef NDEBUG
647 for (int i=0; i<=DCay-1; i++) {
648 for (int j=i+1; j<=DCay-1; j++) {
649 for (int k=j+1; k<=DCay-1; k++) {
650 const int uidx=im3(i,j,k);
651 for (int l=0; l<=DCay-1; l++) {
652 for (int m=l+1; m<=DCay-1; m++) {
653 for (int n=m+1; n<=DCay-1; n++) {
654 int lidx=im3(l,m,n);
655 if (lidx > uidx) continue;
656 pM3[is(uidx,lidx)]=std::numeric_limits<double>::quiet_NaN();
657 }
658 }
659 }
660 }
661 }
662 }
663#endif
664// for (int i=0; i<=0; i++) {
665 {
666 const int i=0;
667 for (int j=i+1; j<=DCay-2; j++) {
668 for (int k=j+1; k<=DCay-1; k++) {
669 const int uidx=im3(i,j,k);
670// for (int l=0; l<=0; l++) {
671 {
672 const int l=0;
673 for (int m=l+1; m<=DCay-2; m++) {
674 for (int n=m+1; n<=DCay-1; n++) {
675 int lidx=im3(l,m,n);
676 if (lidx > uidx) continue;
677
678 int iu[3]={i,j,k};
679 int nu[3];
680 freeidxM3(iu, nu);
681
682 int id[3]={l,m,n};
683 int nd[3];
684 freeidxM3(id, nd);
685
686 int powsign=-2*((i+j+k+l+m+n)&0x1)+1;
687
688 // nu[0]!=0 and nd[0]!=0
689 pM3[is(uidx,lidx)]=powsign*(
690 + (+Kay(nu[0],nd[1])*Kay(nu[1],nd[2])
691 -Kay(nu[0],nd[2])*Kay(nu[1],nd[1]))*Kay(nu[2],nd[0])
692 + (+Kay(nu[0],nd[2])*Kay(nu[1],nd[0])
693 -Kay(nu[0],nd[0])*Kay(nu[1],nd[2]))*Kay(nu[2],nd[1])
694 + (+Kay(nu[0],nd[0])*Kay(nu[1],nd[1])
695 -Kay(nu[0],nd[1])*Kay(nu[1],nd[0]))*Kay(nu[2],nd[2])
696 );
697 }
698 }
699 }
700 for (int l=1; l<=smax; l++) {
701 for (int m=l+1; m<=DCay-2; m++) {
702 for (int n=m+1; n<=DCay-1; n++) {
703 int lidx=im3(l,m,n);
704 if (lidx > uidx) continue;
705
706 int iu[3]={i,j,k};
707 int nu[3];
708 freeidxM3(iu, nu);
709
710 int id[3]={l,m,n};
711 int nd[3];
712 freeidxM3(id, nd);
713
714 int powsign=-2*((i+j+k+l+m+n)&0x1)+1;
715
716 // nu[0]!=0 and nd[0]==0
717 pM3[is(uidx,lidx)]=powsign*(
718 + (+Kay(nu[0],nd[1])*Kay(nu[1],nd[2])
719 -Kay(nu[0],nd[2])*Kay(nu[1],nd[1]))
720 + (+Kay(nu[0],nd[2])
721 -Kay(nu[1],nd[2]))*Kay(nu[2],nd[1])
722 + (+Kay(nu[1],nd[1])
723 -Kay(nu[0],nd[1]))*Kay(nu[2],nd[2])
724 );
725 }
726 }
727 }
728 }
729 }
730 }
731
732 for (int i=1; i<=smax; i++) {
733 for (int j=i+1; j<=DCay-2; j++) {
734 for (int k=j+1; k<=DCay-1; k++) {
735 const int uidx=im3(i,j,k);
736// for (int l=0; l<=0; l++) {
737 {
738 const int l=0;
739 for (int m=l+1; m<=DCay-2; m++) {
740 for (int n=m+1; n<=DCay-1; n++) {
741 int lidx=im3(l,m,n);
742 if (lidx > uidx) continue;
743
744 int iu[3]={i,j,k};
745 int nu[3];
746 freeidxM3(iu, nu);
747
748 int id[3]={l,m,n};
749 int nd[3];
750 freeidxM3(id, nd);
751
752 int powsign=-2*((i+j+k+l+m+n)&0x1)+1;
753
754 // nu[0]==0 and nd[0]!=0
755 pM3[is(uidx,lidx)]=powsign*(
756 + (+Kay(nu[1],nd[2])
757 -Kay(nu[1],nd[1]))*Kay(nu[2],nd[0])
758 + (+Kay(nu[1],nd[0])
759 -Kay(nu[1],nd[2]))*Kay(nu[2],nd[1])
760 + (+Kay(nu[1],nd[1])
761 -Kay(nu[1],nd[0]))*Kay(nu[2],nd[2])
762 );
763 }
764 }
765 }
766 for (int l=1; l<=smax; l++) {
767 for (int m=l+1; m<=DCay-2; m++) {
768 for (int n=m+1; n<=DCay-1; n++) {
769 int lidx=im3(l,m,n);
770 if (lidx > uidx) continue;
771
772 int iu[3]={i,j,k};
773 int nu[3];
774 freeidxM3(iu, nu);
775
776 int id[3]={l,m,n};
777 int nd[3];
778 freeidxM3(id, nd);
779
780 int powsign=-2*((i+j+k+l+m+n)&0x1)+1;
781
782 // nu[0]==0 and nd[0]==0
783 pM3[is(uidx,lidx)]=powsign*(
784 + Kay(nu[1],nd[2])
785 - Kay(nu[1],nd[1])
786 + Kay(nu[2],nd[1])
787 - Kay(nu[2],nd[2])
788 );
789 }
790 }
791 }
792 }
793 }
794 }
795 fEval[E_M3]=true;
796}
797
798/* --------------------------------------------------------
799 *
800 * Plain scalar integrals
801 *
802 * --------------------------------------------------------
803 */
804
805/* --------------------------------------------------------
806 I4s box
807 * --------------------------------------------------------
808 */
809ncomplex Minor5::I4s(int ep, int s)
810{
811 if (not fEval[E_I4s+ep]) {
812 I4sEval(ep);
813 }
814 return pI4s[ep][s-1];
815}
816void Minor5::I4sEval(int ep) // IR-div
817{
818 // Kinematics is in LT notation (for calling qcdloop)
819 double p1=kinem.p1();
820 double p2=kinem.p2();
821 double p3=kinem.p3();
822 double p4=kinem.p4();
823 double p5=kinem.p5();
824 double s12=kinem.s12();
825 double s23=kinem.s23();
826 double s34=kinem.s34();
827 double s45=kinem.s45();
828 double s15=kinem.s15();
829 double m1=kinem.m1();
830 double m2=kinem.m2();
831 double m3=kinem.m3();
832 double m4=kinem.m4();
833 double m5=kinem.m5();
834
835 pI4s[ep][1-1]=ICache::getI4(ep, Kinem4(s12,p3,p4,p5,s45,s34,m5,m2,m3,m4));
836 if (smax==5) {
837 pI4s[ep][2-1]=ICache::getI4(ep, Kinem4(p1,s23,p4,p5,s45,s15,m5,m1,m3,m4));
838 pI4s[ep][3-1]=ICache::getI4(ep, Kinem4(p1,p2,s34,p5,s12,s15,m5,m1,m2,m4));
839 pI4s[ep][4-1]=ICache::getI4(ep, Kinem4(p1,p2,p3,s45,s12,s23,m5,m1,m2,m3));
840 pI4s[ep][5-1]=ICache::getI4(ep, Kinem4(p2,p3,p4,s15,s23,s34,m1,m2,m3,m4));
841 }
842 fEval[E_I4s+ep]=true;
843}
844
845/* --------------------------------------------------------
846 * I3st triangle
847 * --------------------------------------------------------
848 */
849ncomplex Minor5::I3st(int ep, int s, int t) // IR-div
850{
851 assert(s!=t)(static_cast <bool> (s!=t) ? void (0) : __assert_fail (
"s!=t", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 851, __extension__ __PRETTY_FUNCTION__));
852 if (not fEval[E_I3st+ep]) {
853 I3stEval(ep);
854 }
855 int idx = im2(s,t)-5;
856 return pI3st[ep][idx];
857}
858
859void Minor5::I3stEval(int ep)
860{
861 // Kinematics is in LT notation (for calling qcdloop)
862 double p1=kinem.p1();
863 double p2=kinem.p2();
864 double p3=kinem.p3();
865 double p4=kinem.p4();
866 double p5=kinem.p5();
867 double s12=kinem.s12();
868 double s23=kinem.s23();
869 double s34=kinem.s34();
870 double s45=kinem.s45();
871 double s15=kinem.s15();
872 double m1=kinem.m1();
873 double m2=kinem.m2();
874 double m3=kinem.m3();
875 double m4=kinem.m4();
876 double m5=kinem.m5();
877
878 // it is symmetric with zeroes on main diagonal
879 pI3st[ep][im2(1,2)-5]=ICache::getI3(ep, Kinem3(s45,p4, p5, m5,m3,m4));
880 pI3st[ep][im2(1,3)-5]=ICache::getI3(ep, Kinem3(s12,s34,p5, m5,m2,m4));
881 pI3st[ep][im2(1,4)-5]=ICache::getI3(ep, Kinem3(s12,p3, s45,m5,m2,m3));
882 pI3st[ep][im2(1,5)-5]=ICache::getI3(ep, Kinem3(p3, p4, s34,m2,m3,m4));
883 if (smax==5) {
884 pI3st[ep][im2(2,3)-5]=ICache::getI3(ep, Kinem3(p1, s15,p5, m5,m1,m4));
885 pI3st[ep][im2(2,4)-5]=ICache::getI3(ep, Kinem3(p1, s23,s45,m5,m1,m3));
886 pI3st[ep][im2(2,5)-5]=ICache::getI3(ep, Kinem3(s23,p4, s15,m1,m3,m4));
887 pI3st[ep][im2(3,4)-5]=ICache::getI3(ep, Kinem3(p1, p2, s12,m5,m1,m2));
888 pI3st[ep][im2(3,5)-5]=ICache::getI3(ep, Kinem3(p2, s34,s15,m1,m2,m4));
889 pI3st[ep][im2(4,5)-5]=ICache::getI3(ep, Kinem3(p2, p3, s23,m1,m2,m3));
890 }
891 fEval[E_I3st+ep]=true;
892}
893
894/* --------------------------------------------------------
895 * I2stu bubble
896 * --------------------------------------------------------
897 */
898ncomplex Minor5::I2stu(int ep, int s, int t, int u)
899{
900 assert(t!=u && u!=s && s!=t)(static_cast <bool> (t!=u && u!=s && s!=
t) ? void (0) : __assert_fail ("t!=u && u!=s && s!=t"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 900
, __extension__ __PRETTY_FUNCTION__));
901 if (ep>=2) return 0;
902
903 if (not fEval[E_I2stu+ep]) {
904 I2stuEval(ep);
905 }
906 int idx=im3(s,t,u)-10;
907 return pI2stu[ep][idx];
908}
909
910void Minor5::I2stuEval(int ep)
911{
912 // Kinematics is in LT notation (for calling qcdloop)
913 double p1=kinem.p1();
914 double p2=kinem.p2();
915 double p3=kinem.p3();
916 double p4=kinem.p4();
917 double p5=kinem.p5();
918 double s12=kinem.s12();
919 double s23=kinem.s23();
920 double s34=kinem.s34();
921 double s45=kinem.s45();
922 double s15=kinem.s15();
923 double m1=kinem.m1();
924 double m2=kinem.m2();
925 double m3=kinem.m3();
926 double m4=kinem.m4();
927 double m5=kinem.m5();
928
929 // it is symmetric with zeroes on main diagonal
930 pI2stu[ep][im3(1,2,3)-10]=ICache::getI2(ep, Kinem2( p5, m5, m4));
931 pI2stu[ep][im3(1,2,4)-10]=ICache::getI2(ep, Kinem2(s45, m5, m3));
932 pI2stu[ep][im3(1,2,5)-10]=ICache::getI2(ep, Kinem2( p4, m3, m4));
933 pI2stu[ep][im3(1,3,4)-10]=ICache::getI2(ep, Kinem2(s12, m5, m2));
934 pI2stu[ep][im3(1,3,5)-10]=ICache::getI2(ep, Kinem2(s34, m2, m4));
935 pI2stu[ep][im3(1,4,5)-10]=ICache::getI2(ep, Kinem2( p3, m2, m3));
936 if (smax==5) {
937 pI2stu[ep][im3(2,3,4)-10]=ICache::getI2(ep, Kinem2( p1, m5, m1));
938 pI2stu[ep][im3(2,3,5)-10]=ICache::getI2(ep, Kinem2(s15, m1, m4));
939 pI2stu[ep][im3(2,4,5)-10]=ICache::getI2(ep, Kinem2(s23, m1, m3));
940 pI2stu[ep][im3(3,4,5)-10]=ICache::getI2(ep, Kinem2( p2, m1, m2));
941 }
942 fEval[E_I2stu+ep]=true;
943}
944
945/* --------------------------------------------------------
946 *
947 * Rank-1 functions
948 *
949 * --------------------------------------------------------
950 */
951
952/* --------------------------------------------------------
953 * I4Ds box in D+2 dim
954 * --------------------------------------------------------
955 */
956ncomplex Minor5::I4Ds(int ep, int s)
957{
958 if (ep!=0) return 0; // I4Ds is finite
959 if (not fEval[E_I4Ds+ep]) {
960 I4DsEval(ep);
961 }
962 return pI4Ds[ep][s-1];
963}
964
965void Minor5::I4DsEval(const int ep)
966{
967 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 967, __extension__ __PRETTY_FUNCTION__));
968 for (int s=1; s<=smax; s++) {
969 const double dss=M1(s, s);
970 const double d0s0s=M2(0, s, 0, s);
971 ncomplex ivalue=0;
972
973 ncomplex sum1=0;
974 for (int t=1; t<=5; t++) {
975 if (t==s) continue;
976 sum1-=M2(s,t,s,0)*I3st(ep, s, t);
977 }
978 sum1+=d0s0s*I4s(ep, s);
979 ivalue=sum1/dss;
980
981 const double x=dss/d0s0s;
982 if (pmaxS4[s-1] <= deps1) {
983 ncomplex sump;
984 do {
985 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 985, __extension__ __PRETTY_FUNCTION__));
986
987 double dsts0[6];
988 sum1=0;
989 for (int t=1; t<=5; t++) {
990 if (t==s) continue;
991 dsts0[t]=M2(s,t,s,0);
992 sum1+=dsts0[t]*I3Dst(0, s, t);
993 }
994
995 double xn=1;
996 ncomplex dv,s21;
997
998 ncomplex sum[3];
999 sum[0]=sump=sum1;
1000
1001#define stepI4D(n,a,b) \
1002 xn*=x; \
1003 dv=0; \
1004 for (int t=1; t<=5; t++) { \
1005 if (t==s) continue; \
1006 dv+=dsts0[t]*(a*I3D##n##st(0, s, t) - b*I3D##n##st(1, s, t)); \
1007 } \
1008 dv*=xn; \
1009 sum1+=dv;
1010
1011 stepI4D(2,3.,2.)
1012 if ( fabs(sum1.real()*teps)>=fabs(dv.real())
1013 && fabs(sum1.imag()*teps)>=fabs(dv.imag()))
1014 break;
1015 sum[1]=sump=sum1;
1016 s21=sum[1]-sum[0];
1017
1018 stepI4D(3,15.,16.)
1019 sump=sum1;
1020 stepWynn(0)sum[(2+0)%3]=sum1; { const ncomplex s2=sum[(2+0)%3]; const ncomplex
s1=sum[(1+0)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1021 stepI4D(4,105.,142.)
1022 stepWynn(1)sum[(2+1)%3]=sum1; { const ncomplex s2=sum[(2+1)%3]; const ncomplex
s1=sum[(1+1)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1023 stepI4D(5,945.,1488.)
1024 stepWynn(2)sum[(2+2)%3]=sum1; { const ncomplex s2=sum[(2+2)%3]; const ncomplex
s1=sum[(1+2)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1025 stepI4D(6,10395.,18258.)
1026 stepWynn(3)sum[(2+3)%3]=sum1; { const ncomplex s2=sum[(2+3)%3]; const ncomplex
s1=sum[(1+3)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1027 stepI4D(7,135135.,258144.)
1028 stepWynn(4)sum[(2+4)%3]=sum1; { const ncomplex s2=sum[(2+4)%3]; const ncomplex
s1=sum[(1+4)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1029// stepI4D(8,2027025.,4142430.)
1030// stepWynn(5)
1031// stepI4D(9,34459425.,74475360.)
1032// stepWynn(6)
1033#undef stepI4D
1034 } while (0);
1035 ivalue=sump/d0s0s;
1036 }
1037 pI4Ds[ep][s-1]=ivalue;
1038 }
1039 fEval[E_I4Ds+ep]=true;
1040}
1041
1042/* --------------------------------------------------------
1043 * I4Dsi box in D+2 dim with a dot
1044 * --------------------------------------------------------
1045 */
1046ncomplex Minor5::I4Dsi(int ep, int s, int i) // IR-div
1047{
1048 if (s==i) return 0;
1049 if (not fEval[E_I4Dsi+ep]) {
1050 I4DsiEval(ep);
1051 }
1052 return pI4Dsi[ep][i-1][s-1];
1053}
1054
1055void Minor5::I4DsiEval(int ep)
1056{
1057 for (int s=1; s<=smax; s++) {
1058 for (int i=1; i<=CIDX(DCay-2); i++) {
1059 if (s==i) continue;
1060 ncomplex ivalue=0;
1061
1062 if (pmaxS4[s-1] <= deps1 || (fabs(M1(s,s))<1e-11 && fabs(M2(0,s,0,s))>0) ) {
1063 ncomplex sum1=0;
1064 for (int t=1; t<=5; t++) {
1065 if (t==s) continue;
1066 sum1+=M3(0, s, t, 0, s, i)*I3st(ep, s, t);
1067 }
1068 sum1-=M2(0, s, i, s)*I4Ds(ep, s);
1069 ivalue=sum1/M2(0, s, 0, s);
1070 } else {
1071 ncomplex sum1=0;
1072 for (int t=1; t<=5; t++) {
1073 if (t==s) continue;
1074 sum1+=M2(s, t, s, i)*I3st(ep, s, t);
1075 }
1076 sum1-=M2(0, s, i, s)*I4s(ep, s);
1077 ivalue=sum1/M1(s, s);
1078 }
1079 pI4Dsi[ep][i-1][s-1]=ivalue;
1080 }
1081 }
1082 fEval[E_I4Dsi+ep]=true;
1083}
1084
1085/* --------------------------------------------------------
1086 *
1087 * Rank-2 functions
1088 *
1089 * --------------------------------------------------------
1090 */
1091
1092/* --------------------------------------------------------
1093 * I3Dst triangle in D+2 dim
1094 * --------------------------------------------------------
1095 */
1096ncomplex Minor5::I3Dst(int ep, int s, int t)
1097{
1098 assert(s!=t)(static_cast <bool> (s!=t) ? void (0) : __assert_fail (
"s!=t", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1098, __extension__ __PRETTY_FUNCTION__));
1099 if (ep==1) return -0.5;
1100 else if (ep==2) return 0;
1101 if (not fEval[E_I3Dst+ep]) {
1102 I3DstEval(ep);
1103 }
1104 int idx = im2(s,t)-5;
1105 return pI3Dst[ep][idx];
1106}
1107
1108void Minor5::I3DstEval(int ep)
1109{
1110 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1110, __extension__ __PRETTY_FUNCTION__));
1111 for (int s=1; s<=smax; s++) {
1112 for (int t=s+1; t<=5; t++) {
1113 int idx = im2(s,t)-5;
1114 const double dstst=M2(s,t,s,t);
1115 const double d0st0st=M3(0,s,t,0,s,t);
1116 ncomplex ivalue=0;
1117
1118 if ( pmaxT3[idx]!=0 && ( pmaxT3[idx] <= epsir1 && pmaxU3[idx] <= epsir1 ) ) {
1119 // IR
1120 int i=imax3[idx];
1121 int iu[3]={i-1,s-1,t-1};
1122 int tmp;
1123 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
1124 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
1125 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
1126 int nu[3];
1127 freeidxM3(iu, nu);
1128 int u=nu[0]+1;
1129 int v=nu[1]+1;
1130 const double Dii=M4ii(u,v,i);
1131 const double Dui=M4ui(u,v,i);
1132 const double Dvi=M4vi(u,v,i);
1133 ncomplex sum1=+Dii*(I2stu(0, s, t, i)+I2stu(1, s, t, i)) // (i, i)
1134 +Dui*(I2stu(0, s, t, u)+I2stu(1, s, t, u)) // (u, i)
1135 +Dvi*(I2stu(0, s, t, v)+I2stu(1, s, t, v)); // (v, i)
1136 ivalue=0.5*sum1/M3(0,s,t,i,s,t);
1137 } else if (pmaxS3[idx] <= ceps) {
1138 // EXP
1139 const double x=dstst/d0st0st;
1140 ncomplex sump;
1141 do {
1142 double dstust0[6];
1143 ncomplex sum1=0;
1144 for (int u=1; u<=5; u++) {
1145 if (u==t || u==s) continue;
1146 dstust0[u]=M3(s,t,u,s,t,0);
1147 sum1+=dstust0[u]*I2Dstu(0, s, t, u);
1148 }
1149
1150 double xn=1;
1151 ncomplex dv,s21;
1152
1153 ncomplex sum[3];
1154 sum[0]=sump=sum1;
1155
1156#define stepI3D(n,a,b) \
1157 xn*=x; \
1158 dv=0; \
1159 for (int u=1; u<=5; u++) { \
1160 if (u==t || u==s) continue; \
1161 dv+=dstust0[u]*(a*I2D##n##stu(0, s, t, u) - b*I2D##n##stu(1, s, t, u)); \
1162 } \
1163 dv*=xn; \
1164 sum1+=dv;
1165
1166 stepI3D(2,4.,2.)
1167 if ( fabs(sum1.real()*teps)>=fabs(dv.real())
1168 && fabs(sum1.imag()*teps)>=fabs(dv.imag()))
1169 break;
1170 sum[1]=sump=sum1;
1171 s21=sum[1]-sum[0];
1172
1173 stepI3D(3,24.,20.)
1174 sump=sum1;
1175 stepWynn(0)sum[(2+0)%3]=sum1; { const ncomplex s2=sum[(2+0)%3]; const ncomplex
s1=sum[(1+0)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1176 stepI3D(4,192.,208.)
1177 stepWynn(1)sum[(2+1)%3]=sum1; { const ncomplex s2=sum[(2+1)%3]; const ncomplex
s1=sum[(1+1)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1178 stepI3D(5,1920.,2464.)
1179 stepWynn(2)sum[(2+2)%3]=sum1; { const ncomplex s2=sum[(2+2)%3]; const ncomplex
s1=sum[(1+2)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1180 stepI3D(6,23040.,33408.)
1181 stepWynn(3)sum[(2+3)%3]=sum1; { const ncomplex s2=sum[(2+3)%3]; const ncomplex
s1=sum[(1+3)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1182// stepI3D(7,322560.,513792.)
1183// stepWynn(4)
1184#undef stepI3D
1185 } while (0);
1186 ivalue=sump/d0st0st;
1187 } else {
1188 // NORMAL
1189 ncomplex sum1=0;
1190 for (int u=1; u<=5; u++) {
1191 if (u==t || u==s) continue;
1192 sum1-=M3(u,s,t,0,s,t)*I2stu(ep, s, t, u);
1193 }
1194 sum1+=d0st0st*I3st(ep, s, t);
1195 ivalue=sum1/(2*dstst)-0.5; // 2*(-1/2)/2 == -0.5
1196 }
1197 pI3Dst[ep][idx]=ivalue;
1198 }
1199 }
1200 fEval[E_I3Dst+ep]=true;
1201}
1202
1203
1204/* --------------------------------------------------------
1205 * I4D2s box in D+4 dim
1206 * --------------------------------------------------------
1207 */
1208ncomplex Minor5::I4D2s(int ep, int s)
1209{
1210 if (ep==1) return 1./6.;
1211 else if (ep==2) return 0;
1212 if (not fEval[E_I4D2s+ep]) {
1213 I4D2sEval(ep);
1214 }
1215 return pI4D2s[ep][s-1];
1216}
1217
1218void Minor5::I4D2sEval(int ep) {
1219 for (int s=1; s<=smax; s++) {
1220 const double dss=M1(s, s);
1221 const double d0s0s=M2(0, s, 0, s);
1222 ncomplex ivalue=0;
1223
1224 ncomplex sum1=0;
1225 for (int t=1; t<=5; t++) {
1226 if (t==s) continue;
1227 sum1-=M2(s,t,s,0)*I3Dst(ep, s, t);
1228 }
1229 sum1+=d0s0s*I4Ds(ep, s);
1230 ivalue=sum1/(3*dss)+1./9.; // +2*(1/6)/3
1231
1232 const double x=dss/d0s0s;
1233 if (pmaxS4[s-1] <= deps2) {
1234 ncomplex sump;
1235 do {
1236 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1236, __extension__ __PRETTY_FUNCTION__));
1237
1238 double dsts0[6];
1239 sum1=0;
1240 for (int t=1; t<=5; t++) {
1241 if (t==s) continue;
1242 dsts0[t]=M2(s,t,s,0);
1243 sum1+=dsts0[t]*I3D2st(0, s, t);
1244 }
1245
1246 double xn=1;
1247 ncomplex dv,s21;
1248
1249 ncomplex sum[3];
1250 sum[0]=sump=sum1;
1251
1252#define stepI4D(n,a,b) \
1253 xn*=x; \
1254 dv=0; \
1255 for (int t=1; t<=5; t++) { \
1256 if (t==s) continue; \
1257 dv+=dsts0[t]*(a*I3D##n##st(0, s, t) - b*I3D##n##st(1, s, t)); \
1258 } \
1259 dv*=xn; \
1260 sum1+=dv;
1261
1262 stepI4D(3,5.,2.)
1263 if ( fabs(sum1.real()*teps)>=fabs(dv.real())
1264 && fabs(sum1.imag()*teps)>=fabs(dv.imag()))
1265 break;
1266 sum[1]=sump=sum1;
1267 s21=sum[1]-sum[0];
1268
1269 stepI4D(4,35.,24.)
1270 sump=sum1;
1271 stepWynn(0)sum[(2+0)%3]=sum1; { const ncomplex s2=sum[(2+0)%3]; const ncomplex
s1=sum[(1+0)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1272 stepI4D(5,315.,286.)
1273 stepWynn(1)sum[(2+1)%3]=sum1; { const ncomplex s2=sum[(2+1)%3]; const ncomplex
s1=sum[(1+1)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1274 stepI4D(6,3465.,3776.)
1275 stepWynn(2)sum[(2+2)%3]=sum1; { const ncomplex s2=sum[(2+2)%3]; const ncomplex
s1=sum[(1+2)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1276 stepI4D(7,45045.,56018.)
1277 stepWynn(3)sum[(2+3)%3]=sum1; { const ncomplex s2=sum[(2+3)%3]; const ncomplex
s1=sum[(1+3)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1278// stepI4D(8,675675.,930360.)
1279// stepWynn(4)
1280// stepI4D(9,11486475.,17167470.)
1281// stepWynn(5)
1282#undef stepI4D
1283 } while (0);
1284 ivalue=sump/d0s0s;
1285 }
1286 pI4D2s[ep][s-1]=ivalue;
1287 }
1288 fEval[E_I4D2s+ep]=true;
1289}
1290
1291/* --------------------------------------------------------
1292 *
1293 * Rank-3 functions
1294 *
1295 * --------------------------------------------------------
1296 */
1297
1298/* --------------------------------------------------------
1299 * I4D2si box in D+4 dim with a dot
1300 * --------------------------------------------------------
1301 */
1302ncomplex Minor5::I4D2si(int ep, int s, int i)
1303{
1304 if (s==i) return 0;
1305 if (ep!=0) return 0; // I4D2si is finite
1306 if (not fEval[E_I4D2si+ep]) {
1307 I4D2siEval(ep);
1308 }
1309 return pI4D2si[ep][i-1][s-1];
1310}
1311
1312void Minor5::I4D2siEval(int ep)
1313{
1314 for (int s=1; s<=smax; s++) {
1315 for (int i=1; i<=CIDX(DCay-2); i++) {
1316 if (s==i) continue;
1317 ncomplex ivalue=0;
1318
1319 if (pmaxS4[s-1] <= deps2 || (fabs(M1(s,s))<1e-11 && fabs(M2(0,s,0,s))>0) ) {
1320 ncomplex sum1=0;
1321 for (int t=1; t<=5; t++) {
1322 if (t==s) continue;
1323 sum1+=M3(0, s, t, 0, s, i)*I3Dst(ep, s, t);
1324 }
1325 sum1+=M2(0, s, i, s)*(-3.*I4D2s(ep, s)+1./3.); // 1/3 == 2*1/6
1326 ivalue=sum1/M2(0, s, 0, s);
1327 } else {
1328 ncomplex sum1=0;
1329 for (int t=1; t<=5; t++) {
1330 if (t==s) continue;
1331 sum1+=M2(s, t, s, i)*I3Dst(ep, s, t);
1332 }
1333 sum1-=M2(0, s, i, s)*I4Ds(ep, s);
1334 ivalue=sum1/M1(s, s);
1335 }
1336
1337 /* // Test for formula 6.11
1338 const double ds0=M1(s, 0);
1339 ncomplex sum1=M2(0,s,0,i)*I4Ds(ep, s)-3*M1(s,i)*I4D2s(ep,s);
1340 sum1+=(ep == 2) ? 0 : 2*M1(s,i)*I4D2s(ep+1,s);
1341 for (int t=1; t<=5; t++) {
1342 sum1+=-M2(s,t,i,0)*I3Dst(ep, s, t);
1343 }
1344 ncomplex ivalue=sum1/ds0;
1345 */
1346
1347 pI4D2si[ep][i-1][s-1]=ivalue;
1348 }
1349 }
1350 fEval[E_I4D2si+ep]=true;
1351}
1352
1353/* --------------------------------------------------------
1354 * I3Dsti triangle in D+2 dim with a dot
1355 * --------------------------------------------------------
1356 */
1357ncomplex Minor5::I3Dsti(int ep, int s, int t, int i) // IR-div
1358{
1359 assert(s!=t && s!=i && t!=i)(static_cast <bool> (s!=t && s!=i && t!=
i) ? void (0) : __assert_fail ("s!=t && s!=i && t!=i"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 1359
, __extension__ __PRETTY_FUNCTION__));
1360 if (not fEval[E_I3Dsti+ep]) {
1361 I3DstiEval(ep);
1362 }
1363 int idx = im2(s,t)-5;
1364 return pI3Dsti[ep][i-1][idx];
1365}
1366
1367void Minor5::I3DstiEval(int ep)
1368{
1369 for (int i=1; i<=CIDX(DCay-2); i++) {
1370 for (int s=1; s<=smax; s++) { if (i==s) continue;
1371 for (int t=s+1; t<=5; t++) { if (i==t) continue;
1372 int idx = im2(s,t)-5;
1373
1374 const double ds0ts0t=M3(s,0,t,s,0,t);
1375 if (ep!=0 && fabs(ds0ts0t) > m3eps) { // if ds0ts0t!=0 I3Dsti is finite
1376 pI3Dsti[ep][i-1][idx]=0;
1377 continue;
1378 }
1379
1380 ncomplex ivalue=0;
1381
1382 if ( ep!=0 ||
1383 ( (pmaxT3[idx]==0 || (pmaxT3[idx] > epsir2 || pmaxU3[idx] > epsir2))
1384 && pmaxS3[idx] > ceps ) ) {
1385 // Variant with Gram3
1386 ncomplex sum1=0;
1387 for (int u=1; u<=5; u++) {
1388 if (u==t || u==s) continue;
1389 sum1+=M3(u,s,t,i,s,t)*I2stu(ep,s,t,u);
1390 }
1391 sum1-=M3(0,s,t,i,s,t)*I3st(ep, s, t);
1392 ivalue=sum1/M2(s,t,s,t);
1393 } else {
1394 ncomplex sum1=0;
1395 int iu[3]={i-1,s-1,t-1};
1396 int tmp;
1397 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
1398 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
1399 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
1400 int nu[3];
1401 freeidxM3(iu, nu);
1402 int u=nu[0]+1;
1403 int v=nu[1]+1;
1404
1405 if ( pmaxT3[idx] <= epsir2 && pmaxU3[idx] <= epsir2 ) {
1406 // small G3 & C3
1407 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1407, __extension__ __PRETTY_FUNCTION__));
1408 int j=imax3[idx];
1409 sum1=0;
1410 ncomplex const I3term=I3st(ep,s,t)+2.*I3st(ep+1,s,t);
1411 ncomplex const I2Uterm=I2stui(ep,s,t,u,i,v)+2.*I2stui(ep+1,s,t,u,i,v);
1412 ncomplex const I2Vterm=I2stui(ep,s,t,v,i,u)+2.*I2stui(ep+1,s,t,v,i,u);
1413 if (j==i) { // j->i
1414 const double Dii=M4ii(u,v,i);
1415 const double Dui=M4ui(u,v,i);
1416 const double Dvi=M4vi(u,v,i);
1417 sum1+=+Dii*(I3term) // (i, i)
1418 +Dui*(I2Uterm) // (u, i)
1419 +Dvi*(I2Vterm); // (v, i)
1420 } else if (j==u) { // j->u
1421 const double Dui=M4ui(u,v,i);
1422 const double Duu=M4uu(u,v,i);
1423 const double Dvu=M4vu(u,v,i);
1424 sum1+=+Dui*(I3term) // (u, i)
1425 +Duu*(I2Uterm) // (u, u)
1426 +Dvu*(I2Vterm); // (v, u)
1427 } else { assert(j==v)(static_cast <bool> (j==v) ? void (0) : __assert_fail (
"j==v", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1427, __extension__ __PRETTY_FUNCTION__)); // j->v
1428 const double Dvi=M4vi(u,v,i);
1429 const double Dvu=M4vu(u,v,i);
1430 const double Dvv=M4vv(u,v,i);
1431 sum1+=+Dvi*(I3term) // (v, i)
1432 +Dvu*(I2Uterm) // (v, u)
1433 +Dvv*(I2Vterm); // (v, v)
1434 }
1435 ivalue=sum1/(M3(s,0,t,s,j,t));
1436 } else {
1437 // small G3
1438 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1438, __extension__ __PRETTY_FUNCTION__));
1439 const double Dii=M4ii(u,v,i);
1440 const double Dui=M4ui(u,v,i);
1441 const double Dvi=M4vi(u,v,i);
1442 sum1+=Dii*I2stu(ep,s,t,i) // (i, i)
1443 +Dui*I2stu(ep,s,t,u) // (u, i)
1444 +Dvi*I2stu(ep,s,t,v); // (v, i)
1445 sum1+=M3(s,0,t,s,i,t)*(-2.*I3Dst(ep, s, t)-1.); //+2.*I3Dst(ep+1, s, t));
1446 ivalue=sum1/ds0ts0t;
1447 }
1448 }
1449 pI3Dsti[ep][i-1][idx]=ivalue;
1450 }
1451 }
1452 }
1453 fEval[E_I3Dsti+ep]=true;
1454}
1455
1456/* --------------------------------------------------------
1457 * I4D2sij box in D+4 dim with two dots
1458 * --------------------------------------------------------
1459 */
1460ncomplex Minor5::I4D2sij(int ep, int s, int i, int j) // IR-div
1461{
1462 if (s==i || s==j) return 0;
1463 if (not fEval[E_I4D2sij+ep]) {
1464 I4D2sijEval(ep);
1465 }
1466 return pI4D2sij[ep][is(i-1,j-1)][s-1];
1467}
1468
1469void Minor5::I4D2sijEval(int ep)
1470{
1471 for (int s=1; s<=smax; s++) {
1472 // symmetric in 'i,j'
1473 for (int i=1; i<=CIDX(DCay-2); i++) { if (s==i) continue;
1474 for (int j=i; j<=CIDX(DCay-2); j++) { if (s==j) continue;
1475 ncomplex ivalue=0;
1476
1477 if (pmaxS4[s-1] <= deps2 || (fabs(M1(s,s))<1e-11 && fabs(M2(0,s,0,s))>0) ) {
1478 ncomplex sum1=0;
1479 for (int t=1; t<=5; t++) {
1480 if (t==s || t==i) continue;
1481 sum1+=M3(s,0,t,s,0,j)*I3Dsti(ep, s, t, i);
1482 }
1483 sum1+=M3(s,0,i,s,0,j)*I4Ds(ep, s);
1484 sum1+=M2(s,0,s,j)*(-2.*I4D2si(ep, s, i));
1485 ivalue=sum1/M2(0, s, 0, s);
1486 } else {
1487 ncomplex sum1=0;
1488 for (int t=1; t<=5; t++) {
1489 if (t==s || t==i) continue;
1490 sum1+=M2(s,t,s,j)*I3Dsti(ep, s, t, i);
1491 }
1492 sum1+=M2(s,i,s,j)*I4Ds(ep, s);
1493 sum1-=M2(s,0,s,j)*I4Dsi(ep, s, i);
1494 ivalue=sum1/M1(s, s);
1495 }
1496 pI4D2sij[ep][iss(i-1,j-1)][s-1]=ivalue;
1497 }
1498 }
1499 }
1500 fEval[E_I4D2sij+ep]=true;
1501}
1502
1503/* --------------------------------------------------------
1504 *
1505 * Rank-4 functions
1506 *
1507 * --------------------------------------------------------
1508 */
1509
1510/* --------------------------------------------------------
1511 * I2Dstu bubble in D+2 dim
1512 * --------------------------------------------------------
1513 */
1514ncomplex Minor5::I2Dstu(int ep, int s, int t, int u)
1515{
1516 assert(t!=u && u!=s && s!=t)(static_cast <bool> (t!=u && u!=s && s!=
t) ? void (0) : __assert_fail ("t!=u && u!=s && s!=t"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 1516
, __extension__ __PRETTY_FUNCTION__));
1517 if (ep==2) return 0;
1518 if (not fEval[E_I2Dstu+ep]) {
1519 I2DstuEval(0,ep,1,2,3,4,5,kinem.p5());
1520 I2DstuEval(1,ep,1,2,4,3,5,kinem.s45());
1521 I2DstuEval(2,ep,1,2,5,3,4,kinem.p4());
1522
1523 I2DstuEval(3,ep,1,3,4,2,5,kinem.s12());
1524 I2DstuEval(4,ep,1,3,5,2,4,kinem.s34());
1525
1526 I2DstuEval(5,ep,1,4,5,2,3,kinem.p3());
1527
1528 if (smax==5) {
1529 I2DstuEval(6,ep,2,3,4,1,5,kinem.p1());
1530 I2DstuEval(7,ep,2,3,5,1,4,kinem.s15());
1531
1532 I2DstuEval(8,ep,2,4,5,1,3,kinem.s23());
1533
1534
1535 I2DstuEval(9,ep,3,4,5,1,2,kinem.p2());
1536 }
1537
1538 fEval[E_I2Dstu+ep]=true;
1539 }
1540 int idx=im3(s,t,u)-10;
1541 return pI2Dstu[ep][idx];
1542}
1543
1544void Minor5::I2DstuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq)
1545{
1546 ncomplex sum1=0;
1547 if (ep==0) {
1548 const double dstustu=-2*qsq; /*M3(s,t,u,s,t,u);*/
1549
1550 const double msq1=kinem.mass(m);
1551 const double msq2=kinem.mass(n);
1552 const double s_cutoff=seps1*pmaxM2[im2(m,n)-5];
1553
1554 if (fabs(dstustu) <= s_cutoff) {
1555 const double mm12=msq1-msq2;
1556 if (fabs(mm12) < meps) {
1557 sum1=-ICache::getI1(ep, Kinem1(msq1));
1558 }
1559 else {
1560 sum1= - 0.25*( msq1 + msq2 )
1561 + 0.5*( - msq1*ICache::getI1(ep, Kinem1(msq1))
1562 + msq2*ICache::getI1(ep, Kinem1(msq2))
1563 )/(mm12);
1564 }
1565 }
1566 else {
1567 ncomplex sumX=3.*I2stu(ep,s,t,u)+2.*I2stu(ep+1,s,t,u);
1568 ncomplex sumY=3.*ICache::getI1(ep, Kinem1(msq2))+2*msq2;
1569 ncomplex sumZ=3.*ICache::getI1(ep, Kinem1(msq1))+2*msq1;
1570
1571 const double ds0tu=(Cay[nss(m,m)]*Cay[nss(n,n)]-Cay[nss(m,n)]*Cay[nss(m,n)]);
1572 sum1+=sumX*ds0tu;
1573
1574 const double dsvtuY=-(Cay[nss(n,n)]-Cay[nss(m,n)]); /* minus sign of minor v=m */
1575 sum1-=sumY*dsvtuY;
1576
1577 const double dsvtuZ=+(Cay[nss(m,n)]-Cay[nss(m,m)]); /* plus sign of minor v=n */
1578 sum1-=sumZ*dsvtuZ;
1579
1580 sum1/=9*dstustu;
1581 }
1582 }
1583 else { assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1583, __extension__ __PRETTY_FUNCTION__));
1584 sum1=-(Cay[nss(m,m)]+Cay[nss(m,n)]+Cay[nss(n,n)])/6.;
1585 }
1586 pI2Dstu[ep][idx]=sum1;
1587}
1588
1589/* --------------------------------------------------------
1590 * I3D2st triangle in D+4 dim
1591 * --------------------------------------------------------
1592 */
1593ncomplex Minor5::I3D2st(int ep, int s, int t)
1594{
1595 assert(s!=t)(static_cast <bool> (s!=t) ? void (0) : __assert_fail (
"s!=t", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1595, __extension__ __PRETTY_FUNCTION__));
1596 if (ep==2) return 0;
1597 if (not fEval[E_I3D2st+ep]) {
1598 I3D2stEval(ep);
1599 }
1600 int idx = im2(s,t)-5;
1601 return pI3D2st[ep][idx];
1602}
1603
1604void Minor5::I3D2stEval(int ep)
1605{
1606 for (int s=1; s<=smax; s++) {
1607 for (int t=s+1; t<=5; t++) {
1608 int idx = im2(s,t)-5;
1609 ncomplex ivalue=0;
1610
1611 if (ep==0) {
1612 const double dstst=M2(s,t,s,t);
1613 const double d0st0st=M3(0,s,t,0,s,t);
1614
1615 if ( pmaxT3[idx]!=0 && ( pmaxT3[idx] <= epsir1 && pmaxU3[idx] <= epsir1 ) ) {
1616 // IR
1617 int i=imax3[idx];
1618 int iu[3]={i-1,s-1,t-1};
1619 int tmp;
1620 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
1621 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
1622 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
1623 int nu[3];
1624 freeidxM3(iu, nu);
1625 int u=nu[0]+1;
1626 int v=nu[1]+1;
1627 const double Dii=M4ii(u,v,i);
1628 const double Dui=M4ui(u,v,i);
1629 const double Dvi=M4vi(u,v,i);
1630 ncomplex sum1=+Dii*(I2Dstu(0, s, t, i)+0.5*I2Dstu(1, s, t, i)) // (i, i)
1631 +Dui*(I2Dstu(0, s, t, u)+0.5*I2Dstu(1, s, t, u)) // (u, i)
1632 +Dvi*(I2Dstu(0, s, t, v)+0.5*I2Dstu(1, s, t, v)); // (v, i)
1633 ivalue=0.25*sum1/M3(0,s,t,i,s,t);
1634 } else if (pmaxS3[idx] <= ceps) {
1635 // EXP
1636 const double x=dstst/d0st0st;
1637 ncomplex sump;
1638 do {
1639 double dstust0[6];
1640 ncomplex sum1=0;
1641 for (int u=1; u<=5; u++) {
1642 if (u==t || u==s) continue;
1643 dstust0[u]=M3(s,t,u,s,t,0);
1644 sum1+=dstust0[u]*I2D2stu(0, s, t, u);
1645 }
1646
1647 double xn=1;
1648 ncomplex dv,s21;
1649
1650 ncomplex sum[3];
1651 sum[0]=sump=sum1;
1652
1653#define stepI3D(n,a,b) \
1654 xn*=x; \
1655 dv=0; \
1656 for (int u=1; u<=5; u++) { \
1657 if (u==t || u==s) continue; \
1658 dv+=dstust0[u]*(a*I2D##n##stu(0, s, t, u) - b*I2D##n##stu(1, s, t, u)); \
1659 } \
1660 dv*=xn; \
1661 sum1+=dv;
1662
1663 stepI3D(3,6.,2.)
1664 if ( fabs(sum1.real()*teps)>=fabs(dv.real())
1665 && fabs(sum1.imag()*teps)>=fabs(dv.imag()))
1666 break;
1667 sum[1]=sump=sum1;
1668 s21=sum[1]-sum[0];
1669
1670 stepI3D(4,48.,28.)
1671 sump=sum1;
1672 stepWynn(0)sum[(2+0)%3]=sum1; { const ncomplex s2=sum[(2+0)%3]; const ncomplex
s1=sum[(1+0)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1673 stepI3D(5,480.,376.)
1674 stepWynn(1)sum[(2+1)%3]=sum1; { const ncomplex s2=sum[(2+1)%3]; const ncomplex
s1=sum[(1+1)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1675 stepI3D(6,5760.,5472.)
1676 stepWynn(2)sum[(2+2)%3]=sum1; { const ncomplex s2=sum[(2+2)%3]; const ncomplex
s1=sum[(1+2)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1677// stepI3D(7,80640.,88128.)
1678// stepWynn(3)
1679// stepI3D(8,1290240.,1571328.)
1680// stepWynn(4)
1681#undef stepI3D
1682 } while (0);
1683 ivalue=sump/d0st0st;
1684 } else {
1685 // NORMAL
1686 ncomplex sum1=0;
1687 for (int u=1; u<=5; u++) {
1688 if (u==t || u==s) continue;
1689 sum1-=M3(u,s,t,0,s,t)*I2Dstu(ep, s, t, u);
1690 }
1691 sum1+=d0st0st*I3Dst(ep, s, t);
1692 ivalue=sum1/(4*dstst)+I3D2st(ep+1, s, t)*0.5; // 2*x/4 == 0.5*x
1693 }
1694 } else {
1695 assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1695, __extension__ __PRETTY_FUNCTION__));
1696 int iu[3]={0,s,t};
1697 int nu[3];
1698 freeidxM3(iu, nu);
1699 ivalue=(Cay[nss(nu[0],nu[0])]+Cay[nss(nu[1],nu[1])]+Cay[nss(nu[2],nu[2])]
1700 +Cay[nss(nu[0],nu[1])]+Cay[nss(nu[0],nu[2])]+Cay[nss(nu[1],nu[2])])/24.;
1701 }
1702 pI3D2st[ep][idx]=ivalue;
1703 }
1704 }
1705 fEval[E_I3D2st+ep]=true;
1706}
1707
1708/* --------------------------------------------------------
1709 * I4D3s box in D+6 dim
1710 * --------------------------------------------------------
1711 */
1712ncomplex Minor5::I4D3s(int ep, int s)
1713{
1714 if (ep==2) return 0;
1715 if (not fEval[E_I4D3s+ep]) {
1716 I4D3sEval(ep);
1717 }
1718 return pI4D3s[ep][s-1];
1719}
1720void Minor5::I4D3sEval(int ep)
1721{
1722 for (int s=1; s<=smax; s++) {
1723 const double dss=M1(s, s);
1724 const double d0s0s=M2(0, s, 0, s);
1725 ncomplex ivalue=0;
1726
1727 if (ep==0) {
1728 ncomplex sum1=0;
1729 for (int t=1; t<=5; t++) {
1730 if (t==s) continue;
1731 sum1-=M2(s,t,s,0)*I3D2st(ep, s, t);
1732 }
1733 sum1+=d0s0s*I4D2s(ep, s);
1734 ivalue=sum1/(5*dss)+2.*I4D3s(ep+1, s)/5.;
1735
1736 const double x=dss/d0s0s;
1737 if (pmaxS4[s-1] <= deps3) {
1738 ncomplex sump;
1739 do {
1740 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1740, __extension__ __PRETTY_FUNCTION__));
1741
1742 double dsts0[6];
1743 sum1=0;
1744 for (int t=1; t<=5; t++) {
1745 if (t==s) continue;
1746 dsts0[t]=M2(s,t,s,0);
1747 sum1+=dsts0[t]*I3D3st(0, s, t);
1748 }
1749
1750 double xn=1;
1751 ncomplex dv,s21;
1752
1753 ncomplex sum[3];
1754 sum[0]=sump=sum1;
1755
1756#define stepI4D(n,a,b) \
1757 xn*=x; \
1758 dv=0; \
1759 for (int t=1; t<=5; t++) { \
1760 if (t==s) continue; \
1761 dv+=dsts0[t]*(a*I3D##n##st(0, s, t) - b*I3D##n##st(1, s, t)); \
1762 } \
1763 dv*=xn; \
1764 sum1+=dv;
1765
1766 stepI4D(4,7.,2.)
1767 if ( fabs(sum1.real()*teps)>=fabs(dv.real())
1768 && fabs(sum1.imag()*teps)>=fabs(dv.imag()))
1769 break;
1770 sum[1]=sump=sum1;
1771 s21=sum[1]-sum[0];
1772
1773 stepI4D(5,63.,32.)
1774 sump=sum1;
1775 stepWynn(0)sum[(2+0)%3]=sum1; { const ncomplex s2=sum[(2+0)%3]; const ncomplex
s1=sum[(1+0)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1776 stepI4D(6,693.,478.)
1777 stepWynn(1)sum[(2+1)%3]=sum1; { const ncomplex s2=sum[(2+1)%3]; const ncomplex
s1=sum[(1+1)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1778 stepI4D(7,9009.,7600.)
1779 stepWynn(2)sum[(2+2)%3]=sum1; { const ncomplex s2=sum[(2+2)%3]; const ncomplex
s1=sum[(1+2)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
1780// stepI4D(8,135135.,132018.)
1781// stepWynn(3)
1782// stepI4D(9,2297295.,2514576.)
1783// stepWynn(4)
1784// stepI4D(10,43648605.,52371534.)
1785// stepWynn(5)
1786#undef stepI4D
1787 } while (0);
1788 ivalue=sump/d0s0s;
1789 }
1790 }
1791 else {
1792 assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1792, __extension__ __PRETTY_FUNCTION__));
1793 double sum1=0;
1794 for (int i=1; i<=5; i++) { if (i==s) continue;
1795 for (int j=i; j<=5; j++) { if (j==s) continue;
1796 sum1+=Cay[nss(i,j)];
1797 }
1798 }
1799 ivalue=-sum1/120.;
1800 }
1801 pI4D3s[ep][s-1]=ivalue;
1802 }
1803 fEval[E_I4D3s+ep]=true;
1804}
1805
1806/* --------------------------------------------------------
1807 * I3D2sti triangle in D+4 dim with a dot
1808 * --------------------------------------------------------
1809 */
1810ncomplex Minor5::I3D2sti(int ep, int s, int t, int i)
1811{
1812 assert(s!=t && s!=i && t!=i)(static_cast <bool> (s!=t && s!=i && t!=
i) ? void (0) : __assert_fail ("s!=t && s!=i && t!=i"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 1812
, __extension__ __PRETTY_FUNCTION__));
1813 if (ep==1) return 1./6.;
1814 else if (ep==2) return 0.;
1815 if (not fEval[E_I3D2sti+ep]) {
1816 I3D2stiEval(ep);
1817 }
1818 int idx = im2(s,t)-5;
1819 return pI3D2sti[ep][i-1][idx];
1820}
1821
1822void Minor5::I3D2stiEval(int ep)
1823{
1824 for (int i=1; i<=CIDX(DCay-2); i++) {
1825 for (int s=1; s<=smax; s++) { if (i==s) continue;
1826 for (int t=s+1; t<=5; t++) { if (i==t) continue;
1827 int idx = im2(s,t)-5;
1828 ncomplex ivalue=0;
1829
1830 if ( (pmaxT3[idx]==0 || (pmaxT3[idx] > epsir2 || pmaxU3[idx] > epsir2))
1831 && pmaxS3[idx] > ceps ) {
1832 // Variant with Gram3
1833 ncomplex sum1=0;
1834 for (int u=1; u<=5; u++) {
1835 if (u==t || u==s) continue;
1836 sum1+=M3(u,s,t,i,s,t)*I2Dstu(ep,s,t,u);
1837 }
1838 sum1-=M3(0,s,t,i,s,t)*I3Dst(ep, s, t);
1839 ivalue=sum1/M2(s,t,s,t);
1840 } else {
1841 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1841, __extension__ __PRETTY_FUNCTION__));
1842 int iu[3]={i-1,s-1,t-1};
1843 int tmp;
1844 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
1845 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
1846 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
1847 int nu[3];
1848 freeidxM3(iu, nu);
1849 int u=nu[0]+1;
1850 int v=nu[1]+1;
1851
1852 if ( pmaxT3[idx] <= epsir2 && pmaxU3[idx] <= epsir2 ) {
1853 // small G3 & C3
1854 int j=imax3[idx];
1855 ncomplex sum1=0;
1856 ncomplex const I3term=I3Dst(ep,s,t)-1./3.; //+2./3.*I3Dst(ep+1,s,t))
1857 ncomplex const I2Uterm=I2Dstui(ep, s, t, u, i)+2./3.*I2Dstui(ep+1, s, t, u, i);
1858 ncomplex const I2Vterm=I2Dstui(ep, s, t, v, i)+2./3.*I2Dstui(ep+1, s, t, v, i);
1859
1860 if (j==i) { // j->i
1861 const double Dii=M4ii(u,v,i);
1862 const double Dui=M4ui(u,v,i);
1863 const double Dvi=M4vi(u,v,i);
1864 sum1+=+Dii*(I3term) // (i, i)
1865 +Dui*(I2Uterm) // (u, i)
1866 +Dvi*(I2Vterm); // (v, i)
1867 } else if (j==u) { // j->u
1868 const double Dui=M4ui(u,v,i);
1869 const double Duu=M4uu(u,v,i);
1870 const double Dvu=M4vu(u,v,i);
1871 sum1+=+Dui*(I3term) // (u, i)
1872 +Duu*(I2Uterm) // (u, u)
1873 +Dvu*(I2Vterm); // (v, u)
1874 } else { assert(j==v)(static_cast <bool> (j==v) ? void (0) : __assert_fail (
"j==v", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 1874, __extension__ __PRETTY_FUNCTION__)); // j->v
1875 const double Dvi=M4vi(u,v,i);
1876 const double Dvv=M4vv(u,v,i);
1877 const double Dvu=M4vu(u,v,i);
1878 sum1+=+Dvi*(I3term) // (v, i)
1879 +Dvu*(I2Uterm) // (v, u)
1880 +Dvv*(I2Vterm); // (v, v)
1881 }
1882 ivalue=sum1/(3*M3(s,0,t,s,j,t));
1883 } else {
1884 // small G3
1885 const double Dii=M4ii(u,v,i);
1886 const double Dui=M4ui(u,v,i);
1887 const double Dvi=M4vi(u,v,i);
1888 ncomplex sum1=0;
1889 sum1+=Dii*I2Dstu(ep, s, t, i) // (i, i)
1890 +Dui*I2Dstu(ep, s, t, u) // (u, i)
1891 +Dvi*I2Dstu(ep, s, t, v); // (v, i)
1892 sum1+=M3(s,0,t,s,i,t)*(-4.*I3D2st(ep, s, t)+2.*I3D2st(ep+1, s, t));
1893 ivalue=sum1/M3(s,0,t,s,0,t);
1894 }
1895 }
1896 pI3D2sti[ep][i-1][idx]=ivalue;
1897 }
1898 }
1899 }
1900 fEval[E_I3D2sti+ep]=true;
1901}
1902
1903/* --------------------------------------------------------
1904 * I4D3si box in D+6 dim with a dot
1905 * --------------------------------------------------------
1906 */
1907ncomplex Minor5::I4D3si(int ep, int s, int i)
1908{
1909 if (s==i) return 0;
1910 if (ep==1) return -1./24.;
1911 else if (ep==2) return 0;
1912 if (not fEval[E_I4D3si+ep]) {
1913 I4D3siEval(ep);
1914 }
1915 return pI4D3si[ep][i-1][s-1];
1916}
1917
1918void Minor5::I4D3siEval(int ep)
1919{
1920 for (int s=1; s<=smax; s++) {
1921 for (int i=1; i<=CIDX(DCay-2); i++) {
1922 if (s==i) continue;
1923 ncomplex ivalue=0;
1924
1925 if (pmaxS4[s-1] <= deps3) {
1926 ncomplex sum1=0;
1927 for (int t=1; t<=5; t++) {
1928 if (t==s) continue;
1929 sum1+=M3(0, s, t, 0, s, i)*I3D2st(ep, s, t);
1930 }
1931 sum1+=M2(0, s, i, s)*(-5.*I4D3s(ep, s)+2.*I4D3s(ep+1, s));
1932 ivalue=sum1/M2(0, s, 0, s);
1933 } else {
1934 ncomplex sum1=0;
1935 for (int t=1; t<=5; t++) {
1936 if (t==s) continue;
1937 sum1+=M2(s, t, s, i)*I3D2st(ep, s, t);
1938 }
1939 sum1-=M2(0, s, i, s)*I4D2s(ep, s);
1940 ivalue=sum1/M1(s, s);
1941 }
1942 pI4D3si[ep][i-1][s-1]=ivalue;
1943 }
1944 }
1945 fEval[E_I4D3si+ep]=true;
1946}
1947
1948/* --------------------------------------------------------
1949 * I4D3sij box in D+6 dim with two dots
1950 * --------------------------------------------------------
1951 */
1952ncomplex Minor5::I4D3sij(int ep, int s, int i, int j)
1953{
1954 if (s==i || s==j) return 0;
1955 else if (ep!=0) return 0; // I4D3sij is finite
1956 if (not fEval[E_I4D3sij+ep]) {
1957 I4D3sijEval(ep);
1958 }
1959 return pI4D3sij[ep][is(i-1,j-1)][s-1];
1960}
1961
1962void Minor5::I4D3sijEval(int ep)
1963{
1964 for (int s=1; s<=smax; s++) {
1965 // symmetric in 'i,j'
1966 for (int i=1; i<=CIDX(DCay-2); i++) { if (s==i) continue;
1967 for (int j=i; j<=CIDX(DCay-2); j++) { if (s==j) continue;
1968 ncomplex ivalue=0;
1969
1970 if (pmaxS4[s-1] <= deps3) {
1971 ncomplex sum1=0;
1972 for (int t=1; t<=5; t++) {
1973 if (t==s || t==i) continue;
1974 sum1+=M3(s,0,t,s,0,j)*I3D2sti(ep, s, t, i);
1975 }
1976 sum1+=M3(s,0,i,s,0,j)*I4D2s(ep, s);
1977 sum1+=M2(s,0,s,j)*(-4.*I4D3si(ep, s, i)+2.*I4D3si(ep+1, s, i));
1978 ivalue=sum1/M2(0,s,0,s);
1979 } else {
1980 ncomplex sum1=0;
1981 for (int t=1; t<=5; t++) {
1982 if (t==s || t==i) continue;
1983 sum1+=M2(s,t,s,j)*I3D2sti(ep, s, t, i);
1984 }
1985 sum1+=M2(s,i,s,j)*I4D2s(ep, s);
1986 sum1-=M2(s,0,s,j)*I4D2si(ep, s, i);
1987 ivalue=sum1/M1(s,s);
1988 }
1989 pI4D3sij[ep][iss(i-1,j-1)][s-1]=ivalue;
1990 }
1991 }
1992 }
1993 fEval[E_I4D3sij+ep]=true;
1994}
1995
1996
1997/* --------------------------------------------------------
1998 * I2Dstui bubble in D+2 dim with a dot
1999 * --------------------------------------------------------
2000 */
2001ncomplex Minor5::I2Dstui(int ep, int s, int t, int u, int i)
2002{
2003 assert(s!=t && t!=u && u!=s && s!=i && t!=i && u!=i)(static_cast <bool> (s!=t && t!=u && u!=
s && s!=i && t!=i && u!=i) ? void (0)
: __assert_fail ("s!=t && t!=u && u!=s && s!=i && t!=i && u!=i"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 2003
, __extension__ __PRETTY_FUNCTION__));
2004// if (ep==1) return -0.5; // not quite true
2005 if (ep==2) return 0;
2006 if (not fEval[E_I2Dstui+ep]) {
2007 I2DstuiEval(ep,1,4,5,2,3,kinem.p3());
2008 I2DstuiEval(ep,1,3,5,2,4,kinem.s34());
2009 I2DstuiEval(ep,1,3,4,2,5,kinem.s12());
2010 I2DstuiEval(ep,1,4,5,3,2,kinem.p3());
2011 I2DstuiEval(ep,1,2,5,3,4,kinem.p4());
2012 I2DstuiEval(ep,1,2,4,3,5,kinem.s45());
2013 I2DstuiEval(ep,1,3,5,4,2,kinem.s34());
2014 I2DstuiEval(ep,1,2,5,4,3,kinem.p4());
2015 I2DstuiEval(ep,1,2,3,4,5,kinem.p5());
2016#ifdef USE_ZERO_CHORD
2017 I2DstuiEval(ep,1,3,4,5,2,kinem.s12());
2018 I2DstuiEval(ep,1,2,4,5,3,kinem.s45());
2019 I2DstuiEval(ep,1,2,3,5,4,kinem.p5());
2020#endif
2021
2022 if (smax==5) {
2023 I2DstuiEval(ep,3,4,5,1,2,kinem.p2());
2024 I2DstuiEval(ep,2,4,5,1,3,kinem.s23());
2025 I2DstuiEval(ep,2,3,5,1,4,kinem.s15());
2026 I2DstuiEval(ep,2,3,4,1,5,kinem.p1());
2027 I2DstuiEval(ep,3,4,5,2,1,kinem.p2());
2028 I2DstuiEval(ep,2,4,5,3,1,kinem.s23());
2029 I2DstuiEval(ep,2,3,5,4,1,kinem.s15());
2030#ifdef USE_ZERO_CHORD
2031 I2DstuiEval(ep,2,3,4,5,1,kinem.p1());
2032#endif
2033 }
2034
2035 fEval[E_I2Dstui+ep]=true;
2036 }
2037 int ip=15-s-t-u-i;
2038 return pI2Dstui[ep][i-1][ip-1];
2039}
2040
2041void Minor5::I2DstuiEval(int ep, int s, int t, int u, int i, int ip, double qsq)
2042{
2043 ncomplex sum1=0;
2044 if (ep==0) {
2045 const double dstustu=-2*qsq; /*M3(s,t,u,s,t,u);*/
2046 const double msq1=kinem.mass(i);
2047 const double msq2=kinem.mass(ip);
2048 const double s_cutoff=seps1*pmaxM2[im2(i,ip)-5];
2049
2050 if (fabs(dstustu) <= s_cutoff) {
2051 const double mm12=msq1-msq2;
2052 if (fabs(mm12) < meps) {
2053 if (msq1 > meps) {
2054 sum1=(msq1 - ICache::getI1(ep, Kinem1(msq1)))/(2*msq1);
2055 } else {
2056 sum1=0;
2057 }
2058 }
2059 else {
2060 sum1=(msq1 + msq2)/(4*msq1 - 4*msq2)
2061 - ((msq1 - 2*msq2)*ICache::getI1(ep, Kinem1(msq1))
2062 + msq2*ICache::getI1(ep, Kinem1(msq2))
2063 )/(2*mm12*mm12);
2064 }
2065 }
2066 else {
2067 sum1+=+(Cay[nss(ip,ip)]-Cay[ns(i,ip)])*I2stu(ep,s,t,u);
2068 sum1+=+ICache::getI1(ep, Kinem1(msq1));
2069 sum1+=-ICache::getI1(ep, Kinem1(msq2));
2070 sum1/=dstustu;
2071 }
2072 }
2073 else { assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2073, __extension__ __PRETTY_FUNCTION__));
2074 if ( fabs(qsq) < meps
2075 && fabs(kinem.mass(i)) < meps
2076 && fabs(kinem.mass(ip)) < meps ) {
2077 sum1=0;
2078 }
2079 else {
2080 sum1=-0.5;
2081 }
2082 }
2083 pI2Dstui[ep][i-1][ip-1]=sum1;
2084}
2085
2086/* --------------------------------------------------------
2087 * I3D2stij triangle in D+4 dim with two dots
2088 * --------------------------------------------------------
2089 */
2090ncomplex Minor5::I3D2stij(int ep, int s, int t, int i, int j) // IR-div
2091{
2092 assert(s!=t && s!=i && s!=j && t!=i && t!=j)(static_cast <bool> (s!=t && s!=i && s!=
j && t!=i && t!=j) ? void (0) : __assert_fail
("s!=t && s!=i && s!=j && t!=i && t!=j"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 2092
, __extension__ __PRETTY_FUNCTION__));
2093 if (not fEval[E_I3D2stij+ep]) {
2094 I3D2stijEval(ep);
2095 }
2096 int idx = im2(s,t)-5;
2097 return pI3D2stij[ep][is(i-1,j-1)][idx];
2098}
2099
2100void Minor5::I3D2stijEval(int ep)
2101{
2102 for (int s=1; s<=smax; s++) {
2103 for (int t=s+1; t<=5; t++) {
2104 int idx = im2(s,t)-5;
2105
2106 const double ds0ts0t=M3(s,0,t,s,0,t);
2107 if (ep!=0 && fabs(ds0ts0t) > m3eps) { // if ds0ts0t!=0 I3D2stij is finite
2108 for (int ij=iss(1-1,1-1); ij<=iss(CIDX(DCay-2)-1,CIDX(DCay-2)-1); ij++) {
2109 pI3D2stij[ep][ij][idx]=0;
2110 }
2111 continue;
2112 }
2113
2114 const double dstst=M2(s,t,s,t);
2115 // symmetric in 'i,j'
2116 for (int i=1; i<=CIDX(DCay-2); i++) { if (i==s || i==t) continue;
2117 for (int j=i; j<=CIDX(DCay-2); j++) { if (j==s || j==t) continue;
2118 ncomplex ivalue=0;
2119
2120 if (pmaxS3[idx] > ceps || ep!=0) {
2121 // Variant with Gram3
2122 ncomplex sum1=0;
2123 for (int u=1; u<=5; u++) {
2124 if (u==t || u==s || u==i) continue;
2125 sum1+=M3(s,u,t,s,j,t)*I2Dstui(ep, s, t, u, i);
2126 }
2127 sum1+=-M3(s,0,t,s,j,t)*I3Dsti(ep, s, t, i)+M3(s,i,t,s,j,t)*I3Dst(ep, s, t);
2128 ivalue=sum1/dstst;
2129 } else {
2130 ncomplex sum1=0;
2131 int iu[3]={j-1,s-1,t-1};
2132 int tmp;
2133 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
2134 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
2135 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
2136 int nu[3];
2137 freeidxM3(iu, nu);
2138 int u=nu[0]+1;
2139 int v=nu[1]+1;
2140 const double Djj=M4ii(u,v,j);
2141 const double Duj=M4ui(u,v,j);
2142 const double Dvj=M4vi(u,v,j);
2143 if ( fabs(ds0ts0t) > 0. ) {
2144 if (j==i) {
2145 sum1+=+Djj*I3Dst(ep,s,t)
2146 +Duj*I2Dstui(ep, s, t, u, i)
2147 +Dvj*I2Dstui(ep, s, t, v, i);
2148 } else {
2149 sum1+=Djj*I2Dstui(ep, s, t, j, i);
2150 if (i==u) {
2151 sum1+=+Duj*I3Dst(ep,s,t)
2152 +Dvj*I2Dstui(ep, s, t, v, i);
2153 } else {
2154 sum1+=+Dvj*I3Dst(ep,s,t)
2155 +Duj*I2Dstui(ep, s, t, u, i);
2156 }
2157 }
2158 if (ep<2)
2159 sum1+=M3(s,0,t,s,j,t)*(-3.*I3D2sti(ep, s, t, i)+2.*I3D2sti(ep+1, s, t, i));
2160 else
2161 sum1+=M3(s,0,t,s,j,t)*(-3.*I3D2sti(ep, s, t, i));
2162 ivalue=sum1/ds0ts0t;
2163 } else {
2164 ivalue=std::numeric_limits<double>::quiet_NaN();
2165 // TODO add: need I2Dstuij and I2stui
2166 }
2167 }
2168 pI3D2stij[ep][iss(i-1,j-1)][idx]=ivalue;
2169 }
2170 }
2171 }
2172 }
2173 fEval[E_I3D2stij+ep]=true;
2174}
2175
2176/* --------------------------------------------------------
2177 * I4D3sijk box in D+6 dim with three dots
2178 * --------------------------------------------------------
2179 */
2180ncomplex Minor5::I4D3sijk(int ep, int s, int i, int j, int k) // IR-div
2181{
2182 if (s==i || s==j || s==k) return 0;
2183 if (not fEval[E_I4D3sijk+ep]) {
2184 I4D3sijkEval(ep);
2185 }
2186 return pI4D3sijk[ep][is(i-1,j-1,k-1)][s-1];
2187}
2188
2189void Minor5::I4D3sijkEval(int ep)
2190{
2191 for (int s=1; s<=smax; s++) {
2192 // symmetric in 'i,j,k'
2193 for (int i=1; i<=CIDX(DCay-2); i++) { if (i==s) continue;
2194 for (int j=i; j<=CIDX(DCay-2); j++) { if (j==s) continue;
2195 for (int k=j; k<=CIDX(DCay-2); k++) { if (k==s) continue;
2196 ncomplex ivalue=0;
2197
2198 if (pmaxS4[s-1] <= deps3) {
2199 ncomplex sum1=0;
2200 for (int t=1; t<=5; t++) {
2201 if (s==t || t==i || t==j) continue;
2202 sum1+=M3(s,0,t,s,0,k)*I3D2stij(ep,s,t,i,j);
2203 }
2204 sum1+=+M3(s,0,i,s,0,k)*I4D2si(ep, s, j)
2205 +M3(s,0,j,s,0,k)*I4D2si(ep, s, i);
2206 if (ep<2) {
2207 sum1+=M2(s,0,s,k)*(-3.*I4D3sij(ep, s, i, j)+2.*I4D3sij(ep+1, s, i, j));
2208 }
2209 else { // ep==2
2210 sum1+=M2(s,0,s,k)*(-3.*I4D3sij(ep, s, i, j));
2211 }
2212 ivalue=sum1/M2(s,0,s,0);
2213 } else {
2214 ncomplex sum1=0;
2215 for (int t=1; t<=5; t++) {
2216 if (t==s || t==i || t==j) continue;
2217 sum1+=M2(s,t,s,k)*I3D2stij(ep,s,t,i,j);
2218 }
2219 sum1-=M2(s,0,s,k)*I4D2sij(ep,s,i,j);
2220 sum1+=M2(s,i,s,k)*I4D2si(ep,s,j)+M2(s,j,s,k)*I4D2si(ep,s,i);
2221 ivalue=sum1/M1(s,s);
2222 }
2223 pI4D3sijk[ep][iss(i-1,j-1,k-1)][s-1]=ivalue;
2224 }
2225 }
2226 }
2227 }
2228 fEval[E_I4D3sijk+ep]=true;
2229}
2230
2231/* --------------------------------------------------------
2232 *
2233 * Rank-5 functions
2234 *
2235 * --------------------------------------------------------
2236 */
2237
2238/* --------------------------------------------------------
2239 * I2D2stu bubble in D+4 dim
2240 * --------------------------------------------------------
2241 */
2242ncomplex Minor5::I2D2stu(int ep, int s, int t, int u)
2243{
2244 assert(t!=u && u!=s && s!=t)(static_cast <bool> (t!=u && u!=s && s!=
t) ? void (0) : __assert_fail ("t!=u && u!=s && s!=t"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 2244
, __extension__ __PRETTY_FUNCTION__));
2245 if (ep==2) return 0;
2246 if (not fEval[E_I2D2stu+ep]) {
2247 I2D2stuEval(0,ep,1,2,3,4,5,kinem.p5());
2248 I2D2stuEval(1,ep,1,2,4,3,5,kinem.s45());
2249 I2D2stuEval(2,ep,1,2,5,3,4,kinem.p4());
2250
2251 I2D2stuEval(3,ep,1,3,4,2,5,kinem.s12());
2252 I2D2stuEval(4,ep,1,3,5,2,4,kinem.s34());
2253
2254 I2D2stuEval(5,ep,1,4,5,2,3,kinem.p3());
2255
2256 if (smax==5) {
2257 I2D2stuEval(6,ep,2,3,4,1,5,kinem.p1());
2258 I2D2stuEval(7,ep,2,3,5,1,4,kinem.s15());
2259
2260 I2D2stuEval(8,ep,2,4,5,1,3,kinem.s23());
2261
2262
2263 I2D2stuEval(9,ep,3,4,5,1,2,kinem.p2());
2264 }
2265
2266 fEval[E_I2D2stu+ep]=true;
2267 }
2268 int idx=im3(s,t,u)-10;
2269 return pI2D2stu[ep][idx];
2270}
2271
2272void Minor5::I2D2stuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq)
2273{
2274 ncomplex sum1=0;
2275 if (ep==0) {
2276 const double dstustu=-2*qsq; /*M3(s,t,u,s,t,u);*/
2277 const double msq1=kinem.mass(m);
2278 const double msq2=kinem.mass(n);
2279 const double s_cutoff=seps1*pmaxM2[im2(m,n)-5];
2280
2281 if (fabs(dstustu) <= s_cutoff) {
2282 const double mm12=msq1-msq2;
2283 if (fabs(mm12) < meps) {
2284 sum1=0.25*msq1*(msq1 + 2.*ICache::getI1(ep, Kinem1(msq1)));
2285 }
2286 else {
2287 sum1= 5*(msq1*msq1 + msq1*msq2 + msq2*msq2)/36
2288 + ( + msq1*msq1*ICache::getI1(ep, Kinem1(msq1))
2289 - msq2*msq2*ICache::getI1(ep, Kinem1(msq2))
2290 )/(6*mm12);
2291 }
2292 }
2293 else {
2294 ncomplex sumX=5.*I2Dstu(ep,s,t,u)+2.*I2Dstu(ep+1,s,t,u);
2295 ncomplex sumY=-0.25*msq2*(10.*ICache::getI1(ep, Kinem1(msq2))+9*msq2);
2296 ncomplex sumZ=-0.25*msq1*(10.*ICache::getI1(ep, Kinem1(msq1))+9*msq1);
2297
2298 const double ds0tu=(Cay[nss(m,m)]*Cay[nss(n,n)]-Cay[nss(m,n)]*Cay[nss(m,n)]);
2299 sum1+=sumX*ds0tu;
2300
2301 const double dsvtuY=-(Cay[nss(n,n)]-Cay[nss(m,n)]); /* minus sign of minor v=m */
2302 sum1-=sumY*dsvtuY;
2303
2304 const double dsvtuZ=+(Cay[nss(m,n)]-Cay[nss(m,m)]); /* plus sign of minor v=n */
2305 sum1-=sumZ*dsvtuZ;
2306
2307 sum1/=25*dstustu;
2308 }
2309 }
2310 else { assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2310, __extension__ __PRETTY_FUNCTION__));
2311 const double y11=Cay[nss(m,m)];
2312 const double y12=Cay[nss(m,n)];
2313 const double y22=Cay[nss(n,n)];
2314 sum1= ( 3*( y11*(y11+y12)+(y12+y22)*y22)+2*y12*y12+y11*y22 )/120;
2315 }
2316 pI2D2stu[ep][idx]=sum1;
2317}
2318
2319/* --------------------------------------------------------
2320 * I3D3st triangle in D+6 dim
2321 * --------------------------------------------------------
2322 */
2323ncomplex Minor5::I3D3st(int ep, int s, int t)
2324{
2325 assert(s!=t)(static_cast <bool> (s!=t) ? void (0) : __assert_fail (
"s!=t", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2325, __extension__ __PRETTY_FUNCTION__));
2326 if (ep==2) return 0;
2327 if (not fEval[E_I3D3st+ep]) {
2328 I3D3stEval(ep);
2329 }
2330 int idx = im2(s,t)-5;
2331 return pI3D3st[ep][idx];
2332}
2333
2334void Minor5::I3D3stEval(int ep)
2335{
2336 for (int s=1; s<=smax; s++) {
2337 for (int t=s+1; t<=5; t++) {
2338 int idx = im2(s,t)-5;
2339 ncomplex ivalue=0;
2340
2341 if (ep==0) {
2342 const double dstst=M2(s,t,s,t);
2343 const double d0st0st=M3(0,s,t,0,s,t);
2344
2345 if ( pmaxT3[idx]!=0 && ( pmaxT3[idx] <= epsir1 && pmaxU3[idx] <= epsir1 ) ) {
2346 // IR
2347 int i=imax3[idx];
2348 int iu[3]={i-1,s-1,t-1};
2349 int tmp;
2350 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
2351 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
2352 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
2353 int nu[3];
2354 freeidxM3(iu, nu);
2355 int u=nu[0]+1;
2356 int v=nu[1]+1;
2357 const double Dii=M4ii(u,v,i);
2358 const double Dui=M4ui(u,v,i);
2359 const double Dvi=M4vi(u,v,i);
2360 ncomplex sum1=+Dii*(I2D2stu(0, s, t, i)+I2D2stu(1, s, t, i)/3.) // (i, i)
2361 +Dui*(I2D2stu(0, s, t, u)+I2D2stu(1, s, t, u)/3.) // (u, i)
2362 +Dvi*(I2D2stu(0, s, t, v)+I2D2stu(1, s, t, v)/3.); // (v, i)
2363 ivalue=sum1/(6*M3(0,s,t,i,s,t));
2364 } else if (pmaxS3[idx] <= ceps) {
2365 // EXP
2366 const double x=dstst/d0st0st;
2367 ncomplex sump;
2368 do {
2369 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2369, __extension__ __PRETTY_FUNCTION__));
2370
2371 double dstust0[6];
2372 ncomplex sum1=0;
2373 for (int u=1; u<=5; u++) {
2374 if (u==t || u==s) continue;
2375 dstust0[u]=M3(s,t,u,s,t,0);
2376 sum1+=dstust0[u]*I2D3stu(0, s, t, u);
2377 }
2378
2379 double xn=1;
2380 ncomplex dv,s21;
2381
2382 ncomplex sum[3];
2383 sum[0]=sump=sum1;
2384
2385#define stepI3D(n,a,b) \
2386 xn*=x; \
2387 dv=0; \
2388 for (int u=1; u<=5; u++) { \
2389 if (u==t || u==s) continue; \
2390 dv+=dstust0[u]*(a*I2D##n##stu(0, s, t, u) - b*I2D##n##stu(1, s, t, u)); \
2391 } \
2392 dv*=xn; \
2393 sum1+=dv;
2394
2395 stepI3D(4,8.,2.)
2396 if ( fabs(sum1.real()*teps)>=fabs(dv.real())
2397 && fabs(sum1.imag()*teps)>=fabs(dv.imag()))
2398 break;
2399 sum[1]=sump=sum1;
2400 s21=sum[1]-sum[0];
2401
2402 stepI3D(5,80.,36.)
2403 sump=sum1;
2404 stepWynn(0)sum[(2+0)%3]=sum1; { const ncomplex s2=sum[(2+0)%3]; const ncomplex
s1=sum[(1+0)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
2405 stepI3D(6,960.,592.)
2406 stepWynn(1)sum[(2+1)%3]=sum1; { const ncomplex s2=sum[(2+1)%3]; const ncomplex
s1=sum[(1+1)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
2407// stepI3D(7,13440.,10208.)
2408// stepWynn(2)
2409// stepI3D(8,215040.,190208.)
2410// stepWynn(3)
2411// stepI3D(9,3870720.,3853824.)
2412// stepWynn(4)
2413#undef stepI3D
2414 } while (0);
2415 ivalue=sump/d0st0st;
2416 } else {
2417 // NORMAL
2418 ncomplex sum1=0;
2419 for (int u=1; u<=5; u++) {
2420 if (u==t || u==s) continue;
2421 sum1-=M3(u,s,t,0,s,t)*I2D2stu(ep, s, t, u);
2422 }
2423 sum1+=d0st0st*I3D2st(ep, s, t);
2424 ivalue=sum1/(6*dstst)+I3D3st(ep+1, s, t)/3.;
2425 }
2426 } else {
2427 assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2427, __extension__ __PRETTY_FUNCTION__));
2428 int iu[3]={0,s,t};
2429 int nu[3];
2430 freeidxM3(iu, nu);
2431 const double y11=Cay[nss(nu[0],nu[0])];
2432 const double y12=Cay[nss(nu[0],nu[1])];
2433 const double y13=Cay[nss(nu[0],nu[2])];
2434 const double y22=Cay[nss(nu[1],nu[1])];
2435 const double y23=Cay[nss(nu[1],nu[2])];
2436 const double y33=Cay[nss(nu[2],nu[2])];
2437 ivalue=-(+3*(y11*(y11+y12+y13)+y22*(y12+y22+y23)+y33*(y13+y23+y33))
2438 +2*(y12*(y12+y23)+y13*(y12+y13)+y23*(y13+y23))
2439 + (y11*(y22+y23)+y22*(y13+y33)+y33*(y11+y12))
2440 )/720.;
2441 }
2442 pI3D3st[ep][idx]=ivalue;
2443 }
2444 }
2445 fEval[E_I3D3st+ep]=true;
2446}
2447
2448/* --------------------------------------------------------
2449 * I4D4s box in D+8 dim
2450 * --------------------------------------------------------
2451 */
2452ncomplex Minor5::I4D4s(int ep, int s)
2453{
2454 if (ep==2) return 0;
2455 if (not fEval[E_I4D4s+ep]) {
2456 I4D4sEval(ep);
2457 }
2458 return pI4D4s[ep][s-1];
2459}
2460
2461void Minor5::I4D4sEval(int ep)
2462{
2463 for (int s=1; s<=smax; s++) {
2464 const double dss=M1(s, s);
2465 const double d0s0s=M2(0, s, 0, s);
2466 ncomplex ivalue=0;
2467
2468 if (ep==0) {
2469 ncomplex sum1=0;
2470 for (int t=1; t<=5; t++) {
2471 if (t==s) continue;
2472 sum1-=M2(s,t,s,0)*I3D3st(ep, s, t);
2473 }
2474 sum1+=d0s0s*I4D3s(ep, s);
2475 ivalue=sum1/(7*dss)+2.*I4D4s(ep+1, s)/7.;
2476
2477 const double x=dss/d0s0s;
2478 if (pmaxS4[s-1] <= deps3) {
2479 ncomplex sump;
2480 do {
2481 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2481, __extension__ __PRETTY_FUNCTION__));
2482
2483 double dsts0[6];
2484 sum1=0;
2485 for (int t=1; t<=5; t++) {
2486 if (t==s) continue;
2487 dsts0[t]=M2(s,t,s,0);
2488 sum1+=dsts0[t]*I3D4st(0, s, t);
2489 }
2490
2491 double xn=1;
2492 ncomplex dv,s21;
2493
2494 ncomplex sum[3];
2495 sum[0]=sump=sum1;
2496
2497#define stepI4D(n,a,b) \
2498 xn*=x; \
2499 dv=0; \
2500 for (int t=1; t<=5; t++) { \
2501 if (t==s) continue; \
2502 dv+=dsts0[t]*(a*I3D##n##st(0, s, t) - b*I3D##n##st(1, s, t)); \
2503 } \
2504 dv*=xn; \
2505 sum1+=dv;
2506
2507 stepI4D(5,9.,2.)
2508 if ( fabs(sum1.real()*teps)>=fabs(dv.real())
2509 && fabs(sum1.imag()*teps)>=fabs(dv.imag()))
2510 break;
2511 sum[1]=sump=sum1;
2512 s21=sum[1]-sum[0];
2513
2514 stepI4D(6,99.,40.)
2515 sump=sum1;
2516 stepWynn(0)sum[(2+0)%3]=sum1; { const ncomplex s2=sum[(2+0)%3]; const ncomplex
s1=sum[(1+0)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
2517 stepI4D(7,1287.,718.)
2518 stepWynn(1)sum[(2+1)%3]=sum1; { const ncomplex s2=sum[(2+1)%3]; const ncomplex
s1=sum[(1+1)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break;
2519// stepI4D(8,19305.,13344.)
2520// stepWynn(2)
2521// stepI4D(9,328185.,265458.)
2522// stepWynn(3)
2523// stepI4D(10,6235515.,5700072.)
2524// stepWynn(4)
2525// stepI4D(11,130945815.,132172542.)
2526// stepWynn(5)
2527#undef stepI4D
2528
2529 } while (0);
2530 ivalue=sump/d0s0s;
2531 }
2532 }
2533 else {
2534 assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2534, __extension__ __PRETTY_FUNCTION__));
2535 double sum1=0;
2536 for (int i=1; i<=5; i++) { if (i==s) continue;
2537 for (int j=i; j<=5; j++) { if (j==s) continue;
2538 for (int k=j; k<=5; k++) { if (k==s) continue;
2539 for (int l=k; l<=5; l++) { if (l==s) continue;
2540 sum1+=+Cay[nss(i,j)]*Cay[nss(k,l)]
2541 +Cay[nss(i,k)]*Cay[nss(j,l)]
2542 +Cay[nss(i,l)]*Cay[nss(j,k)];
2543 }
2544 }
2545 }
2546 }
2547 ivalue=sum1/5040.;
2548 }
2549 pI4D4s[ep][s-1]=ivalue;
2550 }
2551 fEval[E_I4D4s+ep]=true;
2552}
2553
2554
2555/* --------------------------------------------------------
2556 * I4D4si box in D+8 dim with a dot
2557 * --------------------------------------------------------
2558 */
2559ncomplex Minor5::I4D4si(int ep, int s, int i)
2560{
2561 if (s==i) return 0;
2562 if (ep==2) return 0;
2563 if (not fEval[E_I4D4si+ep]) {
2564 I4D4siEval(ep);
2565 }
2566 return pI4D4si[ep][i-1][s-1];
2567}
2568
2569void Minor5::I4D4siEval(int ep)
2570{
2571 for (int s=1; s<=smax; s++) {
2572 for (int i=1; i<=CIDX(DCay-2); i++) { if (s==i) continue;
2573 ncomplex ivalue=0;
2574
2575 if (ep == 0) {
2576 if (pmaxS4[s-1] <= deps3) {
2577 ncomplex sum1=0;
2578 for (int t=1; t<=5; t++) {
2579 if (t==s) continue;
2580 sum1+=M3(0, s, t, 0, s, i)*I3D3st(ep, s, t);
2581 }
2582 sum1+=M2(0, s, i, s)*(-7.*I4D4s(ep, s)+2.*I4D4s(ep+1, s));
2583 ivalue=sum1/M2(0, s, 0, s);
2584 } else {
2585 ncomplex sum1=0;
2586 for (int t=1; t<=5; t++) {
2587 if (t==s) continue;
2588 sum1+=M2(s,t,s,i)*I3D3st(ep, s, t);
2589 }
2590 sum1-=M2(s,0,s,i)*I4D3s(ep, s);
2591 sum1/=M1(s,s);
2592 ivalue=sum1;
2593 }
2594 } else {
2595 assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2595, __extension__ __PRETTY_FUNCTION__));
2596 double sum1=0;
2597 sum1+=Cay[nss(i,i)];
2598 for (int j=1; j<=5; j++) {
2599 if (j==s) continue;
2600 sum1+=Cay[ns(i,j)];
2601 for (int k=j; k<=5; k++) {
2602 if (k==s) continue;
2603 sum1+=Cay[nss(j,k)];
2604 }
2605 }
2606 ivalue=sum1/720.;
2607 }
2608 pI4D4si[ep][i-1][s-1]=ivalue;
2609 }
2610 }
2611 fEval[E_I4D4si+ep]=true;
2612}
2613
2614/* --------------------------------------------------------
2615 * I3D3sti triangle in D+6 dim with a dot
2616 * --------------------------------------------------------
2617 */
2618ncomplex Minor5::I3D3sti(int ep, int s, int t, int i)
2619{
2620 assert(s!=t && s!=i && t!=i)(static_cast <bool> (s!=t && s!=i && t!=
i) ? void (0) : __assert_fail ("s!=t && s!=i && t!=i"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 2620
, __extension__ __PRETTY_FUNCTION__));
2621 if (ep==2) return 0.;
2622 if (not fEval[E_I3D3sti+ep]) {
2623 I3D3stiEval(ep);
2624 }
2625 int idx = im2(s,t)-5;
2626 return pI3D3sti[ep][i-1][idx];
2627}
2628
2629void Minor5::I3D3stiEval(int ep)
2630{
2631 for (int i=1; i<=CIDX(DCay-2); i++) {
2632 for (int s=1; s<=smax; s++) { if (i==s) continue;
2633 for (int t=s+1; t<=5; t++) { if (i==t) continue;
2634 int idx = im2(s,t)-5;
2635 ncomplex ivalue=0;
2636
2637 if (ep==0) {
2638 if ( (pmaxT3[idx]==0 || (pmaxT3[idx] > epsir2 || pmaxU3[idx] > epsir2))
2639 && pmaxS3[idx] > ceps ) {
2640 // Variant with Gram3
2641 ncomplex sum1=0;
2642 for (int u=1; u<=5; u++) {
2643 if (u==t || u==s) continue;
2644 sum1+=M3(u,s,t,i,s,t)*I2D2stu(ep,s,t,u);
2645 }
2646 sum1-=M3(0,s,t,i,s,t)*I3D2st(ep,s,t);
2647 ivalue=sum1/M2(s,t,s,t);
2648 } else {
2649 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2649, __extension__ __PRETTY_FUNCTION__));
2650 int iu[3]={i-1,s-1,t-1};
2651 int tmp;
2652 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
2653 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
2654 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
2655 int nu[3];
2656 freeidxM3(iu, nu);
2657 int u=nu[0]+1;
2658 int v=nu[1]+1;
2659
2660 if ( pmaxT3[idx] <= epsir2 && pmaxU3[idx] <= epsir2 ) {
2661 // small G3 & C3
2662 int j=imax3[idx];
2663 ncomplex sum1=0;
2664 ncomplex const I3term=I3D2st(ep,s,t)+2./5.*I3D2st(ep+1,s,t);
2665 ncomplex const I2Uterm=I2D2stui(ep, s, t, u, i)+2./5.*I2D2stui(ep+1, s, t, u, i);
2666 ncomplex const I2Vterm=I2D2stui(ep, s, t, v, i)+2./5.*I2D2stui(ep+1, s, t, v, i);
2667
2668 if (j==i) { // j->i
2669 const double Dii=M4ii(u,v,i);
2670 const double Dui=M4ui(u,v,i);
2671 const double Dvi=M4vi(u,v,i);
2672 sum1+=+Dii*(I3term) // (i, i)
2673 +Dui*(I2Uterm) // (u, i)
2674 +Dvi*(I2Vterm); // (v, i)
2675 } else if (j==u) { // j->u
2676 const double Dui=M4ui(u,v,i);
2677 const double Duu=M4uu(u,v,i);
2678 const double Dvu=M4vu(u,v,i);
2679 sum1+=+Dui*(I3term) // (u, i)
2680 +Duu*(I2Uterm) // (u, u)
2681 +Dvu*(I2Vterm); // (v, u)
2682 } else { assert(j==v)(static_cast <bool> (j==v) ? void (0) : __assert_fail (
"j==v", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2682, __extension__ __PRETTY_FUNCTION__)); // j->v
2683 const double Dvi=M4vi(u,v,i);
2684 const double Dvv=M4vv(u,v,i);
2685 const double Dvu=M4vu(u,v,i);
2686 sum1+=+Dvi*(I3term) // (v, i)
2687 +Dvu*(I2Uterm) // (v, u)
2688 +Dvv*(I2Vterm); // (v, v)
2689 }
2690 ivalue=sum1/(5*M3(s,0,t,s,j,t));
2691 } else {
2692 // small G3
2693 const double Dii=M4ii(u,v,i);
2694 const double Dui=M4ui(u,v,i);
2695 const double Dvi=M4vi(u,v,i);
2696 ncomplex sum1=0;
2697 sum1+=Dii*I2D2stu(ep, s, t, i) // (i, i)
2698 +Dui*I2D2stu(ep, s, t, u) // (u, i)
2699 +Dvi*I2D2stu(ep, s, t, v); // (v, i)
2700 sum1+=M3(s,0,t,s,i,t)*(-6.*I3D3st(ep, s, t)+2.*I3D3st(ep+1, s, t));
2701 ivalue=sum1/M3(s,0,t,s,0,t);
2702 }
2703 }
2704 } else {
2705 assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2705, __extension__ __PRETTY_FUNCTION__));
2706 int iu[3]={0,s,t};
2707 int nu[3];
2708 freeidxM3(iu, nu);
2709 int j= nu[1]==i ? nu[0] : nu[1];
2710 int k= nu[2]==i ? nu[0] : nu[2];
2711 ivalue=-( 3*Cay[nss(i,i)]+2*(Cay[ns(i,j)]+Cay[ns(i,k)])
2712 +Cay[nss(j,j)]+Cay[ns(j,k)]+Cay[nss(k,k)]
2713 )/120.;
2714 }
2715 pI3D3sti[ep][i-1][idx]=ivalue;
2716 }
2717 }
2718 }
2719 fEval[E_I3D3sti+ep]=true;
2720}
2721
2722/* --------------------------------------------------------
2723 * I4D4sij box in D+8 dim with two dots
2724 * --------------------------------------------------------
2725 */
2726ncomplex Minor5::I4D4sij(int ep, int s, int i, int j)
2727{
2728 if (s==i || s==j) return 0;
2729 if (ep==1) return ( i==j ? 1./60. : 1./120. );
2730 else if (ep==2) return 0;
2731 if (not fEval[E_I4D4sij+ep]) {
2732 I4D4sijEval(ep);
2733 }
2734 return pI4D4sij[ep][is(i-1,j-1)][s-1];
2735}
2736
2737void Minor5::I4D4sijEval(int ep)
2738{
2739 for (int s=1; s<=smax; s++) {
2740 // symmetric in 'i,j'
2741 for (int i=1; i<=CIDX(DCay-2); i++) { if (s==i) continue;
2742 for (int j=i; j<=CIDX(DCay-2); j++) { if (s==j) continue;
2743 ncomplex ivalue=0;
2744
2745 if (pmaxS4[s-1] <= deps3) {
2746 ncomplex sum1=0;
2747 for (int t=1; t<=5; t++) {
2748 if (t==s || t==i) continue;
2749 sum1+=M3(s,0,t,s,0,j)*I3D3sti(ep, s, t, i);
2750 }
2751 sum1+=M3(s,0,i,s,0,j)*I4D3s(ep, s);
2752 sum1+=M2(s,0,s,j)*(-6.*I4D4si(ep, s, i)+2.*I4D4si(ep+1, s, i));
2753 ivalue=sum1/M2(0,s,0,s);
2754 } else {
2755 ncomplex sum1=0;
2756 for (int t=1; t<=5; t++) {
2757 if (t==s || t==i) continue;
2758 sum1+=M2(s,t,s,j)*I3D3sti(ep, s, t, i);
2759 }
2760 sum1+=M2(s,i,s,j)*I4D3s(ep, s);
2761 sum1-=M2(s,0,s,j)*I4D3si(ep, s, i);
2762 sum1/=M1(s,s);
2763 ivalue=sum1;
2764 }
2765 pI4D4sij[ep][iss(i-1,j-1)][s-1]=ivalue;
2766 }
2767 }
2768 }
2769 fEval[E_I4D4sij+ep]=true;
2770}
2771
2772/* --------------------------------------------------------
2773 * I2D2stui bubble in D+4 dim with a dot
2774 * --------------------------------------------------------
2775 */
2776ncomplex Minor5::I2D2stui(int ep, int s, int t, int u, int i)
2777{
2778 assert(s!=t && t!=u && u!=s && s!=i && t!=i && u!=i)(static_cast <bool> (s!=t && t!=u && u!=
s && s!=i && t!=i && u!=i) ? void (0)
: __assert_fail ("s!=t && t!=u && u!=s && s!=i && t!=i && u!=i"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 2778
, __extension__ __PRETTY_FUNCTION__));
2779 if (ep==2) return 0;
2780 if (not fEval[E_I2D2stui+ep]) {
2781 I2D2stuiEval(ep,1,4,5,2,3,kinem.p3());
2782 I2D2stuiEval(ep,1,3,5,2,4,kinem.s34());
2783 I2D2stuiEval(ep,1,3,4,2,5,kinem.s12());
2784 I2D2stuiEval(ep,1,4,5,3,2,kinem.p3());
2785 I2D2stuiEval(ep,1,2,5,3,4,kinem.p4());
2786 I2D2stuiEval(ep,1,2,4,3,5,kinem.s45());
2787 I2D2stuiEval(ep,1,3,5,4,2,kinem.s34());
2788 I2D2stuiEval(ep,1,2,5,4,3,kinem.p4());
2789 I2D2stuiEval(ep,1,2,3,4,5,kinem.p5());
2790#ifdef USE_ZERO_CHORD
2791 I2D2stuiEval(ep,1,3,4,5,2,kinem.s12());
2792 I2D2stuiEval(ep,1,2,4,5,3,kinem.s45());
2793 I2D2stuiEval(ep,1,2,3,5,4,kinem.p5());
2794#endif
2795
2796 if (smax==5) {
2797 I2D2stuiEval(ep,3,4,5,1,2,kinem.p2());
2798 I2D2stuiEval(ep,2,4,5,1,3,kinem.s23());
2799 I2D2stuiEval(ep,2,3,5,1,4,kinem.s15());
2800 I2D2stuiEval(ep,2,3,4,1,5,kinem.p1());
2801 I2D2stuiEval(ep,3,4,5,2,1,kinem.p2());
2802 I2D2stuiEval(ep,2,4,5,3,1,kinem.s23());
2803 I2D2stuiEval(ep,2,3,5,4,1,kinem.s15());
2804#ifdef USE_ZERO_CHORD
2805 I2D2stuiEval(ep,2,3,4,5,1,kinem.p1());
2806#endif
2807 }
2808
2809 fEval[E_I2D2stui+ep]=true;
2810 }
2811 int ip=15-s-t-u-i; // ip
2812 return pI2D2stui[ep][i-1][ip-1];
2813}
2814
2815void Minor5::I2D2stuiEval(int ep, int s, int t, int u, int i, int ip, double qsq)
2816{
2817 ncomplex sum1=0;
2818 if (ep==0) {
2819 const double dstustu=-2*qsq; /*M3(s,t,u,s,t,u);*/
2820 const double msq1=kinem.mass(i);
2821 const double msq2=kinem.mass(ip);
2822 const double s_cutoff=seps1*pmaxM2[im2(i,ip)-5];
2823
2824 if (fabs(dstustu) <= s_cutoff) {
2825 const double mm12=msq1-msq2;
2826 if (fabs(mm12) < meps) {
2827 sum1=0.5*ICache::getI1(ep, Kinem1(msq1));
2828 }
2829 else { assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2829, __extension__ __PRETTY_FUNCTION__));
2830 sum1=(4*msq1*msq1-5*(msq1+msq2)*msq2)/(36*mm12)
2831 + ( msq1*(2*msq1 - 3*msq2)*ICache::getI1(ep, Kinem1(msq1))
2832 + msq2*msq2*ICache::getI1(ep, Kinem1(msq2))
2833 )/(6*mm12*mm12);
2834 }
2835 }
2836 else {
2837 sum1+=+(Cay[nss(ip,ip)]-Cay[ns(i,ip)])*I2Dstu(ep,s,t,u);
2838 sum1+=-0.5*msq1*(ICache::getI1(ep, Kinem1(msq1))+0.5*msq1);
2839 sum1+=+0.5*msq2*(ICache::getI1(ep, Kinem1(msq2))+0.5*msq2);
2840 sum1/=dstustu;
2841 }
2842 }
2843 else { assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2843, __extension__ __PRETTY_FUNCTION__));
2844 if ( fabs(qsq) < meps
2845 && fabs(kinem.mass(i)) < meps
2846 && fabs(kinem.mass(ip)) < meps ) {
2847 sum1=0;
2848 }
2849 else {
2850 sum1=(3*Cay[nss(i,i)] + 2*Cay[ns(i,ip)] + Cay[nss(ip,ip)])/24.;
2851 }
2852 }
2853 pI2D2stui[ep][i-1][ip-1]=sum1;
2854}
2855
2856
2857/* --------------------------------------------------------
2858 * I3D3stij triangle in D+6 dim with two dots
2859 * --------------------------------------------------------
2860 */
2861ncomplex Minor5::I3D3stij(int ep, int s, int t, int i, int j)
2862{
2863 assert(s!=t && s!=i && s!=j && t!=i && t!=j)(static_cast <bool> (s!=t && s!=i && s!=
j && t!=i && t!=j) ? void (0) : __assert_fail
("s!=t && s!=i && s!=j && t!=i && t!=j"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 2863
, __extension__ __PRETTY_FUNCTION__));
2864 if (ep==1) return ( i==j ? -1./12. : -1./24. ); // -1/12 == -2/24
2865 else if (ep==2) return 0;
2866 if (not fEval[E_I3D3stij+ep]) {
2867 I3D3stijEval(ep);
2868 }
2869 int idx = im2(s,t)-5;
2870 return pI3D3stij[ep][is(i-1,j-1)][idx];
2871}
2872
2873void Minor5::I3D3stijEval(int ep)
2874{
2875 for (int s=1; s<=smax; s++) {
2876 for (int t=s+1; t<=5; t++) {
2877 int idx = im2(s,t)-5;
2878 const double dstst=M2(s,t,s,t);
2879 // symmetric in 'i,j'
2880 for (int i=1; i<=CIDX(DCay-2); i++) { if (i==s || i==t) continue;
2881 for (int j=i; j<=CIDX(DCay-2); j++) { if (j==s || j==t) continue;
2882 ncomplex ivalue=0;
2883
2884 if ( (pmaxT3[idx]==0 || (pmaxT3[idx] > epsir2 || pmaxU3[idx] > epsir2))
2885 && pmaxS3[idx] > ceps ) {
2886 // Variant with Gram3
2887 ncomplex sum1=0;
2888 for (int u=1; u<=5; u++) {
2889 if (u==t || u==s || u==i) continue;
2890 sum1+=M3(s,u,t,s,j,t)*I2D2stui(ep, s, t, u, i);
2891 }
2892 sum1+=-M3(s,0,t,s,j,t)*I3D2sti(ep, s, t, i)+M3(s,i,t,s,j,t)*I3D2st(ep, s, t);
2893 ivalue=sum1/dstst;
2894 } else {
2895 assert(ep==0)(static_cast <bool> (ep==0) ? void (0) : __assert_fail (
"ep==0", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 2895, __extension__ __PRETTY_FUNCTION__));
2896 int iu[3]={j-1,s-1,t-1};
2897 int tmp;
2898 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
2899 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
2900 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
2901 int nu[3];
2902 freeidxM3(iu, nu);
2903 int u=nu[0]+1;
2904 int v=nu[1]+1;
2905
2906 const double Djj=M4ii(u,v,j);
2907 const double Duj=M4ui(u,v,j);
2908 const double Dvj=M4vi(u,v,j);
2909 if ( pmaxT3[idx] <= epsir2 && pmaxU3[idx] <= epsir2 ) {
2910 // small G3 & C3
2911 int k=imax3[idx];
2912 ncomplex sum1=0;
2913 if (i==j) {
2914 if (k==j) {
2915 sum1+=2*Djj*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2916 +Duj*(I2D2stuij(ep,s,t,u,j,j)+0.5*I2D2stuij(ep+1,s,t,u,j,j))
2917 +Dvj*(I2D2stuij(ep,s,t,v,j,j)+0.5*I2D2stuij(ep+1,s,t,v,j,j));
2918 } else if (k==u) {
2919 const double Duu=M4uu(u,v,j);
2920 const double Duv=M4vu(u,v,j);
2921 sum1+=2*Duj*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2922 +Duu*(I2D2stuij(ep,s,t,u,j,j)+0.5*I2D2stuij(ep+1,s,t,u,j,j))
2923 +Duv*(I2D2stuij(ep,s,t,v,j,j)+0.5*I2D2stuij(ep+1,s,t,v,j,j));
2924 } else { // k==v
2925 const double Dvv=M4vv(u,v,j);
2926 const double Dvu=M4vu(u,v,j);
2927 sum1+=2*Dvj*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2928 +Dvu*(I2D2stuij(ep,s,t,u,j,j)+0.5*I2D2stuij(ep+1,s,t,u,j,j))
2929 +Dvv*(I2D2stuij(ep,s,t,v,j,j)+0.5*I2D2stuij(ep+1,s,t,v,j,j));
2930 }
2931 } else if (k==j) {
2932 sum1+=Djj*(I3D2sti(ep,s,t,i)+0.5*I3D2sti(ep+1,s,t,i));
2933 if (i==u) {
2934 sum1+=Duj*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2935 +Dvj*(I2D2stuij(ep,s,t,v,i,j)+0.5*I2D2stuij(ep+1,s,t,v,i,j));
2936 } else { // i==v
2937 sum1+=Dvj*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2938 +Duj*(I2D2stuij(ep,s,t,u,i,j)+0.5*I2D2stuij(ep+1,s,t,u,i,j));
2939 }
2940 } else if (k==i) {
2941 if (k==u) {
2942 const double Duu=M4uu(u,v,j);
2943 const double Duv=M4vu(u,v,j);
2944 sum1+=Duj*(I3D2sti(ep,s,t,i)+0.5*I3D2sti(ep+1,s,t,i))
2945 +Duu*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2946 +Duv*(I2D2stuij(ep,s,t,v,i,j)+0.5*I2D2stuij(ep+1,s,t,v,i,j));
2947 } else { // k==v
2948 const double Dvv=M4vv(u,v,j);
2949 const double Dvu=M4vu(u,v,j);
2950 sum1+=Dvj*(I3D2sti(ep,s,t,i)+0.5*I3D2sti(ep+1,s,t,i))
2951 +Dvv*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2952 +Dvu*(I2D2stuij(ep,s,t,u,i,j)+0.5*I2D2stuij(ep+1,s,t,u,i,j));
2953 }
2954 } else {
2955 if (k==u) { // i==v
2956 const double Duu=M4uu(u,v,j);
2957 const double Duv=M4vu(u,v,j);
2958 sum1+=Duj*(I3D2sti(ep,s,t,i)+0.5*I3D2sti(ep+1,s,t,i))
2959 +Duv*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2960 +Duu*(I2D2stuij(ep,s,t,u,i,j)+0.5*I2D2stuij(ep+1,s,t,u,i,j));
2961 } else { // k==v, i==u
2962 const double Dvv=M4vv(u,v,j);
2963 const double Dvu=M4vu(u,v,j);
2964 sum1+=Dvj*(I3D2sti(ep,s,t,i)+0.5*I3D2sti(ep+1,s,t,i))
2965 +Dvu*(I3D2sti(ep,s,t,j)+0.5*I3D2sti(ep+1,s,t,j))
2966 +Dvv*(I2D2stuij(ep,s,t,v,i,j)+0.5*I2D2stuij(ep+1,s,t,v,i,j));
2967 }
2968 }
2969 ivalue=sum1/(4*M3(s,0,t,s,k,t));
2970 } else {
2971 // small G3
2972 ncomplex sum1=0;
2973 if (j==i) {
2974 sum1+=+Djj*I3D2st(ep,s,t)
2975 +Duj*I2D2stui(ep, s, t, u, i)
2976 +Dvj*I2D2stui(ep, s, t, v, i);
2977 } else {
2978 sum1+=Djj*I2D2stui(ep, s, t, j, i);
2979 if (i==u) {
2980 sum1+=+Duj*I3D2st(ep,s,t)
2981 +Dvj*I2D2stui(ep, s, t, v, i);
2982 } else {
2983 sum1+=+Dvj*I3D2st(ep,s,t)
2984 +Duj*I2D2stui(ep, s, t, u, i);
2985 }
2986 }
2987 sum1+=M3(s,0,t,s,j,t)*(-5.*I3D3sti(ep, s, t, i)+2.*I3D3sti(ep+1, s, t, i));
2988 ivalue=sum1/M3(s,0,t,s,0,t);
2989 }
2990 }
2991 pI3D3stij[ep][iss(i-1,j-1)][idx]=ivalue;
2992 }
2993 }
2994 }
2995 }
2996 fEval[E_I3D3stij+ep]=true;
2997}
2998
2999/* --------------------------------------------------------
3000 * I4D4sijk box in D+8 dim with three dots
3001 * --------------------------------------------------------
3002 */
3003ncomplex Minor5::I4D4sijk(int ep, int s, int i, int j, int k)
3004{
3005 if (s==i || s==j || s==k) return 0;
3006 if (ep==2) return 0; // I4D4sijk finite
3007 if (not fEval[E_I4D4sijk+ep]) {
3008 I4D4sijkEval(ep);
3009 }
3010 return pI4D4sijk[ep][is(i-1,j-1,k-1)][s-1];
3011}
3012
3013void Minor5::I4D4sijkEval(int ep)
3014{
3015 for (int s=1; s<=smax; s++) {
3016 // symmetric in 'i,j,k'
3017 for (int i=1; i<=CIDX(DCay-2); i++) { if (i==s) continue;
3018 for (int j=i; j<=CIDX(DCay-2); j++) { if (j==s) continue;
3019 for (int k=j; k<=CIDX(DCay-2); k++) { if (k==s) continue;
3020 ncomplex ivalue=0;
3021
3022 if (pmaxS4[s-1] <= deps3) {
3023 ncomplex sum1=0;
3024 for (int t=1; t<=5; t++) {
3025 if (s==t || t==i || t==j) continue;
3026 sum1+=M3(s,0,t,s,0,k)*I3D3stij(ep,s,t,i,j);
3027 }
3028 sum1+=+M3(s,0,i,s,0,k)*I4D3si(ep, s, j)
3029 +M3(s,0,j,s,0,k)*I4D3si(ep, s, i);
3030
3031 sum1+=M2(s,0,s,k)*(-5.*I4D4sij(ep, s, i, j)+2.*I4D4sij(ep+1, s, i, j));
3032
3033 ivalue=sum1/M2(s,0,s,0);
3034 } else {
3035 ncomplex sum1=0;
3036 for (int t=1; t<=5; t++) {
3037 if (t==s || t==i || t==j) continue;
3038 sum1+=M2(s,t,s,k)*I3D3stij(ep,s,t,i,j);
3039 }
3040 sum1-=M2(s,0,s,k)*I4D3sij(ep,s,i,j);
3041 sum1+=M2(s,i,s,k)*I4D3si(ep,s,j)+M2(s,j,s,k)*I4D3si(ep,s,i);
3042 ivalue=sum1/M1(s,s);
3043 }
3044 pI4D4sijk[ep][iss(i-1,j-1,k-1)][s-1]=ivalue;
3045 }
3046 }
3047 }
3048 }
3049 fEval[E_I4D4sijk+ep]=true;
3050}
3051
3052/* --------------------------------------------------------
3053 * I2D2stuij bubble in D+4 dim with two dots
3054 * --------------------------------------------------------
3055 */
3056ncomplex Minor5::I2D2stuij(int ep, int s, int t, int u, int i, int j)
3057{
3058 assert(s!=t && t!=u && u!=s && s!=i && t!=i && u!=i && s!=j && t!=j && u!=j)(static_cast <bool> (s!=t && t!=u && u!=
s && s!=i && t!=i && u!=i && s
!=j && t!=j && u!=j) ? void (0) : __assert_fail
("s!=t && t!=u && u!=s && s!=i && t!=i && u!=i && s!=j && t!=j && u!=j"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 3058
, __extension__ __PRETTY_FUNCTION__));
3059 if (ep==2) return 0;
3060 if (not fEval[E_I2D2stuij+ep]) {
3061 I2D2stuijEval(ep,1,2,3,4,5,kinem.p5());
3062 I2D2stuijEval(ep,1,2,4,3,5,kinem.s45());
3063 I2D2stuijEval(ep,1,2,5,3,4,kinem.p4());
3064 I2D2stuijEval(ep,1,2,5,4,3,kinem.p4());
3065
3066 I2D2stuijEval(ep,1,3,4,2,5,kinem.s12());
3067 I2D2stuijEval(ep,1,3,5,2,4,kinem.s34());
3068 I2D2stuijEval(ep,1,3,5,4,2,kinem.s34());
3069
3070 I2D2stuijEval(ep,1,4,5,2,3,kinem.p3());
3071 I2D2stuijEval(ep,1,4,5,3,2,kinem.p3());
3072
3073#ifdef USE_ZERO_CHORD
3074 I2D2stuijEval(ep,1,2,3,5,4,kinem.p5());
3075 I2D2stuijEval(ep,1,2,4,5,3,kinem.s45());
3076 I2D2stuijEval(ep,1,3,4,5,2,kinem.s12());
3077#endif
3078 if (smax==5) {
3079 I2D2stuijEval(ep,2,3,4,1,5,kinem.p1());
3080 I2D2stuijEval(ep,2,3,5,1,4,kinem.s15());
3081 I2D2stuijEval(ep,2,3,5,4,1,kinem.s15());
3082 I2D2stuijEval(ep,2,4,5,1,3,kinem.s23());
3083 I2D2stuijEval(ep,2,4,5,3,1,kinem.s23());
3084 I2D2stuijEval(ep,3,4,5,1,2,kinem.p2());
3085 I2D2stuijEval(ep,3,4,5,2,1,kinem.p2());
3086#ifdef USE_ZERO_CHORD
3087 I2D2stuijEval(ep,2,3,4,5,1,kinem.p1());
3088#endif
3089 }
3090
3091 fEval[E_I2D2stuij+ep]=true;
3092 }
3093 int ip=15-s-t-u-i; // ip
3094 return pI2D2stuij[ep][i-1][ip-1][i==j ? 0 : 1];
3095}
3096
3097void Minor5::I2D2stuijEval(int ep, int s, int t, int u, int i, int ip, double qsq)
3098{
3099 ncomplex sum0=0;
3100 ncomplex sum1=0;
3101 if (ep==0) {
3102 const double dstustu=-2*qsq; /*M3(s,t,u,s,t,u);*/
3103 const double msq1=kinem.mass(i);
3104 const double msq2=kinem.mass(ip);
3105 const double s_cutoff=seps2*pmaxM2[im2(i,ip)-5];
3106
3107 if (fabs(dstustu) <= s_cutoff) {
3108 const double mm12=msq1-msq2;
3109 if (fabs(mm12) < meps) {
3110 if (msq1 > meps) {
3111 sum1=( ICache::getI1(ep, Kinem1(msq1)) - msq1 )/(6*msq1);
3112 } else {
3113 sum1=0;
3114 }
3115 sum0=2.*sum1;
3116 }
3117 else {
3118 sum0=2.*( (-4*msq1*msq1 + 5*msq1*msq2 + 5*msq2*msq2)/6.
3119 + ( (msq1*msq1 - 3*msq1*msq2 + 3*msq2*msq2)*ICache::getI1(ep, Kinem1(msq1))
3120 - msq2*msq2*ICache::getI1(ep, Kinem1(msq2))
3121 )/mm12
3122 )/(6.*mm12*mm12);
3123 sum1=(-(msq1*msq1 + 10*msq1*msq2 + msq2*msq2)/6. +
3124 ( msq1*(msq1 - 3*msq2)*ICache::getI1(ep, Kinem1(msq1))
3125 + (3*msq1 - msq2)*msq2*ICache::getI1(ep, Kinem1(msq2))
3126 )/mm12
3127 )/(6.*mm12*mm12);
3128 }
3129 }
3130 else {
3131 sum0+=+(Cay[nss(ip,ip)]-Cay[ns(i,ip)])*I2Dstui(ep,s,t,u,i);
3132 sum0+=-ICache::getI1(ep, Kinem1(msq1));
3133 sum0+=-I2Dstu(ep,s,t,u);
3134 sum0/=dstustu;
3135
3136 sum1+=(Cay[nss(i ,i )]-Cay[ns(i,ip)])*I2Dstui(ep,s,t,u,i);
3137 sum1+=ICache::getI1(ep, Kinem1(msq1));
3138 /* Symmetrization is not needed */
3139// sum1+=(Cay[nss(ip,ip)]-Cay[ns(ip,i)])*I2Dstui(ep,s,t,u,ip);
3140// sum1+=ICache::getI1(ep, Kinem1(msq2));
3141// sum1/=2.0;
3142 sum1+=I2Dstu(ep,s,t,u);
3143 sum1/=dstustu;
3144 }
3145 }
3146 else { assert(ep==1)(static_cast <bool> (ep==1) ? void (0) : __assert_fail (
"ep==1", "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp"
, 3146, __extension__ __PRETTY_FUNCTION__));
3147 if ( fabs(qsq) < meps
3148 && fabs(kinem.mass(i)) < meps
3149 && fabs(kinem.mass(ip)) < meps ) {
3150 sum0=0;
3151 sum1=0;
3152 }
3153 else {
3154 sum0=2./6.;
3155 sum1=1./6.;
3156 }
3157 }
3158 pI2D2stuij[ep][i-1][ip-1][0]=sum0;
3159 pI2D2stuij[ep][i-1][ip-1][1]=sum1;
3160}
3161
3162/* --------------------------------------------------------
3163 * I3D3stijk triangle in D+6 dim with three dots
3164 * --------------------------------------------------------
3165 */
3166ncomplex Minor5::I3D3stijk(int ep, int s, int t, int i, int j, int k) // IR-div
3167{
3168 assert(s!=t && s!=i && s!=j && s!=k && t!=i && t!=j && t!=k)(static_cast <bool> (s!=t && s!=i && s!=
j && s!=k && t!=i && t!=j && t
!=k) ? void (0) : __assert_fail ("s!=t && s!=i && s!=j && s!=k && t!=i && t!=j && t!=k"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.cpp", 3168
, __extension__ __PRETTY_FUNCTION__));
3169 if (not fEval[E_I3D3stijk+ep]) {
3170 I3D3stijkEval(ep);
3171 }
3172 int idx = im2(s,t)-5;
3173 return pI3D3stijk[ep][is(i-1,j-1,k-1)][idx];
3174}
3175
3176void Minor5::I3D3stijkEval(int ep)
3177{
3178 for (int s=1; s<=smax; s++) {
3179 for (int t=s+1; t<=5; t++) {
3180 int idx = im2(s,t)-5;
3181
3182 const double ds0ts0t=M3(s,0,t,s,0,t);
3183 if (ep!=0 && fabs(ds0ts0t) > m3eps) { // if ds0ts0t!=0 I3D3stijk is finite
3184 for (int ijk=iss(1-1,1-1,1-1); ijk<=iss(CIDX(DCay-2)-1,CIDX(DCay-2)-1,CIDX(DCay-2)-1); ijk++) {
3185 pI3D3stijk[ep][ijk][idx]=0;
3186 }
3187 continue;
3188 }
3189
3190 const double dstst=M2(s,t,s,t);
3191 for (int i=1; i<=CIDX(DCay-2); i++) { if (i==s || i==t) continue;
3192 for (int j=i; j<=CIDX(DCay-2); j++) { if (j==s || j==t) continue;
3193 for (int k=j; k<=CIDX(DCay-2); k++) { if (k==s || k==t) continue;
3194 ncomplex ivalue=0;
3195
3196 if (pmaxS3[idx] > ceps || ep!=0) {
3197 // Variant with Gram3
3198 ncomplex sum1=0;
3199 for (int u=1; u<=5; u++) {
3200 if (u==s || u==t || u==i || u==j) continue;
3201 sum1+=M3(s,u,t,s,k,t)*I2D2stuij(ep, s, t, u, i, j);
3202 }
3203 sum1-=M3(s,0,t,s,k,t)*I3D2stij(ep, s, t, i, j);
3204 sum1+=M3(s,i,t,s,k,t)*I3D2sti(ep, s, t, j)+M3(s,j,t,s,k,t)*I3D2sti(ep, s, t, i);
3205 ivalue=sum1/dstst;
3206 } else {
3207 ncomplex sum1=0;
3208 int iu[3]={k-1,s-1,t-1};
3209 int tmp;
3210 tswap(iu[0],iu[2],tmp)if (iu[0] > iu[2]) { tmp=iu[2]; iu[2]=iu[0]; iu[0]=tmp; };
3211 tswap(iu[1],iu[2],tmp)if (iu[1] > iu[2]) { tmp=iu[2]; iu[2]=iu[1]; iu[1]=tmp; };
3212 tswap(iu[0],iu[1],tmp)if (iu[0] > iu[1]) { tmp=iu[1]; iu[1]=iu[0]; iu[0]=tmp; };
3213 int nu[3];
3214 freeidxM3(iu, nu);
3215 int u=nu[0]+1;
3216 int v=nu[1]+1;
3217 const double Dkk=M4ii(u,v,k);
3218 const double Duk=M4ui(u,v,k);
3219 const double Dvk=M4vi(u,v,k);
3220 if ( fabs(ds0ts0t) > 0. ) {
3221 if (j==i) {
3222 if (j==k) {
3223 sum1+=+2*Dkk*I3D2sti(ep,s,t,j)
3224 +Duk*I2D2stuij(ep, s, t, u, j, j)
3225 +Dvk*I2D2stuij(ep, s, t, v, j, j);
3226 } else {
3227 sum1+=Dkk*I2D2stuij(ep, s, t, k, j, j);
3228 if (j==u) {
3229 sum1+=+2*Duk*I3D2sti(ep,s,t,j)
3230 +Dvk*I2D2stuij(ep, s, t, v, j, j);
3231 } else {
3232 sum1+=+2*Dvk*I3D2sti(ep,s,t,j)
3233 +Duk*I2D2stuij(ep, s, t, u, j, j);
3234 }
3235 }
3236 } else {
3237 if (j==k) {
3238 sum1+=+Dkk*I3D2sti(ep,s,t,i);
3239 if (i==u) {
3240 sum1+=+Duk*I3D2sti(ep,s,t,j)
3241 +Dvk*I2D2stuij(ep, s, t, v, i, j);
3242 } else {
3243 sum1+=+Dvk*I3D2sti(ep,s,t,j)
3244 +Duk*I2D2stuij(ep, s, t, u, i, j);
3245 }
3246 } else {
3247 sum1+=+Duk*I3D2sti(ep,s,t,v)
3248 +Dvk*I3D2sti(ep,s,t,u)
3249 +Dkk*I2D2stuij(ep, s, t, k, i, j);
3250 }
3251 }
3252 if (ep<2)
3253 sum1+=M3(s,0,t,s,k,t)*(-4.*I3D3stij(ep, s, t, i, j)+2.*I3D3stij(ep+1, s, t, i, j));
3254 else
3255 sum1+=M3(s,0,t,s,k,t)*(-4.*I3D3stij(ep, s, t, i, j));
3256 ivalue=sum1/ds0ts0t;
3257 } else {
3258 ivalue=std::numeric_limits<double>::quiet_NaN();
3259 // TODO add and check, needs I2D2stuijk
3260 }
3261 }
3262 pI3D3stijk[ep][iss(i-1,j-1,k-1)][idx]=ivalue;
3263 }
3264 }
3265 }
3266 }
3267 }
3268 fEval[E_I3D3stijk+ep]=true;
3269}
3270
3271/* --------------------------------------------------------
3272 * I4D4sijkl box in D+8 dim with four dots
3273 * --------------------------------------------------------
3274 */
3275ncomplex Minor5::I4D4sijkl(int ep, int s, int i, int j, int k, int l) // IR-div
3276{
3277 if (s==i || s==j || s==k || s==l) return 0;
3278 if (not fEval[E_I4D4sijkl+ep]) {
3279 I4D4sijklEval(ep);
3280 }
3281 return pI4D4sijkl[ep][is(i-1,j-1,k-1,l-1)][s-1];
3282}
3283
3284void Minor5::I4D4sijklEval(int ep)
3285{
3286 for (int s=1; s<=smax; s++) {
3287 // symmetric in 'i,j,k,l'
3288 for (int i=1; i<=CIDX(DCay-2); i++) { if (s==i) continue;
3289 for (int j=i; j<=CIDX(DCay-2); j++) { if (s==j) continue;
3290 for (int k=j; k<=CIDX(DCay-2); k++) { if (s==k) continue;
3291 for (int l=k; l<=CIDX(DCay-2); l++) { if (s==l) continue;
3292 ncomplex ivalue=0;
3293
3294 if (pmaxS4[s-1] <= deps3) {
3295 ncomplex sum1=0;
3296 for (int t=1; t<=5; t++) {
3297 if (s==t || t==i || t==j || t==k) continue;
3298 sum1+=M3(s,0,t,s,0,l)*I3D3stijk(ep,s,t,i,j,k);
3299 }
3300 sum1+=+M3(s,0,i,s,0,l)*I4D3sij(ep, s, j, k)
3301 +M3(s,0,j,s,0,l)*I4D3sij(ep, s, i, k)
3302 +M3(s,0,k,s,0,l)*I4D3sij(ep, s, i, j);
3303 if (ep<2) {
3304 sum1+=M2(s,0,s,l)*(-4.*I4D4sijk(ep, s, i, j, k)+2.*I4D4sijk(ep+1, s, i, j, k));
3305 }
3306 else { // ep==2
3307 sum1+=M2(s,0,s,l)*(-4.*I4D4sijk(ep, s, i, j, k));
3308 }
3309 ivalue=sum1/M2(s,0,s,0);
3310 } else {
3311 ncomplex sum1=0;
3312 for (int t=1; t<=5; t++) {
3313 if (t==s || t==i || t==j || t==k) continue;
3314 sum1+=M2(s,t,s,l)*I3D3stijk(ep,s,t,i,j,k);
3315 }
3316 sum1-=M2(s,0,s,l)*I4D3sijk(ep,s,i,j,k);
3317 sum1+=M2(s,i,s,l)*I4D3sij(ep,s,j,k)
3318 +M2(s,j,s,l)*I4D3sij(ep,s,i,k)
3319 +M2(s,k,s,l)*I4D3sij(ep,s,i,j);
3320 ivalue=sum1/M1(s,s);
3321 }
3322 pI4D4sijkl[ep][iss(i-1,j-1,k-1,l-1)][s-1]=ivalue;
3323 }
3324 }
3325 }
3326 }
3327 }
3328 fEval[E_I4D4sijkl+ep]=true;
3329}
←

generators/phokhara/phokhara/eemmg-lib/src/minor.h

→
1/*
2 * minor.h - signed minors classes
3 *
4 * this file is part of PJFry library
5 * Copyright 2011 Valery Yundin
6 */
7
8#ifndef QUL_MINOR_H
9#define QUL_MINOR_H
10
11#include "common.h"
12#include "kinem.h"
13#include "pointer.h"
14#include <bitset>
15
16// one step of Wynn epsilon improvement
17#define stepWynn(n)sum[(2+n)%3]=sum1; { const ncomplex s2=sum[(2+n)%3]; const ncomplex
s1=sum[(1+n)%3]; const ncomplex s10=s21; s21=s2-s1; if ( s21
==s10 || ( fabs(s2.real()*heps)>=fabs(s21.real()) &&
fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) break; dv=sump
; sump=s1+1./(1./s21-1./s10); } if ( fabs(sump.real()*teps)>=
fabs(sump.real()-dv.real()) && fabs(sump.imag()*teps)
>=fabs(sump.imag()-dv.imag()) ) break; \
18 sum[(2+n)%3]=sum1; \
19 { \
20 const ncomplex s2=sum[(2+n)%3]; \
21 const ncomplex s1=sum[(1+n)%3]; \
22 const ncomplex s10=s21; \
23 s21=s2-s1; \
24 if ( s21==s10 \
25 || ( fabs(s2.real()*heps)>=fabs(s21.real()) \
26 && fabs(s2.imag()*heps)>=fabs(s21.imag()) ) ) \
27 break; \
28 dv=sump; \
29 sump=s1+1./(1./s21-1./s10); \
30 } \
31 if ( fabs(sump.real()*teps)>=fabs(sump.real()-dv.real()) \
32 && fabs(sump.imag()*teps)>=fabs(sump.imag()-dv.imag()) ) \
33 break;
34// -------------
35
36#define tswap(x,y,t)if (x > y) { t=y; y=x; x=t; } \
37 if (x > y) { \
38 t=y; \
39 y=x; \
40 x=t; \
41 }
42
43#ifndef USE_ZERO_CHORD
44/* if zero chord is not used - calculate only 4 form factors for 5-point */
45# define CIDX(DCay-2) (DCay-2)
46#else
47/* else calculate all (including coefficient of 0'th chord) */
48# define CIDX(DCay-2) (DCay-1)
49#endif
50
51class MinorBase : public SRefCnt {
52public:
53 MinorBase() {}
54
55#ifdef USE_GOLEM_MODE
56#define EVALUNDEF return std::numeric_limits<double>::signaling_NaN();
57 virtual ncomplex A(int ep) { EVALUNDEF }
58 virtual ncomplex A(int ep, int i) { EVALUNDEF }
59 virtual ncomplex A(int ep, int i, int j) { EVALUNDEF }
60 virtual ncomplex A(int ep, int i, int j, int k) { EVALUNDEF }
61 virtual ncomplex A(int ep, int i, int j, int k, int l) { EVALUNDEF }
62 virtual ncomplex A(int ep, int i, int j, int k, int l, int m) { EVALUNDEF }
63 virtual ncomplex B(int ep) { EVALUNDEF }
64 virtual ncomplex B(int ep, int i) { EVALUNDEF; }
65 virtual ncomplex B(int ep, int i, int j) { EVALUNDEF }
66 virtual ncomplex B(int ep, int i, int j, int k) { EVALUNDEF }
67 virtual ncomplex C(int ep) { EVALUNDEF }
68 virtual ncomplex C(int ep, int i) { EVALUNDEF }
69#undef EVALUNDEF
70#endif /* USE_GOLEM_MODE */
71
72 // Symmetric index - start from 1
73 inline static int ns(int i, int j) CONST__attribute__ ((const)) {
74 return (i <= j ? (i - 1) + ((j - 1) * j) / 2 : (j - 1) + ((i - 1) * i) / 2);
75 }
76
77 inline static int nss(int i, int j) CONST__attribute__ ((const)) { // ordered
78 return (i - 1) + ((j - 1) * j) / 2;
79 }
80
81 // Symmetric index - generic
82 inline static int is(int i, int j) CONST__attribute__ ((const)) {
83 return (i <= j ? i + j * (j + 1) / 2 : j + i * (i + 1) / 2);
84 }
85
86 inline static int is(int i, int j, int k) CONST__attribute__ ((const)) {
87 if (i <= j)
88 {
89 return (j <= k ? i + ti2[j] + ti3[k] : is(i, k) + ti3[j]);
90 } else {
91 return (i > k ? is(j, k) + ti3[i] : j + ti2[i] + ti3[k]);
92 }
93 }
94
95 inline static int is(int i, int j, int k, int l) CONST__attribute__ ((const)) {
96 if (i <= j)
97 {
98 if (j <= k) {
99 return (k <= l ? i + ti2[j] + ti3[k] + ti4[l]
100 : is(i, j, l) + ti4[k]);
101 } else {
102 return (j > l ? is(i, k, l) + ti4[j]
103 : is(i, k) + ti3[j] + ti4[l]);
104 }
105 } else {
106 if (i > k)
107 {
108 return (i > l ? is(j, k, l) + ti4[i]
109 : is(j, k) + ti3[i] + ti4[l]);
110 } else {
111 return (k <= l ? j + ti2[i] + ti3[k] + ti4[l]
112 : is(i, j, l) + ti4[k]);
113 }
114 }
115 }
116
117 inline static int iss(int i, int j) CONST__attribute__ ((const)) { // ordered
118 assert(i <= j)(static_cast <bool> (i <= j) ? void (0) : __assert_fail
("i <= j", "generators/phokhara/phokhara/eemmg-lib/src/minor.h"
, 118, __extension__ __PRETTY_FUNCTION__));
119 return i + j * (j + 1) / 2;
120 }
121
122 inline static int iss(int i, int j, int k) CONST__attribute__ ((const)) { // ordered
123 assert(i <= j&& j <= k)(static_cast <bool> (i <= j&& j <= k) ? void
(0) : __assert_fail ("i <= j&& j <= k", "generators/phokhara/phokhara/eemmg-lib/src/minor.h"
, 123, __extension__ __PRETTY_FUNCTION__));
124 return i + ti2[j] + ti3[k];
125 }
126
127 inline static int iss(int i, int j, int k, int l) CONST__attribute__ ((const)) { // ordered
128 assert(i <= j&& j <= k&& k <= l)(static_cast <bool> (i <= j&& j <= k&&
k <= l) ? void (0) : __assert_fail ("i <= j&& j <= k&& k <= l"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.h", 128, __extension__
__PRETTY_FUNCTION__));
129 return i + ti2[j] + ti3[k] + ti4[l];
130 }
131
132 inline static int iss(int i, int j, int k, int l, int m) CONST__attribute__ ((const)) { // ordered
133 assert(i <= j&& j <= k&& k <= l&& l <= m)(static_cast <bool> (i <= j&& j <= k&&
k <= l&& l <= m) ? void (0) : __assert_fail ("i <= j&& j <= k&& k <= l&& l <= m"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.h", 133, __extension__
__PRETTY_FUNCTION__));
134 return i + ti2[j] + ti3[k] + ti4[l] + ti5[m];
135 }
136
137 inline static double getmeps()
138 {
139 return meps;
140 }
141
142 // Utility functions
143 static int im3(int i, int j, int k) CONST__attribute__ ((const)); // Antisymmetric index for "3 out of 6" minors
144 static int im2(int i, int j) CONST__attribute__ ((const)); // Antisymmetric index for "2 out of 6" minors
145 static int signM3ud(int i, int j, int k, int l, int m, int n) CONST__attribute__ ((const)); // Signature[{i,j,k}]*Signature[{l,m,n}]
146 static int signM2ud(int i, int j, int l, int m) CONST__attribute__ ((const)); // Signature[{i,j}]*Signature[{l,m}]
147
148 // fill 'free' array with indices which are not occupied by 'set' array
149 static void freeidxM3(int set[], int free[]);
150
151 static void Rescale(double factor);
152
153private:
154 static const unsigned char ti2[8];
155 static const unsigned char ti3[8];
156 static const unsigned char ti4[8];
157 static const unsigned char ti5[8];
158
159protected:
160 static const unsigned char idxtbl[64];
161
162 static const double teps; // expansion target accuracy
163 static const double heps; // near double precision eps
164
165 static const double ceps;
166
167 static const double deps1;
168 static const double deps2;
169 static const double deps3;
170
171 static const double seps1; // two-point cutoff
172 static const double seps2; // two-point cutoff
173
174 static const double epsir1; // 3-point IR cutoff
175 static const double epsir2; // 3-point IR cutoff
176
177 static double deps;
178 static double meps; // onshell cutoff
179 static double m3eps; // I3 IR cutoff
180};
181
182template <int N>
183class Minor : public MinorBase {
184public:
185 Minor() {}
186
187 inline double Kay(int i, int j) PURE__attribute__ ((pure)) {
188 if (i == 0)
189 {
190 return j == 0 ? 0 : 1;
191 } else {
192 return j == 0 ? 1 : Cay[ns(i, j)];
193 }
194 }
195
196protected:
197 // Cayley matrix (symmetric)
198 static const int DCay = N + 1;
199 double Cay[(DCay - 1) * (DCay) / 2];
200};
201
202class Minor5 : public Minor<5> {
203public:
204 friend class SPtr<Minor5>;
205 typedef SPtr<Minor5> Ptr;
206 static Ptr create(const Kinem5& k) { return Ptr(new Minor5(k)); }
207 static Ptr create(const Kinem4& k) { return Ptr(new Minor5(k)); }
208
209 ncomplex evalE(int ep);
210 ncomplex evalE(int ep, int i);
211 ncomplex evalE(int ep, int i, int j);
212 ncomplex evalE(int ep, int i, int j, int k);
213 ncomplex evalE(int ep, int i, int j, int k, int l);
214 ncomplex evalE(int ep, int i, int j, int k, int l, int m);
215
216#ifdef USE_GOLEM_MODE
217#ifdef USE_GOLEM_MODE_6
218 int psix;
219#define ix(i) i<psix ? i : i-1
220#else
221#define ix(i) i
222#endif
223 virtual ncomplex A(int ep) { return evalE(ep); }
224 virtual ncomplex A(int ep, int i) { return evalE(ep, ix(i)); }
225 virtual ncomplex A(int ep, int i, int j) { return evalE(ep, ix(i), ix(j)); }
226 virtual ncomplex A(int ep, int i, int j, int k) { return evalE(ep, ix(i), ix(j), ix(k)); }
227 virtual ncomplex A(int ep, int i, int j, int k, int l) { return evalE(ep, ix(i), ix(j), ix(k), ix(l)); }
228 virtual ncomplex A(int ep, int i, int j, int k, int l, int m) { return evalE(ep, ix(i), ix(j), ix(k), ix(l), ix(m)); }
229 virtual ncomplex B(int ep) { return evalE(ep, 0, 0); }
230 virtual ncomplex B(int ep, int i) { return evalE(ep, 0, 0, ix(i)); }
231 virtual ncomplex B(int ep, int i, int j) { return evalE(ep, 0, 0, ix(i), ix(j)); }
232 virtual ncomplex B(int ep, int i, int j, int k) { return evalE(ep, 0, 0, ix(i), ix(j), ix(k)); }
233 virtual ncomplex C(int ep) { return evalE(ep, 0, 0, 0, 0); }
234 virtual ncomplex C(int ep, int i) { return evalE(ep, 0, 0, 0, 0, ix(i)); }
235#undef ix
236#endif /* USE_GOLEM_MODE */
237
238 ncomplex I4s(int ep, int s);
239
240 ncomplex I3st(int ep, int s, int t);
241 ncomplex I4Ds(int ep, int s);
242 ncomplex I4Dsi(int ep, int s, int i);
243
244 ncomplex I2stu(int ep, int s, int t, int u);
245 ncomplex I3Dst(int ep, int s, int t);
246 ncomplex I4D2s(int ep, int s);
247
248 ncomplex I4D2si(int ep, int s, int i);
249 ncomplex I3Dsti(int ep, int s, int t, int i);
250 ncomplex I3D2st(int ep, int s, int t);
251 ncomplex I4D3s(int ep, int s);
252 ncomplex I4D2sij(int ep, int s, int i, int j);
253
254 ncomplex I2Dstu(int ep, int s, int t, int u);
255 ncomplex I3D2sti(int ep, int s, int t, int i);
256 ncomplex I4D3si(int ep, int s, int i);
257 ncomplex I4D3sij(int ep, int s, int i, int j);
258
259 ncomplex I2Dstui(int ep, int s, int t, int u, int i);
260 ncomplex I3D2stij(int ep, int s, int t, int i, int j);
261 ncomplex I4D3sijk(int ep, int s, int i, int j, int k);
262
263 ncomplex I4D4s(int ep, int s);
264 ncomplex I4D4si(int ep, int s, int i);
265 ncomplex I3D3sti(int ep, int s, int t, int i);
266 ncomplex I4D4sij(int ep, int s, int i, int j);
267 ncomplex I2D2stui(int ep, int s, int t, int u, int i);
268 ncomplex I3D3stij(int ep, int s, int t, int i, int j);
269 ncomplex I4D4sijk(int ep, int s, int i, int j, int k);
270 ncomplex I2D2stuij(int ep, int s, int t, int u, int i, int j);
271 ncomplex I3D3stijk(int ep, int s, int t, int i, int j, int k);
272 ncomplex I4D4sijkl(int ep, int s, int i, int j, int k, int l);
273
274 // Gram4
275 ncomplex I2D2stu(int ep, int s, int t, int u);
276 ncomplex I3D3st(int ep, int s, int t);
277 ncomplex I2D3stu(int ep, int s, int t, int u);
278 ncomplex I3D4st(int ep, int s, int t);
279 ncomplex I2D4stu(int ep, int s, int t, int u);
280 ncomplex I3D5st(int ep, int s, int t);
281 ncomplex I2D5stu(int ep, int s, int t, int u);
282 ncomplex I3D6st(int ep, int s, int t);
283 ncomplex I2D6stu(int ep, int s, int t, int u);
284 ncomplex I3D7st(int ep, int s, int t);
285
286 //Aux
287
288 double M1(int i, int l) PURE__attribute__ ((pure));
289 double M2(int i, int j, int l, int m) PURE__attribute__ ((pure));
290 double M3(int i, int j, int k, int l, int m, int n) PURE__attribute__ ((pure));
291
292 double gram3(double p1, double p2, double p3) PURE__attribute__ ((pure));
293
294private:
295 Minor5(const Kinem5& k); // prevent direct creation
296 Minor5(const Kinem4& k); // prevent direct creation
297 Minor5(const Minor5& m) : smax(0) { assert(0)(static_cast <bool> (0) ? void (0) : __assert_fail ("0"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.h", 297, __extension__
__PRETTY_FUNCTION__)); } // prevent copy-constructing
298 Minor5& operator= (const Minor5& m) { assert(0)(static_cast <bool> (0) ? void (0) : __assert_fail ("0"
, "generators/phokhara/phokhara/eemmg-lib/src/minor.h", 298, __extension__
__PRETTY_FUNCTION__)); } // prevent reassignment
299
300 Kinem5 kinem;
301 const int smax;
302
303 // Maximal elements of sub-Cayley's
304 double pmaxS4[DCay - 1]; // s{1..5}
305
306 double pmaxS3[10]; // s,t - symm 5x5 w/o diag els
307 double pmaxT3[10]; // s,t - symm 5x5 w/o diag els
308 double pmaxU3[10]; // s,t - symm 5x5 w/o diag els
309 int imax3[10];
310
311 double pmaxM2[10]; // i,ip - symm 5x5 w/o diag els
312
313 // flags marking evuation steps
314 enum EvalM {E_None = 0,
315 E_M1, E_M2, E_M3,
316 E_I4D4sijkl,
317 E_I4D4sijk = E_I4D4sijkl + 3,
318 E_I3D3stij = E_I4D4sijk + 3,
319 E_I4D4sij = E_I3D3stij + 3,
320 E_I3D3sti = E_I4D4sij + 3,
321 E_I4D4si = E_I3D3sti + 3,
322 E_I4D4s = E_I4D4si + 3,
323 E_I2D6stu = E_I4D4s + 3,
324 E_I3D7st = E_I2D6stu + 3,
325 E_I2D5stu = E_I3D7st + 3,
326 E_I3D6st = E_I2D5stu + 3,
327 E_I2D4stu = E_I3D6st + 3,
328 E_I3D5st = E_I2D4stu + 3,
329 E_I2D3stu = E_I3D5st + 3,
330 E_I3D4st = E_I2D3stu + 3,
331 E_I3D3stijk = E_I3D4st + 3,
332 E_I2D2stui = E_I3D3stijk + 3,
333 E_I2D2stuij = E_I2D2stui + 3,
334 E_I2D2stu = E_I2D2stuij + 3,
335 E_I3D3st = E_I2D2stu + 3,
336 E_I4D3sijk = E_I3D3st + 3,
337 E_I3D2stij = E_I4D3sijk + 3,
338 E_I2Dstui = E_I3D2stij + 3,
339 E_I3D2st = E_I2Dstui + 2,
340 E_I4D3s = E_I3D2st + 3,
341 E_I4D3sij = E_I4D3s + 3,
342 E_I4D3si = E_I4D3sij + 3,
343 E_I3D2sti = E_I4D3si + 3,
344 E_I2Dstu = E_I3D2sti + 3,
345 E_I3Dsti = E_I2Dstu + 3,
346 E_I4D2sij = E_I3Dsti + 3,
347 E_I2stu = E_I4D2sij + 3,
348 E_I4D2s = E_I2stu + 3,
349 E_I4D2si = E_I4D2s + 3,
350 E_I3Dst = E_I4D2si + 3,
351 E_I4Dsi = E_I3Dst + 3,
352 E_I4Ds = E_I4Dsi + 3,
353 E_I4s = E_I4Ds + 3,
354 E_I3st = E_I4s + 3,
355 E_DUM = E_I3st + 3,
356 E_LEN
357 };
358 std::bitset<E_LEN> fEval;
359
360 // number of unique ways you can scratch out 3 rows in 6x6 matrix
361 static const int DM3 = 20; // Binomial[6,3] = 6*5*4/3! = 20
362 double pM3[DM3 * (DM3 + 1) / 2];
363 static const int DM2 = 15; // Binomial[6,2] = 15
364 double pM2[DM2 * (DM2 + 1) / 2];
365 static const int DM1 = 6; // Binomial[6,1] = 6
366 double pM1[DM1 * (DM1 + 1) / 2];
367
368 // Integral values
369 ncomplex pI4s[3][DCay - 1]; // s{1..5}
370
371 ncomplex pI3st[3][10]; // symm 5x5 w/o diag els
372 ncomplex pI4Ds[1][DCay - 1]; // s{1..5} // finite
373 ncomplex pI4Dsi[3][CIDX(DCay-2)][DCay - 1]; // i{1..4}, s{1..5}
374
375 ncomplex pI2stu[2][10]; // symm 5x5x5 w/o diag els // (0,1) parts only
376 ncomplex pI3Dst[1][10]; // symm 5x5 w/o diag els 5*6-5=10 // finite part only
377 ncomplex pI4D2s[1][DCay - 1]; // s{1..5} // finite part only
378
379 ncomplex pI4D2si[1][CIDX(DCay-2)][DCay - 1]; // i{1..4} s{1..5} // finite
380 ncomplex pI3Dsti[3][CIDX(DCay-2)][10]; // i{1..4} + symm 5x5 w/o diag els
381 ncomplex pI4D2sij[3][CIDX(DCay-2) * (CIDX(DCay-2) + 1) / 2][DCay - 1]; // symm 4x4, s{1..5}
382
383 ncomplex pI2Dstu[2][10]; // symm 5x5x5 w/o diag els // (0,1) parts only
384 ncomplex pI3D2st[2][10]; // symm 5x5 w/o diag els // (0,1) parts only
385 ncomplex pI4D3s[2][DCay - 1]; // s{1..5} // (0,1) parts only
386 ncomplex pI3D2sti[2][CIDX(DCay-2)][10]; // i{1..4} + symm 5x5 w/o diag els // (0,1) parts only
387 ncomplex pI4D3si[2][CIDX(DCay-2)][DCay - 1]; // i{1..4} s{1..5} // (0,1) parts only
388 ncomplex pI4D3sij[1][CIDX(DCay-2) * (CIDX(DCay-2) + 1) / 2][DCay - 1]; // symm 4x4, s{1..5} // finite
389
390 ncomplex pI2Dstui[2][CIDX(DCay-2)][DCay - 1]; // ~~ 16 elements // finite part only
391 ncomplex pI3D2stij[3][CIDX(DCay-2) * (CIDX(DCay-2) + 1) / 2][10]; // 'symm 4x4' x 'symm 5x5 w/o diag els'
392 ncomplex pI4D3sijk[3][CIDX(DCay-2) * (CIDX(DCay-2) + 1) * (CIDX(DCay-2) + 2) / 6][DCay - 1]; // 4x4x4 symm(20 els), s{1..5}
393
394 ncomplex pI2D2stu[2][10]; // symm 5x5x5 w/o diag els // (0,1) parts only
395 ncomplex pI3D3st[2][10]; // symm 5x5 w/o diag els // (0,1) parts only
396 ncomplex pI4D4s[2][DCay - 1]; // s{1..5} // (0,1) parts only
397 ncomplex pI4D4si[2][CIDX(DCay-2)][DCay - 1]; // i{1..4} s{1..5} // (0,1) parts only
398 ncomplex pI3D3sti[2][CIDX(DCay-2)][10]; // i{1..4} + symm 5x5 w/o diag els // (0,1) parts only
399 ncomplex pI4D4sij[2][CIDX(DCay-2) * (CIDX(DCay-2) + 1) / 2][DCay - 1]; // symm 4x4, s{1..5} // (0,1) parts only
400 ncomplex pI2D2stui[2][CIDX(DCay-2)][DCay - 1]; // ~~ 16 elements // (0,1) parts only
401 ncomplex pI3D3stij[2][CIDX(DCay-2) * (CIDX(DCay-2) + 1) / 2][10]; // 'symm 4x4' x 'symm 5x5 w/o diag els' // (0,1) parts only
402 ncomplex pI4D4sijk[2][CIDX(DCay-2) * (CIDX(DCay-2) + 1) * (CIDX(DCay-2) + 2) / 6][DCay - 1]; // 4x4x4 symm(20 els), s{1..5} // finite
403 ncomplex pI2D2stuij[2][CIDX(DCay-2)][DCay - 1][2]; // as stui but with 0:i==j and 1:i!=j
404 ncomplex pI3D3stijk[3][CIDX(DCay-2) * (CIDX(DCay-2) + 1) * (CIDX(DCay-2) + 2) / 6][10]; // 4x4x4 symm(20 els), x 'symm 5x5 w/o diag els'
405 ncomplex pI4D4sijkl[3][CIDX(DCay-2) * (CIDX(DCay-2) + 1) * (CIDX(DCay-2) + 2) * (CIDX(DCay-2) + 3) / 24][DCay - 1]; // 4x4x4x4 symm(35 els), s{1..5} // ????
406
407 // Gram4
408 // TODO cleanup these funcs
409 ncomplex pI2D3stu[2][10]; // symm 5x5x5 w/o diag els // (0,1) parts only
410 ncomplex pI3D4st[2][10]; // symm 5x5 w/o diag els // (0,1) parts only
411 ncomplex pI2D4stu[2][10]; // symm 5x5x5 w/o diag els // (0,1) parts only
412 ncomplex pI3D5st[2][10]; // symm 5x5 w/o diag els // (0,1) parts only
413 ncomplex pI2D5stu[2][10]; // symm 5x5x5 w/o diag els // (0,1) parts only
414 ncomplex pI3D6st[2][10]; // symm 5x5 w/o diag els // (0,1) parts only
415 ncomplex pI2D6stu[2][10]; // symm 5x5x5 w/o diag els // (0,1) parts only
416 ncomplex pI3D7st[2][10]; // symm 5x5 w/o diag els // (0,1) parts only
417
418 // Aux
419
420 // evaluate and save for future use M1,M2,M3 minors
421 void evalM1();
422 void evalM2();
423 void evalM3();
424
425 // find and save maximal elements of subcayley matrices
426 void maxCay();
427
428 // evaluate and save for the future use scratched integrals
429 void I4sEval(int ep);
430
431 void I3stEval(int ep);
432 void I4DsEval(int ep);
433 void I4DsiEval(int ep);
434
435 void I2stuEval(int ep);
436 void I3DstEval(int ep);
437 void I4D2sEval(int ep);
438
439 void I4D2siEval(int ep);
440 void I3DstiEval(int ep);
441 void I3D2stEval(int ep);
442 void I4D3sEval(int ep);
443 void I4D2sijEval(int ep);
444
445 void I2DstuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq);
446 void I3D2stiEval(int ep);
447 void I4D3siEval(int ep);
448 void I4D3sijEval(int ep);
449
450 void I2DstuiEval(int ep, int s, int t, int u, int i, int ip, double qsq);
451 void I3D2stijEval(int ep);
452 void I4D3sijkEval(int ep);
453
454 void I4D4sEval(int ep);
455 void I4D4siEval(int ep);
456 void I3D3stiEval(int ep);
457 void I4D4sijEval(int ep);
458 void I2D2stuiEval(int ep, int s, int t, int u, int i, int ip, double qsq);
459 void I3D3stijEval(int ep);
460 void I4D4sijkEval(int ep);
461 void I2D2stuijEval(int ep, int s, int t, int u, int i, int ip, double qsq);
462 void I3D3stijkEval(int ep);
463 void I4D4sijklEval(int ep);
464
465 // Gram4
466 void I2D2stuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq);
467 void I3D3stEval(int ep);
468 void I2D3stuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq);
469 void I3D4stEval(int ep);
470 void I2D4stuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq);
471 void I3D5stEval(int ep);
472 void I2D5stuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq);
473 void I3D6stEval(int ep);
474 void I2D6stuEval(int idx, int ep, int s, int t, int u, int m, int n, double qsq);
475 void I3D7stEval(int ep);
476
477 // IR
478 ncomplex I2mDstu(int ep, int s, int t, int u, int m, int n);
479 ncomplex I2stui(int ep, int s, int t, int u, int i, int ip);
480
481 double M4ii(int u, int v, int i) PURE__attribute__ ((pure));
482 double M4ui(int u, int v, int i) PURE__attribute__ ((pure));
483 double M4vi(int u, int v, int i) PURE__attribute__ ((pure));
484 double M4uu(int u, int v, int i) PURE__attribute__ ((pure));
485 double M4vu(int u, int v, int i) PURE__attribute__ ((pure));
486 double M4vv(int u, int v, int i) PURE__attribute__ ((pure));
487};
488
489class Minor4 : public Minor<4> {
490public:
491 friend class SPtr<Minor4>;
492 typedef SPtr<Minor4> Ptr;
493 static Ptr create(const Kinem4& k, Minor5::Ptr mptr5, int s, int is)
494 {
495 return Ptr(new Minor4(k, mptr5, s, is));
2
←
Calling '~SPtr'
→
6
←
Returning from '~SPtr'
→
496 }
497
498 ncomplex evalD(int ep);
499 ncomplex evalD(int ep, int i);
500 ncomplex evalD(int ep, int i, int j);
501 ncomplex evalD(int ep, int i, int j, int k);
502 ncomplex evalD(int ep, int i, int j, int k, int l);
503
504#ifdef USE_GOLEM_MODE
505 virtual ncomplex A(int ep);
506 virtual ncomplex A(int ep, int i);
507 virtual ncomplex A(int ep, int i, int j);
508 virtual ncomplex A(int ep, int i, int j, int k);
509 virtual ncomplex A(int ep, int i, int j, int k, int l);
510 virtual ncomplex B(int ep);
511 virtual ncomplex B(int ep, int i);
512 virtual ncomplex B(int ep, int i, int j);
513 virtual ncomplex C(int ep);
514#endif /* USE_GOLEM_MODE */
515
516private:
517 Minor4(const Kinem4& k, Minor5::Ptr mptr, int s, int is);
518
519 Kinem4 kinem;
520
521 const Minor5::Ptr pm5;
522 const int ps, offs;
523};
524
525class Minor3 : public Minor<3> {
526public:
527 friend class SPtr<Minor3>;
528 typedef SPtr<Minor3> Ptr;
529 static Ptr create(const Kinem3& k, Minor5::Ptr mptr5, int s, int t, int is)
530 {
531 return Ptr(new Minor3(k, mptr5, s, t, is));
532 }
533
534 ncomplex evalC(int ep);
535 ncomplex evalC(int ep, int i);
536 ncomplex evalC(int ep, int i, int j);
537 ncomplex evalC(int ep, int i, int j, int k);
538
539#ifdef USE_GOLEM_MODE
540 virtual ncomplex A(int ep);
541 virtual ncomplex A(int ep, int i);
542 virtual ncomplex A(int ep, int i, int j);
543 virtual ncomplex A(int ep, int i, int j, int k);
544 virtual ncomplex B(int ep);
545 virtual ncomplex B(int ep, int i);
546#endif /* USE_GOLEM_MODE */
547
548private:
549 Minor3(const Kinem3& k);
550 Minor3(const Kinem3& k, Minor5::Ptr mptr5, int s, int t, int is);
551
552 Kinem3 kinem;
553
554 const Minor5::Ptr pm5;
555 const int ps, pt, offs;
556};
557
558class Minor2 : public Minor<2> {
559public:
560 friend class SPtr<Minor2>;
561 typedef SPtr<Minor2> Ptr;
562 static Ptr create(const Kinem2& k, Minor5::Ptr mptr5, int s, int t, int u, int is)
563 {
564 return Ptr(new Minor2(k, mptr5, s, t, u, is));
565 }
566
567 ncomplex evalB(int ep);
568 ncomplex evalB(int ep, int i);
569 ncomplex evalB(int ep, int i, int j);
570
571#ifdef USE_GOLEM_MODE
572 virtual ncomplex A(int ep);
573 virtual ncomplex A(int ep, int i);
574 virtual ncomplex A(int ep, int i, int j);
575 virtual ncomplex B(int ep);
576#endif /* USE_GOLEM_MODE */
577
578private:
579 Minor2(const Kinem2& k);
580 Minor2(const Kinem2& k, Minor5::Ptr mptr5, int s, int t, int u, int is);
581
582 Kinem2 kinem;
583
584 const Minor5::Ptr pm5;
585 const int ps, pt, pu, offs;
586};
587
588
589/* ===============================================
590 *
591 * Utility functions
592 *
593 * ===============================================
594 */
595
596// Completely antysymmetric i,j,k size 6 matrix index
597inline
598int MinorBase::im3(int i, int j, int k)
599{
600 return idxtbl[(1 << i) | (1 << j) | (1 << k)];
601}
602
603// Completely antysymmetric i,j size 6 matrix index
604inline
605int MinorBase::im2(int i, int j)
606{
607 return idxtbl[(1 << i) | (1 << j)];
608}
609
610// Signature[{i,j,k}]*Signature[{l,m,n}]
611inline
612int MinorBase::signM3ud(int i, int j, int k, int l, int m, int n)
613{
614 int t = (j - i) * (k - j) * (k - i) * (m - l) * (n - m) * (n - l);
615 return t == 0 ? 0 : 2 * (t >> (sizeof(int) * 8 - 1)) + 1;
616}
617
618// Signature[{i,j}]*Signature[{l,m}]
619inline
620int MinorBase::signM2ud(int i, int j, int l, int m)
621{
622 int t = (j - i) * (m - l);
623 return t == 0 ? 0 : 2 * (t >> (sizeof(int) * 8 - 1)) + 1;
624}
625
626#endif /* QUL_MINOR_H */
←

generators/phokhara/phokhara/eemmg-lib/src/pointer.h

1/*
2 * pointer.h - smart pointer and array classes
3 *
4 * this file is part of PJFry library
5 * Copyright 2011 Valery Yundin
6 */
7
8#ifndef QUL_POINTER_H
9#define QUL_POINTER_H
10
11#include <memory>
12#include <cstring>
13
14class SRefCnt {
15public:
16 SRefCnt() : count(0) { }
17
18protected:
19 int count;
20};
21
22template <typename T>
23class SPtr {
24public:
25 T* operator-> () const { return pObj; }
26 bool operator== (const T* pobj) const { return pobj == pObj; }
27 bool operator!= (const T* pobj) const { return pobj != pObj; }
28
29 bool operator== (const SPtr<T>& spobj) const { return spobj.pObj == pObj; }
30 bool operator!= (const SPtr<T>& spobj) const { return spobj.pObj != pObj; }
31
32 SPtr(T* pobj = 0) : pObj(pobj)
33 {
34 if (pObj) { pObj->count++; }
35 }
36 SPtr(const SPtr& ptr) : pObj(ptr.pObj)
37 {
38 if (pObj) { pObj->count++; }
39 }
40
41 SPtr& operator= (const SPtr& ptr)
42 {
43 if (this == &ptr) { return *this; }
44 if (pObj && --(pObj->count) == 0) { delete pObj; }
45 if ((pObj = ptr.pObj)) { pObj->count++; }
46 return *this;
47 }
48
49 ~SPtr()
50 {
51 if (pObj
2.1
Field 'pObj' is non-null
8.1
Field 'pObj' is non-null
2.1
Field 'pObj' is non-null
8.1
Field 'pObj' is non-null
2.1
Field 'pObj' is non-null
8.1
Field 'pObj' is non-null
 && --(pObj->count) == 0) { delete pObj; }
3
←
Assuming the condition is true
→
4
←
Taking true branch
→
5
←
Memory is released
→
9
←
Use of memory after it is freed
52 }
53
54private:
55 T* pObj;
56};
57
58template <typename T, int N> class DArray;
59// Iterator for DArray class
60template <typename T, int N>
61class NIter {
62 friend class DArray<T, N>;
63public:
64 inline T& operator* () { return ptr[idx % N]; }
65 inline T* operator-> () { return &ptr[idx % N]; }
66 inline NIter& operator++ () { idx++; return *this;}
67 inline NIter& operator+= (int n) { idx += n; return *this;}
68
69 inline bool operator== (const NIter& iter) { return idx == iter.idx && ptr == iter.ptr; }
70 inline bool operator!= (const NIter& iter) { return idx != iter.idx || ptr != iter.ptr; }
71
72 NIter(T* begin, int last) : ptr(begin), idx(last) {}
73private:
74 T* ptr;
75 int idx;
76};
77
78// DArray - static array with stack-like iterator
79template <typename T, int N>
80class DArray {
81public:
82 DArray() : last(N), len(0) {}
83
84 typedef NIter<T, N> iterator;
85 iterator begin() { return iterator(elems, last); }
86 iterator end() { return iterator(elems, last + len); }
87
88 T& insert(const T& el)
89 {
90 len = (len == N ? len : len + 1);
91 last = ((last - 1) + N) % N;
92 elems[last] = el;
93 return elems[last];
94 }
95
96#ifdef USE_SMART_INSERT"1"
97 void remove(iterator& it)
98 {
99// assert(it.ptr==elems);
100 int i = it.idx % N;
101 elems[i] = T();
102 if (i >= last) {
103 memmove(&elems[last + 1], &elems[last], (i - last)*sizeof(T));
104 memset(&elems[last], 0, sizeof(T));
105 last = (last + 1) % N;
106 } else {
107 memmove(&elems[i], &elems[i + 1], (last - i - 1)*sizeof(T));
108 memset(&elems[last - 1], 0, sizeof(T));
109 }
110 len = len - 1;
111 }
112#endif
113
114 void reset()
115 {
116#ifndef USE_DIRTY_RESET"1"
117 for (int i = 0; i < len; i++) {
118 elems[i] = T();
119 }
120#endif
121 last = N;
122 len = 0;
123 }
124
125 static const int size = N;
126 const T& operator [](const int idx) const { return elems[idx]; }
127 T& operator [](const int idx) { return elems[idx]; }
128
129private:
130 T elems[N];
131 int last;
132 int len;
133};
134
135// CArray - simple array with trivial iterator
136template <typename T, int N>
137class CArray {
138public:
139 typedef std::auto_ptr<CArray> Ptr;
140
141 CArray(T dval = T())
142 {
143 for (iterator i = begin(); i != end(); ++i) {
144 *i = dval;
145 }
146 }
147
148 typedef T* iterator;
149 iterator begin() { return &elems[0]; }
150 iterator end() { return &elems[N]; }
151
152
153 static const int size = N;
154
155 const T& operator [](const int idx) const { return elems[idx]; }
156 T& operator [](const int idx) { return elems[idx]; }
157
158private:
159 T elems[N];
160};
161
162#endif /* QUL_POINTER_H */
163
164/*
165class Test: public SRefCnt
166{
167 public:
168 friend class SPtr<Test>;
169 typedef SPtr<Test> Ptr;
170
171 static Ptr create(int x) { return Ptr(new Test(x)); };
172
173 ~Test() { };
174
175 private:
176 Test() { };
177};
178
179*/