Bug Summary

File:generators/phokhara/phokhara/eemmg-lib/src/pointer.h
Warning:line 52, column 3
Potential leak of memory pointed to by field 'pObj'

Annotated Source Code

Press '?' to see keyboard shortcuts

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

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

→
1/*
2 * cache.h - cache classes
3 *
4 * this file is part of PJFry library
5 * Copyright 2011 Valery Yundin
6 */
7
8#include "cache.h"
9#include "integral.h"
10
11/* ------------------------------------------------------------
12 * ------------------------------------------------------------
13 * ICache section
14 * ------------------------------------------------------------
15 * ------------------------------------------------------------
16*/
17
18double ICache::mu2=1;
19
20double ICache::getMu2()
21{
22 return mu2;
23}
24
25double ICache::setMu2(const double newmu2)
26{
27 if (newmu2 > 0 && newmu2 != mu2) {
28 MinorBase::Rescale(newmu2/mu2);
29 mu2=newmu2;
30#ifdef USE_ONELOOP
31 double mu=sqrt(mu2);
32 F77_FUNC_(avh_olo_mu_set,AVH_OLO_MU_SET)avh_olo_mu_set_(&mu);
33#endif
34 ICache::Clear();
35 MCache::Clear();
36 }
37 return mu2;
38}
39
40void ICache::Clear()
41{
42 ClearCC();
43 ClearIC();
44}
45
46/* clear coefficient cache */
47void ICache::ClearCC()
48{
49 ic5[0].reset();
50 ic5[1].reset();
51 ic4[0].reset();
52 ic4[1].reset();
53 ic3[0].reset();
54 ic3[1].reset();
55 ic2[0].reset();
56 ic2[1].reset();
57}
58
59/* clear scalar integral cache */
60void ICache::ClearIC()
61{
62 ics4.reset();
63 ics3.reset();
64 ics2.reset();
65 ics1.reset();
66}
67
68/* clear minor cache */
69void MCache::Clear()
70{
71 cm5.reset();
72 cm4.reset();
73 cm3.reset();
74 cm2.reset();
75}
76
77// const double ICache::sNan=std::numeric_limits<double>::signaling_NaN();
78// const int64_t ICache::sNAN=0x7ffc0000BA13BA13LL;
79// const double ICache::sNAN.d64=reinterpret_cast<const double&>(sNAN); // Bad because breaks strict aliasing
80// const double ICache::sNAN.d64=((const ICache::ID64&)sNAN).dbl;
81
82const ICache::ID64 ICache::sNAN={ 0x7ffc0000BA13BA13LL };
83
84ICache::Array5 ICache::ic5[3];
85ICache::Array4 ICache::ic4[3];
86ICache::Array3 ICache::ic3[3];
87ICache::Array2 ICache::ic2[3];
88
89ICache::ArrayS4 ICache::ics4;
90ICache::ArrayS3 ICache::ics3;
91ICache::ArrayS2 ICache::ics2;
92ICache::ArrayS1 ICache::ics1;
93
94/* ===========================================================
95 *
96 * PENTAGON: 5 point coefficients
97 *
98 * ===========================================================
99 */
100
101/* --------------------------------------------------------
102 Rank-5 PENTAGON
103 * --------------------------------------------------------
104 */
105ncomplex ICache::getE(int ep, int i, int j, int k, int l, int m, const Kinem5 &kin)
106{
107#ifdef USE_CACHE_HIGH"1"
108 assert( (i==0 && j==0 && ( (k==0 && l==0 && m!=0) || (k!=0 && l!=0 && m!=0) ))(static_cast <bool> ((i==0 && j==0 && (
(k==0 && l==0 && m!=0) || (k!=0 && l
!=0 && m!=0) )) || (i!=0 && j!=0 && k
!=0 && l!=0 && m!=0)) ? void (0) : __assert_fail
("(i==0 && j==0 && ( (k==0 && l==0 && m!=0) || (k!=0 && l!=0 && m!=0) )) || (i!=0 && j!=0 && k!=0 && l!=0 && m!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 109
, __extension__ __PRETTY_FUNCTION__))
109 || (i!=0 && j!=0 && k!=0 && l!=0 && m!=0) )(static_cast <bool> ((i==0 && j==0 && (
(k==0 && l==0 && m!=0) || (k!=0 && l
!=0 && m!=0) )) || (i!=0 && j!=0 && k
!=0 && l!=0 && m!=0)) ? void (0) : __assert_fail
("(i==0 && j==0 && ( (k==0 && l==0 && m!=0) || (k!=0 && l!=0 && m!=0) )) || (i!=0 && j!=0 && k!=0 && l!=0 && m!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 109
, __extension__ __PRETTY_FUNCTION__));
110 int coefn=0;
111 if (i==0 && j==0) {
112 if (k==0 && l==0)
113 coefn=(m-1)+ee00001;
114 else
115 coefn=Minor5::is(k-1,l-1,m-1)+ee00111;
116 }
117 else {
118 coefn=Minor5::iss(i-1,j-1,k-1,l-1,m-1)+ee11111;
119 }
120 Save5 *s5=getS5(ep, kin, coefn);
121
122 ncomplex ivalue=(*s5)[coefn];
123 if (ivalue.real() == sNAN) {
124 ivalue=MCache::getMinor5(kin)->evalE(ep, i, j, k, l, m);
125 (*s5)[coefn]=ivalue;
126 }
127#else
128 ncomplex ivalue=MCache::getMinor5(kin)->evalE(ep, i, j, k, l, m);
129#endif /* USE_CACHE_HIGH */
130 return ivalue;
131}
132
133/* --------------------------------------------------------
134 Rank-4 PENTAGON
135 * --------------------------------------------------------
136 */
137ncomplex ICache::getE(int ep, int i, int j, int k, int l, const Kinem5 &kin)
138{
139#ifdef USE_CACHE_HIGH"1"
140 assert( (i==0 && j==0 && ( (k==0 && l==0) || (k!=0 && l!=0) )) || (i!=0 && j!=0 && k!=0 && l!=0) )(static_cast <bool> ((i==0 && j==0 && (
(k==0 && l==0) || (k!=0 && l!=0) )) || (i!=0
&& j!=0 && k!=0 && l!=0)) ? void (0)
: __assert_fail ("(i==0 && j==0 && ( (k==0 && l==0) || (k!=0 && l!=0) )) || (i!=0 && j!=0 && k!=0 && l!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 140
, __extension__ __PRETTY_FUNCTION__));
141 int coefn=0;
142 if (i==0 && j==0) {
143 if (k==0 && l==0)
144 coefn=ee0000;
145 else
146 coefn=Minor5::is(k-1,l-1)+ee0011;
147 }
148 else {
149 coefn=Minor5::is(i-1,j-1,k-1,l-1)+ee1111;
150 }
151 Save5 *s5=getS5(ep, kin, coefn);
152
153 ncomplex ivalue=(*s5)[coefn];
154 if (ivalue.real() == sNAN) {
155 ivalue=MCache::getMinor5(kin)->evalE(ep, i, j, k, l);
156 (*s5)[coefn]=ivalue;
157 }
158#else
159 ncomplex ivalue=MCache::getMinor5(kin)->evalE(ep, i, j, k, l);
160#endif /* USE_CACHE_HIGH */
161 return ivalue;
162}
163
164/* --------------------------------------------------------
165 Rank-3 PENTAGON
166 * --------------------------------------------------------
167 */
168ncomplex ICache::getE(int ep, int i, int j, int k, const Kinem5 &kin)
169{
170#ifdef USE_CACHE_HIGH"1"
171 assert( (i==0 && j==0 && k!=0) || (i!=0 && j!=0 && k!=0) )(static_cast <bool> ((i==0 && j==0 && k
!=0) || (i!=0 && j!=0 && k!=0)) ? void (0) : __assert_fail
("(i==0 && j==0 && k!=0) || (i!=0 && j!=0 && k!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 171
, __extension__ __PRETTY_FUNCTION__));
172 int coefn=(i==0 && j==0) ? (k-1)+ee001 : Minor5::is(i-1,j-1,k-1)+ee111;
173 Save5 *s5=getS5(ep, kin, coefn);
174
175 ncomplex ivalue=(*s5)[coefn];
176 if (ivalue.real() == sNAN) {
177 ivalue=MCache::getMinor5(kin)->evalE(ep, i, j, k);
178 (*s5)[coefn]=ivalue;
179 }
180#else
181 ncomplex ivalue=MCache::getMinor5(kin)->evalE(ep, i, j, k);
182#endif /* USE_CACHE_HIGH */
183 return ivalue;
184}
185
186/* --------------------------------------------------------
187 Rank-2 PENTAGON
188 * --------------------------------------------------------
189 */
190ncomplex ICache::getE(int ep, int i, int j, const Kinem5 &kin)
191{
192#ifdef USE_CACHE_HIGH"1"
193 assert( (i==0 && j==0) || (i!=0 && j!=0) )(static_cast <bool> ((i==0 && j==0) || (i!=0 &&
j!=0)) ? void (0) : __assert_fail ("(i==0 && j==0) || (i!=0 && j!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 193
, __extension__ __PRETTY_FUNCTION__));
194 int coefn=(i==0 && j==0) ? ee00 : Minor5::is(i-1,j-1)+ee11;
195 Save5 *s5=getS5(ep, kin, coefn);
196
197 ncomplex ivalue=(*s5)[coefn];
198 if (ivalue.real() == sNAN) {
199 ivalue=MCache::getMinor5(kin)->evalE(ep, i, j);
200 (*s5)[coefn]=ivalue;
201 }
202#else
203 ncomplex ivalue=MCache::getMinor5(kin)->evalE(ep, i, j);
204#endif /* USE_CACHE_HIGH */
205 return ivalue;
206}
207
208/* --------------------------------------------------------
209 Rank-1 PENTAGON
210 * --------------------------------------------------------
211 */
212ncomplex ICache::getE(int ep, int i, const Kinem5 &kin)
213{
214#ifdef USE_CACHE_HIGH"1"
215 int coefn=(i-1)+ee1;
216 Save5 *s5=getS5(ep, kin, coefn);
217
218 ncomplex ivalue=(*s5)[coefn];
219 if (ivalue.real() == sNAN) {
220 ivalue=MCache::getMinor5(kin)->evalE(ep, i);
221 (*s5)[coefn]=ivalue;
222 }
223#else
224 ncomplex ivalue=MCache::getMinor5(kin)->evalE(ep, i);
225#endif /* USE_CACHE_HIGH */
226 return ivalue;
227}
228
229/* --------------------------------------------------------
230 Rank-0 PENTAGON
231 * --------------------------------------------------------
232 */
233ncomplex ICache::getE(int ep, const Kinem5 &kin)
234{
235#ifdef USE_CACHE_HIGH"1"
236 int coefn=ee0;
237 Save5 *s5=getS5(ep, kin, coefn);
238
239 ncomplex ivalue=(*s5)[coefn];
240 if (ivalue.real() == sNAN) {
241 ivalue=MCache::getMinor5(kin)->evalE(ep);
242 (*s5)[coefn]=ivalue;
243 }
244#else
245 ncomplex ivalue=MCache::getMinor5(kin)->evalE(ep);
246#endif /* USE_CACHE_HIGH */
247 return ivalue;
248}
249
250/* ===========================================================
251 *
252 * BOX: 4 point coefficients
253 *
254 * ===========================================================
255 */
256
257/* --------------------------------------------------------
258 Rank-4 BOX
259 * --------------------------------------------------------
260 */
261ncomplex ICache::getD(int ep, int i, int j, int k, int l, const Kinem4 &kin)
262{
263#ifdef USE_CACHE_LOW
264 assert( (i==0 && j==0 && ( (k==0 && l==0) || (k!=0 && l!=0) )) || (i!=0 && j!=0 && k!=0 && l!=0) )(static_cast <bool> ((i==0 && j==0 && (
(k==0 && l==0) || (k!=0 && l!=0) )) || (i!=0
&& j!=0 && k!=0 && l!=0)) ? void (0)
: __assert_fail ("(i==0 && j==0 && ( (k==0 && l==0) || (k!=0 && l!=0) )) || (i!=0 && j!=0 && k!=0 && l!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 264
, __extension__ __PRETTY_FUNCTION__));
265 int coefn=0;
266 if (i==0 && j==0) {
267 if (k==0 && l==0)
268 coefn=dd0000;
269 else
270 coefn=Minor4::is(k-1,l-1)+dd0011;
271 }
272 else {
273 coefn=Minor5::is(i-1,j-1,k-1,l-1)+dd1111;
274 }
275 Save4 *s4=getS4(ep, kin, coefn);
276
277 ncomplex ivalue=(*s4)[coefn];
278 if (ivalue.real() == sNAN) {
279 ivalue=MCache::getMinor4(kin)->evalD(ep, i, j, k, l);
280 (*s4)[coefn]=ivalue;
281 }
282#else
283 ncomplex ivalue=MCache::getMinor4(kin)->evalD(ep, i, j, k, l);
284#endif /* USE_CACHE_LOW */
285 return ivalue;
286}
287
288/* --------------------------------------------------------
289 Rank-3 BOX
290 * --------------------------------------------------------
291 */
292ncomplex ICache::getD(int ep, int i, int j, int k, const Kinem4 &kin)
293{
294#ifdef USE_CACHE_LOW
295 assert( (i==0 && j==0 && k!=0) || (i!=0 && j!=0 && k!=0) )(static_cast <bool> ((i==0 && j==0 && k
!=0) || (i!=0 && j!=0 && k!=0)) ? void (0) : __assert_fail
("(i==0 && j==0 && k!=0) || (i!=0 && j!=0 && k!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 295
, __extension__ __PRETTY_FUNCTION__));
296 int coefn=(i==0 && j==0) ? (k-1)+dd001 : Minor4::is(i-1,j-1,k-1)+dd111;
297 Save4 *s4=getS4(ep, kin, coefn);
298
299 ncomplex ivalue=(*s4)[coefn];
300 if (ivalue.real() == sNAN) {
301 ivalue=MCache::getMinor4(kin)->evalD(ep, i, j, k);
302 (*s4)[coefn]=ivalue;
303 }
304#else
305 ncomplex ivalue=MCache::getMinor4(kin)->evalD(ep, i, j, k);
306#endif /* USE_CACHE_LOW */
307 return ivalue;
308}
309
310/* --------------------------------------------------------
311 Rank-2 BOX
312 * --------------------------------------------------------
313 */
314ncomplex ICache::getD(int ep, int i, int j, const Kinem4 &kin)
315{
316#ifdef USE_CACHE_LOW
317 assert( (i==0 && j==0) || (i!=0 && j!=0) )(static_cast <bool> ((i==0 && j==0) || (i!=0 &&
j!=0)) ? void (0) : __assert_fail ("(i==0 && j==0) || (i!=0 && j!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 317
, __extension__ __PRETTY_FUNCTION__));
318 int coefn=(i==0 && j==0) ? dd00 : Minor4::is(i-1,j-1)+dd11;
319 Save4 *s4=getS4(ep, kin, coefn);
320
321 ncomplex ivalue=(*s4)[coefn];
322 if (ivalue.real() == sNAN) {
323 ivalue=MCache::getMinor4(kin)->evalD(ep, i, j);
324 (*s4)[coefn]=ivalue;
325 }
326#else
327 ncomplex ivalue=MCache::getMinor4(kin)->evalD(ep, i, j);
328#endif /* USE_CACHE_LOW */
329 return ivalue;
330}
331
332/* --------------------------------------------------------
333 Rank-1 BOX
334 * --------------------------------------------------------
335 */
336ncomplex ICache::getD(int ep, int i, const Kinem4 &kin)
337{
338#ifdef USE_CACHE_LOW
339 int coefn=(i-1)+dd1;
340 Save4 *s4=getS4(ep, kin, coefn);
341
342 ncomplex ivalue=(*s4)[coefn];
343 if (ivalue.real() == sNAN) {
344 ivalue=MCache::getMinor4(kin)->evalD(ep, i);
345 (*s4)[coefn]=ivalue;
346 }
347#else
348 ncomplex ivalue=MCache::getMinor4(kin)->evalD(ep, i);
1
Calling 'MCache::getMinor4'
→
349#endif /* USE_CACHE_LOW */
350 return ivalue;
351}
352
353/* ===========================================================
354 *
355 * TRIANGLE: 3 point coefficients
356 *
357 * ===========================================================
358 */
359
360
361/* --------------------------------------------------------
362 Rank-3 TRIANGLE
363 * --------------------------------------------------------
364 */
365ncomplex ICache::getC(int ep, int i, int j, int k, const Kinem3 &kin)
366{
367#ifdef USE_CACHE_LOW
368 assert( (i==0 && j==0 && k!=0) || (i!=0 && j!=0 && k!=0) )(static_cast <bool> ((i==0 && j==0 && k
!=0) || (i!=0 && j!=0 && k!=0)) ? void (0) : __assert_fail
("(i==0 && j==0 && k!=0) || (i!=0 && j!=0 && k!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 368
, __extension__ __PRETTY_FUNCTION__));
369 int coefn=(i==0 && j==0) ? (k-1)+cc001 : Minor3::is(i-1,j-1,k-1)+cc111;
370 Save3 *s3=getS3(ep, kin, coefn);
371
372 ncomplex ivalue=(*s3)[coefn];
373 Minor3::Ptr pm3;
374 if (ivalue.real() == sNAN && (pm3=MCache::getMinor3(kin))!=0) {
375 ivalue=pm3->evalC(ep, i, j, k);
376 (*s3)[coefn]=ivalue;
377 }
378#else
379 Minor3::Ptr pm3=MCache::getMinor3(kin);
380 ncomplex ivalue= ( pm3!=0 ? pm3->evalC(ep, i, j, k) : sNAN.d64 );
381#endif /* USE_CACHE_LOW */
382#ifndef NDEBUG
383 if (pm3==0) printf("C%d%d%d(%.18e,%.18e,%.18e,%e,%e,%e)=%f\n",i,j,k,kin.p1(),kin.p2(),kin.p3(),kin.m1(),kin.m2(),kin.m3(),ivalue.real());
384#endif
385 return ivalue;
386}
387
388/* --------------------------------------------------------
389 Rank-2 TRIANGLE
390 * --------------------------------------------------------
391 */
392ncomplex ICache::getC(int ep, int i, int j, const Kinem3 &kin)
393{
394#ifdef USE_CACHE_LOW
395 assert( (i==0 && j==0) || (i!=0 && j!=0) )(static_cast <bool> ((i==0 && j==0) || (i!=0 &&
j!=0)) ? void (0) : __assert_fail ("(i==0 && j==0) || (i!=0 && j!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 395
, __extension__ __PRETTY_FUNCTION__));
396 int coefn=(i==0 && j==0) ? cc00 : Minor3::is(i-1,j-1)+cc11;
397 Save3 *s3=getS3(ep, kin, coefn);
398
399 ncomplex ivalue=(*s3)[coefn];
400 Minor3::Ptr pm3;
401 if (ivalue.real() == sNAN && (pm3=MCache::getMinor3(kin))!=0) {
402 ivalue=pm3->evalC(ep, i, j);
403 (*s3)[coefn]=ivalue;
404 }
405#else
406 Minor3::Ptr pm3=MCache::getMinor3(kin);
407 ncomplex ivalue= ( pm3!=0 ? pm3->evalC(ep, i, j) : sNAN.d64 );
408#endif /* USE_CACHE_LOW */
409#ifndef NDEBUG
410 if (pm3==0) printf("C%d%d(%.18e,%.18e,%.18e,%e,%e,%e)=%f\n",i,j,kin.p1(),kin.p2(),kin.p3(),kin.m1(),kin.m2(),kin.m3(),ivalue.real());
411#endif
412 return ivalue;
413}
414
415/* --------------------------------------------------------
416 Rank-1 TRIANGLE
417 * --------------------------------------------------------
418 */
419ncomplex ICache::getC(int ep, int i, const Kinem3 &kin)
420{
421#ifdef USE_CACHE_LOW
422 int coefn=(i-1)+cc1;
423 Save3 *s3=getS3(ep, kin, coefn);
424
425 ncomplex ivalue=(*s3)[coefn];
426 Minor3::Ptr pm3;
427 if (ivalue.real() == sNAN && (pm3=MCache::getMinor3(kin))!=0) {
428 ivalue=pm3->evalC(ep, i);
429 (*s3)[coefn]=ivalue;
430 }
431#else
432 Minor3::Ptr pm3=MCache::getMinor3(kin);
433 ncomplex ivalue= ( pm3!=0 ? pm3->evalC(ep, i) : sNAN.d64 );
434#endif /* USE_CACHE_LOW */
435#ifndef NDEBUG
436 if (pm3==0) printf("C%d(%.18e,%.18e,%.18e,%e,%e,%e)=%f\n",i,kin.p1(),kin.p2(),kin.p3(),kin.m1(),kin.m2(),kin.m3(),ivalue.real());
437#endif
438 return ivalue;
439}
440
441/* ===========================================================
442 *
443 * BUBBLE: 2 point coefficients
444 *
445 * ===========================================================
446 */
447
448/* --------------------------------------------------------
449 Rank-2 BUBBLE
450 * --------------------------------------------------------
451 */
452ncomplex ICache::getB(int ep, int i, int j, const Kinem2 &kin)
453{
454#ifdef USE_CACHE_LOW
455 assert( (i==0 && j==0) || (i!=0 && j!=0) )(static_cast <bool> ((i==0 && j==0) || (i!=0 &&
j!=0)) ? void (0) : __assert_fail ("(i==0 && j==0) || (i!=0 && j!=0)"
, "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp", 455
, __extension__ __PRETTY_FUNCTION__));
456 int coefn=(i==0 && j==0) ? bb00 : Minor2::is(i-1,j-1)+bb11;
457 Save2 *s2=getS2(ep, kin, coefn);
458
459 ncomplex ivalue=(*s2)[coefn];
460 Minor2::Ptr pm2;
461 if (ivalue.real() == sNAN && (pm2=MCache::getMinor2(kin))!=0) {
462 ivalue=pm2->evalB(ep, i, j);
463 (*s2)[coefn]=ivalue;
464 }
465#else
466 Minor2::Ptr pm2=MCache::getMinor2(kin);
467 ncomplex ivalue= ( pm2!=0 ? pm2->evalB(ep, i, j) : sNAN.d64 );
468#endif /* USE_CACHE_LOW */
469#ifndef NDEBUG
470 if (pm2==0) printf("B%d%d(%.18e,%.18e,%.18e)=%f\n",i,j,kin.p1(),kin.m1(),kin.m2(),ivalue.real());
471#endif
472 return ivalue;
473}
474
475/* --------------------------------------------------------
476 Rank-1 BUBBLE
477 * --------------------------------------------------------
478 */
479ncomplex ICache::getB(int ep, int i, const Kinem2 &kin)
480{
481#ifdef USE_CACHE_LOW
482 int coefn=(i-1)+bb1;
483 Save2 *s2=getS2(ep, kin, coefn);
484
485 ncomplex ivalue=(*s2)[coefn];
486 Minor2::Ptr pm2;
487 if (ivalue.real() == sNAN && (pm2=MCache::getMinor2(kin))!=0) {
488 ivalue=pm2->evalB(ep, i);
489 (*s2)[coefn]=ivalue;
490 }
491#else
492 Minor2::Ptr pm2=MCache::getMinor2(kin);
493 ncomplex ivalue= ( pm2!=0 ? pm2->evalB(ep, i) : sNAN.d64 );
494#endif /* USE_CACHE_LOW */
495#ifndef NDEBUG
496 if (pm2==0) printf("B%d(%.18e,%.18e,%.18e)=%f\n",i,kin.p1(),kin.m1(),kin.m2(),ivalue.real());
497#endif
498 return ivalue;
499}
500
501/* ===========================================================
502 *
503 * Get saved value
504 *
505 * ===========================================================
506 */
507
508// ICache::Save5* ICache::getS5(int ep, const Kinem5 &kin, int coefn)
509// {
510// Save5 *s5=0;
511// for (Array5::iterator it5=ic5[ep].begin(); it5 != ic5[ep].end(); ++it5) {
512// if (it5->key == kin) {
513// s5=it5->val.get();
514// break;
515// }
516// }
517// if (s5 == 0) {
518// Save5::Ptr sptr(new Save5(ncomplex(sNAN.d64, 0)));
519// s5=sptr.get();
520// ic5[ep].insert(Entry5(kin, sptr));
521// }
522// return s5;
523// }
524
525#define getSave(n) \
526ICache::Save##n * ICache::getS##n(int ep, const Kinem##n &kin, int coefn) \
527{ \
528 Save##n *s##n=0; \
529 for (Array##n::iterator it##n=ic##n[ep].begin(); it##n != ic##n[ep].end(); ++it##n) { \
530 if (it##n->key == kin) { \
531 s##n=it##n->val.get(); \
532 break; \
533 } \
534 } \
535 if (s##n == 0) { \
536 Save##n::Ptr sptr(new Save##n(ncomplex(sNAN.d64, 0))); \
537 s##n=sptr.get(); \
538 ic##n[ep].insert(Entry##n(kin, sptr)); \
539 } \
540 return s##n; \
541} \
542
543getSave(5)
544getSave(4)
545getSave(3)
546getSave(2)
547
548#undef getSave
549
550// ncomplex ICache::getI2(int ep, const Kinem2 &k)
551// {
552// ncomplex ivalue(sNAN.d64, 0);
553// for (ArrayS2::iterator it2=ics2[ep].begin(); it2 != ics2[ep].end(); ++it2) {
554// if (it2->key == k) {
555// ivalue=it2->val;
556// break;
557// }
558// }
559// if (ivalue.real() == sNAN) {
560// ivalue=qlI2(k.p1(),
561// k.m1(), k.m2(),
562// -ep);
563// ics2[ep].insert(EntryS2(k,ivalue));
564// }
565// return ivalue;
566// }
567
568#define getIN(n) \
569ncomplex ICache::getI##n(int ep, const Kinem##n &k) \
570{ \
571 Ival ivalue; \
572 bool found=false; \
573 for (ArrayS##n::iterator it##n=ics##n.begin(); it##n != ics##n.end(); ++it##n) { \
574 if (it##n->key == k) { \
575 ivalue=it##n->val; \
576 found=true; \
577 break; \
578 } \
579 } \
580 if ( ! found ) { \
581 ivalue=qlI##n(k); \
582 ics##n.insert(EntryS##n(k,ivalue)); \
583 } \
584 return ivalue.val[ep]; \
585}
586
587getIN(1)
588getIN(2)
589getIN(3)
590getIN(4)
591
592#undef getIN
593
594/* ------------------------------------------------------------
595 * ------------------------------------------------------------
596 * MCache section
597 * ------------------------------------------------------------
598 * ------------------------------------------------------------
599*/
600
601MCache::Array5 MCache::cm5;
602MCache::Array4 MCache::cm4;
603MCache::Array3 MCache::cm3;
604MCache::Array2 MCache::cm2;
605
606// Minor5::Ptr MCache::getMinor5(const Kinem5 &k)
607// {
608// Minor5::Ptr minor;
609// for (Array5::iterator it5=cm5.begin(); it5!=cm5.end(); ++it5) {
610// if (it5->key == k) {
611// minor=it5->val;
612// break;
613// }
614// }
615// if (minor==0) {
616// minor=Minor5::create(k);
617// cm5.insert(Entry5(k,minor));
618// }
619// return minor;
620// }
621
622#define getMinorN(n) \
623Minor##n::Ptr MCache::getMinor##n(const Kinem##n &k) \
624{ \
625 Minor##n::Ptr minor; \
626 for (Array##n::iterator it##n=cm##n.begin(); it##n!=cm##n.end(); ++it##n) { \
627 if (it##n->key == k) { \
628 minor=it##n->val; \
629 break; \
630 } \
631 } \
632 /* if (minor==0) { \
633 minor=Minor##n::create(k); \
634 cm##n.insert(Entry##n(k,minor)); \
635 } \
636 assert(minor!=0); */ \
637 return minor; \
638}
639
640getMinorN(3)
641getMinorN(2)
642
643#undef getMinorN
644
645Minor5::Ptr MCache::getMinor5(const Kinem5 &k)
646{
647 Minor5::Ptr minor;
648 for (Array5::iterator it5=cm5.begin(); it5!=cm5.end(); ++it5) {
649 if (it5->key == k) {
650 minor=it5->val;
651 break;
652 }
653 }
654 if (minor==0) {
655 minor=Minor5::create(k);
656 cm5.insert(Entry5(k,minor));
657 }
658 assert(minor!=0)(static_cast <bool> (minor!=0) ? void (0) : __assert_fail
("minor!=0", "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp"
, 658, __extension__ __PRETTY_FUNCTION__));
659 return minor;
660}
661
662Minor4::Ptr MCache::getMinor4(const Kinem4 &k)
663{
664 Minor4::Ptr minor;
665 for (Array4::iterator it4=cm4.begin(); it4!=cm4.end(); ++it4) {
2
←
Loop condition is false. Execution continues on line 671
→
666 if (it4->key == k) {
667 minor=it4->val;
668 break;
669 }
670 }
671 if (minor==0) {
3
←
Taking true branch
→
672 Minor5::create(k);
4
←
Calling 'Minor5::create'
→
6
←
Returned allocated memory
→
7
←
Calling '~SPtr'
→
673 minor=cm4.begin()->val;
674 cm4.insert(Entry4(k,minor));
675 }
676 assert(minor!=0)(static_cast <bool> (minor!=0) ? void (0) : __assert_fail
("minor!=0", "generators/phokhara/phokhara/eemmg-lib/src/cache.cpp"
, 676, __extension__ __PRETTY_FUNCTION__));
677 return minor;
678}
679
680
681#ifdef USE_SMART_INSERT"1"
682
683void MCache::smartinsertMinor3(const Kinem3 &k, Minor3::Ptr &m)
684{
685 for (Array3::iterator it3=cm3.begin(); it3!=cm3.end(); ++it3) {
686 if (it3->key == k) {
687 cm3.remove(it3);
688 break;
689 }
690 }
691 insertMinor3(k,m);
692}
693
694void MCache::smartinsertMinor2(const Kinem2 &k, Minor2::Ptr &m)
695{
696 for (Array2::iterator it2=cm2.begin(); it2!=cm2.end(); ++it2) {
697 if (it2->key == k) {
698 cm2.remove(it2);
699 break;
700 }
701 }
702 insertMinor2(k,m);
703}
704
705#endif
706
707// -------------------------------------------------------
←

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

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

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

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