Belle II Software  release-06-01-15
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) \
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) \
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)
46 #else
47 /* else calculate all (including coefficient of 0'th chord) */
48 # define CIDX (DCay-1)
49 #endif
50 
51 class MinorBase : public SRefCnt {
52 public:
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 {
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 { // 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 {
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 {
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 {
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 { // ordered
118  assert(i <= j);
119  return i + j * (j + 1) / 2;
120  }
121 
122  inline static int iss(int i, int j, int k) CONST { // ordered
123  assert(i <= j&& j <= k);
124  return i + ti2[j] + ti3[k];
125  }
126 
127  inline static int iss(int i, int j, int k, int l) CONST { // ordered
128  assert(i <= j&& j <= k&& k <= l);
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 { // ordered
133  assert(i <= j&& j <= k&& k <= l&& l <= m);
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; // Antisymmetric index for "3 out of 6" minors
144  static int im2(int i, int j) 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; // Signature[{i,j,k}]*Signature[{l,m,n}]
146  static int signM2ud(int i, int j, int l, int m) 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 
153 private:
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 
159 protected:
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 
182 template <int N>
183 class Minor : public MinorBase {
184 public:
185  Minor() {}
186 
187  inline double Kay(int i, int j) 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 
196 protected:
197  // Cayley matrix (symmetric)
198  static const int DCay = N + 1;
199  double Cay[(DCay - 1) * (DCay) / 2];
200 };
201 
202 class Minor5 : public Minor<5> {
203 public:
204  friend class SPtr<Minor5>;
205  typedef SPtr<Minor5> Ptr;
206  static Ptr create(const Kinem5& k) { return Ptr(new Minor5(k)); }
207  static Ptr create(const Kinem4& k) { return Ptr(new Minor5(k)); }
208 
209  ncomplex evalE(int ep);
210  ncomplex evalE(int ep, int i);
211  ncomplex evalE(int ep, int i, int j);
212  ncomplex evalE(int ep, int i, int j, int k);
213  ncomplex evalE(int ep, int i, int j, int k, int l);
214  ncomplex evalE(int ep, int i, int j, int k, int l, int m);
215 
216 #ifdef USE_GOLEM_MODE
217 #ifdef USE_GOLEM_MODE_6
218  int psix;
219 #define ix(i) i<psix ? i : i-1
220 #else
221 #define ix(i) i
222 #endif
223  virtual ncomplex A(int ep) { return evalE(ep); }
224  virtual ncomplex A(int ep, int i) { return evalE(ep, ix(i)); }
225  virtual ncomplex A(int ep, int i, int j) { return evalE(ep, ix(i), ix(j)); }
226  virtual ncomplex A(int ep, int i, int j, int k) { return evalE(ep, ix(i), ix(j), ix(k)); }
227  virtual ncomplex A(int ep, int i, int j, int k, int l) { return evalE(ep, ix(i), ix(j), ix(k), ix(l)); }
228  virtual ncomplex A(int ep, int i, int j, int k, int l, int m) { return evalE(ep, ix(i), ix(j), ix(k), ix(l), ix(m)); }
229  virtual ncomplex B(int ep) { return evalE(ep, 0, 0); }
230  virtual ncomplex B(int ep, int i) { return evalE(ep, 0, 0, ix(i)); }
231  virtual ncomplex B(int ep, int i, int j) { return evalE(ep, 0, 0, ix(i), ix(j)); }
232  virtual ncomplex B(int ep, int i, int j, int k) { return evalE(ep, 0, 0, ix(i), ix(j), ix(k)); }
233  virtual ncomplex C(int ep) { return evalE(ep, 0, 0, 0, 0); }
234  virtual ncomplex C(int ep, int i) { return evalE(ep, 0, 0, 0, 0, ix(i)); }
235 #undef ix
236 #endif /* USE_GOLEM_MODE */
237 
238  ncomplex I4s(int ep, int s);
239 
240  ncomplex I3st(int ep, int s, int t);
241  ncomplex I4Ds(int ep, int s);
242  ncomplex I4Dsi(int ep, int s, int i);
243 
244  ncomplex I2stu(int ep, int s, int t, int u);
245  ncomplex I3Dst(int ep, int s, int t);
246  ncomplex I4D2s(int ep, int s);
247 
248  ncomplex I4D2si(int ep, int s, int i);
249  ncomplex I3Dsti(int ep, int s, int t, int i);
250  ncomplex I3D2st(int ep, int s, int t);
251  ncomplex I4D3s(int ep, int s);
252  ncomplex I4D2sij(int ep, int s, int i, int j);
253 
254  ncomplex I2Dstu(int ep, int s, int t, int u);
255  ncomplex I3D2sti(int ep, int s, int t, int i);
256  ncomplex I4D3si(int ep, int s, int i);
257  ncomplex I4D3sij(int ep, int s, int i, int j);
258 
259  ncomplex I2Dstui(int ep, int s, int t, int u, int i);
260  ncomplex I3D2stij(int ep, int s, int t, int i, int j);
261  ncomplex I4D3sijk(int ep, int s, int i, int j, int k);
262 
263  ncomplex I4D4s(int ep, int s);
264  ncomplex I4D4si(int ep, int s, int i);
265  ncomplex I3D3sti(int ep, int s, int t, int i);
266  ncomplex I4D4sij(int ep, int s, int i, int j);
267  ncomplex I2D2stui(int ep, int s, int t, int u, int i);
268  ncomplex I3D3stij(int ep, int s, int t, int i, int j);
269  ncomplex I4D4sijk(int ep, int s, int i, int j, int k);
270  ncomplex I2D2stuij(int ep, int s, int t, int u, int i, int j);
271  ncomplex I3D3stijk(int ep, int s, int t, int i, int j, int k);
272  ncomplex I4D4sijkl(int ep, int s, int i, int j, int k, int l);
273 
274  // Gram4
275  ncomplex I2D2stu(int ep, int s, int t, int u);
276  ncomplex I3D3st(int ep, int s, int t);
277  ncomplex I2D3stu(int ep, int s, int t, int u);
278  ncomplex I3D4st(int ep, int s, int t);
279  ncomplex I2D4stu(int ep, int s, int t, int u);
280  ncomplex I3D5st(int ep, int s, int t);
281  ncomplex I2D5stu(int ep, int s, int t, int u);
282  ncomplex I3D6st(int ep, int s, int t);
283  ncomplex I2D6stu(int ep, int s, int t, int u);
284  ncomplex I3D7st(int ep, int s, int t);
285 
286  //Aux
287 
288  double M1(int i, int l) PURE;
289  double M2(int i, int j, int l, int m) PURE;
290  double M3(int i, int j, int k, int l, int m, int n) PURE;
291 
292  double gram3(double p1, double p2, double p3) PURE;
293 
294 private:
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); } // prevent copy-constructing
298  Minor5& operator= (const Minor5& m) { assert(0); } // 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 - 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 - 1]; // i{1..4} s{1..5} // finite
380  ncomplex pI3Dsti[3][CIDX][10]; // i{1..4} + symm 5x5 w/o diag els
381  ncomplex pI4D2sij[3][CIDX * (CIDX + 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][10]; // i{1..4} + symm 5x5 w/o diag els // (0,1) parts only
387  ncomplex pI4D3si[2][CIDX][DCay - 1]; // i{1..4} s{1..5} // (0,1) parts only
388  ncomplex pI4D3sij[1][CIDX * (CIDX + 1) / 2][DCay - 1]; // symm 4x4, s{1..5} // finite
389 
390  ncomplex pI2Dstui[2][CIDX][DCay - 1]; // ~~ 16 elements // finite part only
391  ncomplex pI3D2stij[3][CIDX * (CIDX + 1) / 2][10]; // 'symm 4x4' x 'symm 5x5 w/o diag els'
392  ncomplex pI4D3sijk[3][CIDX * (CIDX + 1) * (CIDX + 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 - 1]; // i{1..4} s{1..5} // (0,1) parts only
398  ncomplex pI3D3sti[2][CIDX][10]; // i{1..4} + symm 5x5 w/o diag els // (0,1) parts only
399  ncomplex pI4D4sij[2][CIDX * (CIDX + 1) / 2][DCay - 1]; // symm 4x4, s{1..5} // (0,1) parts only
400  ncomplex pI2D2stui[2][CIDX][DCay - 1]; // ~~ 16 elements // (0,1) parts only
401  ncomplex pI3D3stij[2][CIDX * (CIDX + 1) / 2][10]; // 'symm 4x4' x 'symm 5x5 w/o diag els' // (0,1) parts only
402  ncomplex pI4D4sijk[2][CIDX * (CIDX + 1) * (CIDX + 2) / 6][DCay - 1]; // 4x4x4 symm(20 els), s{1..5} // finite
403  ncomplex pI2D2stuij[2][CIDX][DCay - 1][2]; // as stui but with 0:i==j and 1:i!=j
404  ncomplex pI3D3stijk[3][CIDX * (CIDX + 1) * (CIDX + 2) / 6][10]; // 4x4x4 symm(20 els), x 'symm 5x5 w/o diag els'
405  ncomplex pI4D4sijkl[3][CIDX * (CIDX + 1) * (CIDX + 2) * (CIDX + 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;
482  double M4ui(int u, int v, int i) PURE;
483  double M4vi(int u, int v, int i) PURE;
484  double M4uu(int u, int v, int i) PURE;
485  double M4vu(int u, int v, int i) PURE;
486  double M4vv(int u, int v, int i) PURE;
487 };
488 
489 class Minor4 : public Minor<4> {
490 public:
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 
516 private:
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 
525 class Minor3 : public Minor<3> {
526 public:
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 
548 private:
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 
558 class Minor2 : public Minor<2> {
559 public:
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 
578 private:
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
597 inline
598 int 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
604 inline
605 int 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}]
611 inline
612 int 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}]
619 inline
620 int 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 */
Kinem2
Definition: kinem.h:104
Kinem3
Definition: kinem.h:121
Kinem4
Definition: kinem.h:144
Kinem5
Definition: kinem.h:176
Minor2
Definition: minor.h:558
Minor3
Definition: minor.h:525
Minor4
Definition: minor.h:489
Minor5
Definition: minor.h:202
Minor
Definition: minor.h:183
SPtr
Definition: pointer.h:23