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