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