Logo ROOT   6.08/07
Reference Guide
RooMath.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 // -- CLASS DESCRIPTION [MISC] --
18 // RooMath is a singleton class implementing various mathematical
19 // functions not found in TMath, mostly involving complex algebra
20 //
21 //
22 
23 #include <complex>
24 #include <cmath>
25 #include <algorithm>
26 
27 #include "RooFit.h"
28 
29 #include "RooMath.h"
30 #include "Riostream.h"
31 #include "RooMsgService.h"
32 
33 using namespace std;
34 
36 ;
37 
38 
39 namespace faddeeva_impl {
40  static inline void cexp(double& re, double& im)
41  {
42  // with gcc on unix machines and on x86_64, we can gain by hand-coding
43  // exp(z) for the x87 coprocessor; other platforms have the default
44  // routines as fallback implementation, and compilers other than gcc on
45  // x86_64 generate better code with the default routines; also avoid
46  // the inline assembly code when the copiler is not optimising code, or
47  // is optimising for code size
48  // (we insist on __unix__ here, since the assemblers on other OSs
49  // typically do not speak AT&T syntax as gas does...)
50 #if !(defined(__GNUC__) || defined(__clang__)) || \
51  !defined(__unix__) || !defined(__x86_64__) || \
52  !defined(__OPTIMIZE__) || defined(__OPTIMIZE_SIZE__) || \
53  defined(__INTEL_COMPILER) || \
54  defined(__OPEN64__) || defined(__PATHSCALE__)
55  const double e = std::exp(re);
56  re = e * std::cos(im);
57  im = e * std::sin(im);
58 #else
59  __asm__ (
60  "fxam\n\t" // examine st(0): NaN? Inf?
61  "fstsw %%ax\n\t"
62  "movb $0x45,%%dh\n\t"
63  "andb %%ah,%%dh\n\t"
64  "cmpb $0x05,%%dh\n\t"
65  "jz 1f\n\t" // have NaN or infinity, handle below
66  "fldl2e\n\t" // load log2(e)
67  "fmulp\n\t" // re * log2(e)
68  "fld %%st(0)\n\t" // duplicate re * log2(e)
69  "frndint\n\t" // int(re * log2(e))
70  "fsubr %%st,%%st(1)\n\t" // st(1) = x = frac(re * log2(e))
71  "fxch\n\t" // swap st(0), st(1)
72  "f2xm1\n\t" // 2^x - 1
73  "fld1\n\t" // st(0) = 1
74  "faddp\n\t" // st(0) = 2^x
75  "fscale\n\t" // 2 ^ (int(re * log2(e)) + x)
76  "fstp %%st(1)\n\t" // pop st(1)
77  "jmp 2f\n\t"
78  "1:\n\t" // handle NaN, Inf...
79  "testl $0x200, %%eax\n\t"// -infinity?
80  "jz 2f\n\t"
81  "fstp %%st\n\t" // -Inf, so pop st(0)
82  "fldz\n\t" // st(0) = 0
83  "2:\n\t" // here. we have st(0) == exp(re)
84  "fxch\n\t" // st(0) = im, st(1) = exp(re)
85  "fsincos\n\t" // st(0) = cos(im), st(1) = sin(im)
86  "fnstsw %%ax\n\t"
87  "testl $0x400,%%eax\n\t"
88  "jz 4f\n\t" // |im| too large for fsincos?
89  "fldpi\n\t" // st(0) = pi
90  "fadd %%st(0)\n\t" // st(0) *= 2;
91  "fxch %%st(1)\n\t" // st(0) = im, st(1) = 2 * pi
92  "3:\n\t"
93  "fprem1\n\t" // st(0) = fmod(im, 2 * pi)
94  "fnstsw %%ax\n\t"
95  "testl $0x400,%%eax\n\t"
96  "jnz 3b\n\t" // fmod done?
97  "fstp %%st(1)\n\t" // yes, pop st(1) == 2 * pi
98  "fsincos\n\t" // st(0) = cos(im), st(1) = sin(im)
99  "4:\n\t" // all fine, fsincos succeeded
100  "fmul %%st(2)\n\t" // st(0) *= st(2)
101  "fxch %%st(2)\n\t" // st(2)=exp(re)*cos(im),st(0)=exp(im)
102  "fmulp %%st(1)\n\t" // st(1)=exp(re)*sin(im), pop st(0)
103  : "=t" (im), "=u" (re): "0" (re), "1" (im) :
104  "eax", "dh", "cc"
105 #ifndef __clang__
106  // normal compilers (like gcc) want to be told that we
107  // clobber x87 registers, even if we pop them afterwards
108  // (so they can make sure they don't save anything there)
109  , "st(5)", "st(6)", "st(7)"
110 #else // __clang__
111  // clang produces an error message with the clobber list
112  // from above - not sure why; it seems harmless to leave
113  // the popped x87 registers out of the clobber list for
114  // clang, and that is in fact what seems to be recommended
115  // here:
116  // http://lists.cs.uiuc.edu/pipermail/cfe-dev/2012-May/021715.html
117 #endif // __clang__
118  );
119 #endif
120  }
121 
122  template <class T, unsigned N, unsigned NTAYLOR, unsigned NCF>
123  static inline std::complex<T> faddeeva_smabmq_impl(
124  T zre, T zim, const T tm,
125  const T (&a)[N], const T (&npi)[N],
126  const T (&taylorarr)[N * NTAYLOR * 2])
127  {
128  // catch singularities in the Fourier representation At
129  // z = n pi / tm, and provide a Taylor series expansion in those
130  // points, and only use it when we're close enough to the real axis
131  // that there is a chance we need it
132  const T zim2 = zim * zim;
133  const T maxnorm = T(9) / T(1000000);
134  if (zim2 < maxnorm) {
135  // we're close enough to the real axis that we need to worry about
136  // singularities
137  const T dnsing = tm * zre / npi[1];
138  const T dnsingmax2 = (T(N) - T(1) / T(2)) * (T(N) - T(1) / T(2));
139  if (dnsing * dnsing < dnsingmax2) {
140  // we're in the interesting range of the real axis as well...
141  // deal with Re(z) < 0 so we only need N different Taylor
142  // expansions; use w(-x+iy) = conj(w(x+iy))
143  const bool negrez = zre < T(0);
144  // figure out closest singularity
145  const int nsing = int(std::abs(dnsing) + T(1) / T(2));
146  // and calculate just how far we are from it
147  const T zmnpire = std::abs(zre) - npi[nsing];
148  const T zmnpinorm = zmnpire * zmnpire + zim2;
149  // close enough to one of the singularities?
150  if (zmnpinorm < maxnorm) {
151  const T* coeffs = &taylorarr[nsing * NTAYLOR * 2];
152  // calculate value of taylor expansion...
153  // (note: there's no chance to vectorize this one, since
154  // the value of the next iteration depend on the ones from
155  // the previous iteration)
156  T sumre = coeffs[0], sumim = coeffs[1];
157  for (unsigned i = 1; i < NTAYLOR; ++i) {
158  const T re = sumre * zmnpire - sumim * zim;
159  const T im = sumim * zmnpire + sumre * zim;
160  sumre = re + coeffs[2 * i + 0];
161  sumim = im + coeffs[2 * i + 1];
162  }
163  // undo the flip in real part of z if needed
164  if (negrez) return std::complex<T>(sumre, -sumim);
165  else return std::complex<T>(sumre, sumim);
166  }
167  }
168  }
169  // negative Im(z) is treated by calculating for -z, and using the
170  // symmetry properties of erfc(z)
171  const bool negimz = zim < T(0);
172  if (negimz) {
173  zre = -zre;
174  zim = -zim;
175  }
176  const T znorm = zre * zre + zim2;
177  if (znorm > tm * tm) {
178  // use continued fraction approximation for |z| large
179  const T isqrtpi = 5.64189583547756287e-01;
180  const T z2re = (zre + zim) * (zre - zim);
181  const T z2im = T(2) * zre * zim;
182  T cfre = T(1), cfim = T(0), cfnorm = T(1);
183  for (unsigned k = NCF; k; --k) {
184  cfre = +(T(k) / T(2)) * cfre / cfnorm;
185  cfim = -(T(k) / T(2)) * cfim / cfnorm;
186  if (k & 1) cfre -= z2re, cfim -= z2im;
187  else cfre += T(1);
188  cfnorm = cfre * cfre + cfim * cfim;
189  }
190  T sumre = (zim * cfre - zre * cfim) * isqrtpi / cfnorm;
191  T sumim = -(zre * cfre + zim * cfim) * isqrtpi / cfnorm;
192  if (negimz) {
193  // use erfc(-z) = 2 - erfc(z) to get good accuracy for
194  // Im(z) < 0: 2 / exp(z^2) - w(z)
195  T ez2re = -z2re, ez2im = -z2im;
196  faddeeva_impl::cexp(ez2re, ez2im);
197  return std::complex<T>(T(2) * ez2re - sumre,
198  T(2) * ez2im - sumim);
199  } else {
200  return std::complex<T>(sumre, sumim);
201  }
202  }
203  const T twosqrtpi = 3.54490770181103205e+00;
204  const T tmzre = tm * zre, tmzim = tm * zim;
205  // calculate exp(i tm z)
206  T eitmzre = -tmzim, eitmzim = tmzre;
207  faddeeva_impl::cexp(eitmzre, eitmzim);
208  // form 1 +/- exp (i tm z)
209  const T numerarr[4] = {
210  T(1) - eitmzre, -eitmzim, T(1) + eitmzre, +eitmzim
211  };
212  // form tm z * (1 +/- exp(i tm z))
213  const T numertmz[4] = {
214  tmzre * numerarr[0] - tmzim * numerarr[1],
215  tmzre * numerarr[1] + tmzim * numerarr[0],
216  tmzre * numerarr[2] - tmzim * numerarr[3],
217  tmzre * numerarr[3] + tmzim * numerarr[2]
218  };
219  // common subexpressions for use inside the loop
220  const T reimtmzm2 = T(-2) * tmzre * tmzim;
221  const T imtmz2 = tmzim * tmzim;
222  const T reimtmzm22 = reimtmzm2 * reimtmzm2;
223  // on non-x86_64 architectures, when the compiler is producing
224  // unoptimised code and when optimising for code size, we use the
225  // straightforward implementation, but for x86_64, we use the
226  // brainf*cked code below that the gcc vectorizer likes to gain a few
227  // clock cycles; non-gcc compilers also get the normal code, since they
228  // usually do a better job with the default code (and yes, it's a pain
229  // that they're all pretending to be gcc)
230 #if (!defined(__x86_64__)) || !defined(__OPTIMIZE__) || \
231  defined(__OPTIMIZE_SIZE__) || defined(__INTEL_COMPILER) || \
232  defined(__clang__) || defined(__OPEN64__) || \
233  defined(__PATHSCALE__) || !defined(__GNUC__)
234  T sumre = (-a[0] / znorm) * (numerarr[0] * zre + numerarr[1] * zim);
235  T sumim = (-a[0] / znorm) * (numerarr[1] * zre - numerarr[0] * zim);
236  for (unsigned i = 0; i < N; ++i) {
237  const unsigned j = (i << 1) & 2;
238  // denominator
239  const T wk = imtmz2 + (npi[i] + tmzre) * (npi[i] - tmzre);
240  // norm of denominator
241  const T norm = wk * wk + reimtmzm22;
242  const T f = T(2) * tm * a[i] / norm;
243  // sum += a[i] * numer / wk
244  sumre -= f * (numertmz[j] * wk + numertmz[j + 1] * reimtmzm2);
245  sumim -= f * (numertmz[j + 1] * wk - numertmz[j] * reimtmzm2);
246  }
247 #else
248  // BEGIN fully vectorisable code - enjoy reading... ;)
249  T tmp[2 * N];
250  for (unsigned i = 0; i < N; ++i) {
251  const T wk = imtmz2 + (npi[i] + tmzre) * (npi[i] - tmzre);
252  tmp[2 * i + 0] = wk;
253  tmp[2 * i + 1] = T(2) * tm * a[i] / (wk * wk + reimtmzm22);
254  }
255  for (unsigned i = 0; i < N / 2; ++i) {
256  T wk = tmp[4 * i + 0], f = tmp[4 * i + 1];
257  tmp[4 * i + 0] = -f * (numertmz[0] * wk + numertmz[1] * reimtmzm2);
258  tmp[4 * i + 1] = -f * (numertmz[1] * wk - numertmz[0] * reimtmzm2);
259  wk = tmp[4 * i + 2], f = tmp[4 * i + 3];
260  tmp[4 * i + 2] = -f * (numertmz[2] * wk + numertmz[3] * reimtmzm2);
261  tmp[4 * i + 3] = -f * (numertmz[3] * wk - numertmz[2] * reimtmzm2);
262  }
263  if (N & 1) {
264  // we may have missed one element in the last loop; if so, process
265  // it now...
266  const T wk = tmp[2 * N - 2], f = tmp[2 * N - 1];
267  tmp[2 * (N - 1) + 0] = -f * (numertmz[0] * wk + numertmz[1] * reimtmzm2);
268  tmp[2 * (N - 1) + 1] = -f * (numertmz[1] * wk - numertmz[0] * reimtmzm2);
269  }
270  T sumre = (-a[0] / znorm) * (numerarr[0] * zre + numerarr[1] * zim);
271  T sumim = (-a[0] / znorm) * (numerarr[1] * zre - numerarr[0] * zim);
272  for (unsigned i = 0; i < N; ++i) {
273  sumre += tmp[2 * i + 0];
274  sumim += tmp[2 * i + 1];
275  }
276  // END fully vectorisable code
277 #endif
278  // prepare the result
279  if (negimz) {
280  // use erfc(-z) = 2 - erfc(z) to get good accuracy for
281  // Im(z) < 0: 2 / exp(z^2) - w(z)
282  const T z2im = -T(2) * zre * zim;
283  const T z2re = -(zre + zim) * (zre - zim);
284  T ez2re = z2re, ez2im = z2im;
285  faddeeva_impl::cexp(ez2re, ez2im);
286  return std::complex<T>(T(2) * ez2re + sumim / twosqrtpi,
287  T(2) * ez2im - sumre / twosqrtpi);
288  } else {
289  return std::complex<T>(-sumim / twosqrtpi, sumre / twosqrtpi);
290  }
291  }
292 
293  static const double npi24[24] = { // precomputed values n * pi
294  0.00000000000000000e+00, 3.14159265358979324e+00, 6.28318530717958648e+00,
295  9.42477796076937972e+00, 1.25663706143591730e+01, 1.57079632679489662e+01,
296  1.88495559215387594e+01, 2.19911485751285527e+01, 2.51327412287183459e+01,
297  2.82743338823081391e+01, 3.14159265358979324e+01, 3.45575191894877256e+01,
298  3.76991118430775189e+01, 4.08407044966673121e+01, 4.39822971502571053e+01,
299  4.71238898038468986e+01, 5.02654824574366918e+01, 5.34070751110264851e+01,
300  5.65486677646162783e+01, 5.96902604182060715e+01, 6.28318530717958648e+01,
301  6.59734457253856580e+01, 6.91150383789754512e+01, 7.22566310325652445e+01,
302  };
303  static const double a24[24] = { // precomputed Fourier coefficient prefactors
304  2.95408975150919338e-01, 2.75840233292177084e-01, 2.24573955224615866e-01,
305  1.59414938273911723e-01, 9.86657664154541891e-02, 5.32441407876394120e-02,
306  2.50521500053936484e-02, 1.02774656705395362e-02, 3.67616433284484706e-03,
307  1.14649364124223317e-03, 3.11757015046197600e-04, 7.39143342960301488e-05,
308  1.52794934280083635e-05, 2.75395660822107093e-06, 4.32785878190124505e-07,
309  5.93003040874588103e-08, 7.08449030774820423e-09, 7.37952063581678038e-10,
310  6.70217160600200763e-11, 5.30726516347079017e-12, 3.66432411346763916e-13,
311  2.20589494494103134e-14, 1.15782686262855879e-15, 5.29871142946730482e-17,
312  };
313  static const double taylorarr24[24 * 12] = {
314  // real part imaginary part, low order coefficients last
315  // nsing = 0
316  0.00000000000000000e-00, 3.00901111225470020e-01,
317  5.00000000000000000e-01, 0.00000000000000000e-00,
318  0.00000000000000000e-00, -7.52252778063675049e-01,
319  -1.00000000000000000e-00, 0.00000000000000000e-00,
320  0.00000000000000000e-00, 1.12837916709551257e+00,
321  1.00000000000000000e-00, 0.00000000000000000e-00,
322  // nsing = 1
323  -2.22423508493755319e-01, 1.87966717746229718e-01,
324  3.41805419240637628e-01, 3.42752593807919263e-01,
325  4.66574321730757753e-01, -5.59649213591058097e-01,
326  -8.05759710273191021e-01, -5.38989366115424093e-01,
327  -4.88914083733395200e-01, 9.80580906465856792e-01,
328  9.33757118080975970e-01, 2.82273885115127769e-01,
329  // nsing = 2
330  -2.60522586513312894e-01, -4.26259455096092786e-02,
331  1.36549702008863349e-03, 4.39243227763478846e-01,
332  6.50591493715480700e-01, -1.23422352472779046e-01,
333  -3.43379903564271318e-01, -8.13862662890748911e-01,
334  -7.96093943501906645e-01, 6.11271022503935772e-01,
335  7.60213717643090957e-01, 4.93801903948967945e-01,
336  // nsing = 3
337  -1.18249853727020186e-01, -1.90471659765411376e-01,
338  -2.59044664869706839e-01, 2.69333898502392004e-01,
339  4.99077838344125714e-01, 2.64644800189075006e-01,
340  1.26114512111568737e-01, -7.46519337025968199e-01,
341  -8.47666863706379907e-01, 1.89347715957263646e-01,
342  5.39641485816297176e-01, 5.97805988669631615e-01,
343  // nsing = 4
344  4.94825297066481491e-02, -1.71428212158876197e-01,
345  -2.97766677111471585e-01, 1.60773286596649656e-02,
346  1.88114210832460682e-01, 4.11734391195006462e-01,
347  3.98540613293909842e-01, -4.63321903522162715e-01,
348  -6.99522070542463639e-01, -1.32412024008354582e-01,
349  3.33997185986131785e-01, 6.01983450812696742e-01,
350  // nsing = 5
351  1.18367078448232332e-01, -6.09533063579086850e-02,
352  -1.74762998833038991e-01, -1.39098099222000187e-01,
353  -6.71534655984154549e-02, 3.34462251996496680e-01,
354  4.37429678577360024e-01, -1.59613865629038012e-01,
355  -4.71863911886034656e-01, -2.92759316465055762e-01,
356  1.80238737704018306e-01, 5.42834914744283253e-01,
357  // nsing = 6
358  8.87698096005701290e-02, 2.84339354980994902e-02,
359  -3.18943083830766399e-02, -1.53946887977045862e-01,
360  -1.71825061547624858e-01, 1.70734367410600348e-01,
361  3.33690792296469441e-01, 3.97048587678703930e-02,
362  -2.66422678503135697e-01, -3.18469797424381480e-01,
363  8.48049724711137773e-02, 4.60546329221462864e-01,
364  // nsing = 7
365  2.99767046276705077e-02, 5.34659695701718247e-02,
366  4.53131030251822568e-02, -9.37915401977138648e-02,
367  -1.57982359988083777e-01, 3.82170507060760740e-02,
368  1.98891589845251706e-01, 1.17546677047049354e-01,
369  -1.27514335237079297e-01, -2.72741112680307074e-01,
370  3.47906344595283767e-02, 3.82277517244493224e-01,
371  // nsing = 8
372  -7.35922494437203395e-03, 3.72011290318534610e-02,
373  5.66783220847204687e-02, -3.21015398169199501e-02,
374  -1.00308737825172555e-01, -2.57695148077963515e-02,
375  9.67294850588435368e-02, 1.18174625238337507e-01,
376  -5.21266530264988508e-02, -2.08850084114630861e-01,
377  1.24443217440050976e-02, 3.19239968065752286e-01,
378  // nsing = 9
379  -1.66126772808035320e-02, 1.46180329587665321e-02,
380  3.85927576915247303e-02, 1.18910471133003227e-03,
381  -4.94003498320899806e-02, -3.93468443660139110e-02,
382  3.92113167048952835e-02, 9.03306084789976219e-02,
383  -1.82889636251263500e-02, -1.53816215444915245e-01,
384  3.88103861995563741e-03, 2.72090310854550347e-01,
385  // nsing = 10
386  -1.21245068916826880e-02, 1.59080224420074489e-03,
387  1.91116222508366035e-02, 1.05879549199053302e-02,
388  -1.97228428219695318e-02, -3.16962067712639397e-02,
389  1.34110372628315158e-02, 6.18045654429108837e-02,
390  -5.52574921865441838e-03, -1.14259663804569455e-01,
391  1.05534036292203489e-03, 2.37326534898818288e-01,
392  // nsing = 11
393  -5.96835002183177493e-03, -2.42594931567031205e-03,
394  7.44753817476594184e-03, 9.33450807578394386e-03,
395  -6.52649522783026481e-03, -2.08165802069352019e-02,
396  3.89988065678848650e-03, 4.12784313451549132e-02,
397  -1.44110721106127920e-03, -8.76484782997757425e-02,
398  2.50210184908121337e-04, 2.11131066219336647e-01,
399  // nsing = 12
400  -2.24505212235034193e-03, -2.38114524227619446e-03,
401  2.36375918970809340e-03, 5.97324040603806266e-03,
402  -1.81333819936645381e-03, -1.28126250720444051e-02,
403  9.69251586187208358e-04, 2.83055679874589732e-02,
404  -3.24986363596307374e-04, -6.97056268370209313e-02,
405  5.17231862038123061e-05, 1.90681117197597520e-01,
406  // nsing = 13
407  -6.76887607549779069e-04, -1.48589685249767064e-03,
408  6.22548369472046953e-04, 3.43871156746448680e-03,
409  -4.26557147166379929e-04, -7.98854145009655400e-03,
410  2.06644460919535524e-04, 2.03107152586353217e-02,
411  -6.34563929410856987e-05, -5.71425144910115832e-02,
412  9.32252179140502456e-06, 1.74167663785025829e-01,
413  // nsing = 14
414  -1.67596437777156162e-04, -8.05384193869903178e-04,
415  1.37627277777023791e-04, 1.97652692602724093e-03,
416  -8.54392244879459717e-05, -5.23088906415977167e-03,
417  3.78965577556493513e-05, 1.52191559129376333e-02,
418  -1.07393019498185646e-05, -4.79347862153366295e-02,
419  1.46503970628861795e-06, 1.60471011683477685e-01,
420  // nsing = 15
421  -3.45715760630978778e-05, -4.31089554210205493e-04,
422  2.57350138106549737e-05, 1.19449262097417514e-03,
423  -1.46322227517372253e-05, -3.61303766799909378e-03,
424  5.99057675687392260e-06, 1.17993805017130890e-02,
425  -1.57660578509526722e-06, -4.09165023743669707e-02,
426  2.00739683204152177e-07, 1.48879348585662670e-01,
427  // nsing = 16
428  -5.99735188857573424e-06, -2.42949218855805052e-04,
429  4.09249090936269722e-06, 7.67400152727128171e-04,
430  -2.14920611287648034e-06, -2.60710519575546230e-03,
431  8.17591694958640978e-07, 9.38581640137393053e-03,
432  -2.00910914042737743e-07, -3.54045580123653803e-02,
433  2.39819738182594508e-08, 1.38916449405613711e-01,
434  // nsing = 17
435  -8.80708505155966658e-07, -1.46479474515521504e-04,
436  5.55693207391871904e-07, 5.19165587844615415e-04,
437  -2.71391142598826750e-07, -1.94439427580099576e-03,
438  9.64641799864928425e-08, 7.61536975207357980e-03,
439  -2.22357616069432967e-08, -3.09762939485679078e-02,
440  2.49806920458212581e-09, 1.30247401712293206e-01,
441  // nsing = 18
442  -1.10007111030476390e-07, -9.35886150886691786e-05,
443  6.46244096997824390e-08, 3.65267193418479043e-04,
444  -2.95175785569292542e-08, -1.48730955943961081e-03,
445  9.84949251974795537e-09, 6.27824679148707177e-03,
446  -2.13827217704781576e-09, -2.73545766571797965e-02,
447  2.26877724435352177e-10, 1.22627158810895267e-01,
448  // nsing = 19
449  -1.17302439957657553e-08, -6.24890956722053332e-05,
450  6.45231881609786173e-09, 2.64799907072561543e-04,
451  -2.76943921343331654e-09, -1.16094187847598385e-03,
452  8.71074689656480749e-10, 5.24514377390761210e-03,
453  -1.78730768958639407e-10, -2.43489203319091538e-02,
454  1.79658223341365988e-11, 1.15870972518909888e-01,
455  // nsing = 20
456  -1.07084502471985403e-09, -4.31515421260633319e-05,
457  5.54152563270547927e-10, 1.96606443937168357e-04,
458  -2.24423474431542338e-10, -9.21550077887211094e-04,
459  6.67734377376211580e-11, 4.43201203646827019e-03,
460  -1.29896907717633162e-11, -2.18236356404862774e-02,
461  1.24042409733678516e-12, 1.09836276968151848e-01,
462  // nsing = 21
463  -8.38816525569060600e-11, -3.06091807093959821e-05,
464  4.10033961556230842e-11, 1.48895624771753491e-04,
465  -1.57238128435253905e-11, -7.42073499862065649e-04,
466  4.43938379112418832e-12, 3.78197089773957382e-03,
467  -8.21067867869285873e-13, -1.96793607299577220e-02,
468  7.46725770201828754e-14, 1.04410965521273064e-01,
469  // nsing = 22
470  -5.64848922712870507e-12, -2.22021942382507691e-05,
471  2.61729537775838587e-12, 1.14683068921649992e-04,
472  -9.53316139085394895e-13, -6.05021573565916914e-04,
473  2.56116039498542220e-13, 3.25530796858307225e-03,
474  -4.51482829896525004e-14, -1.78416955716514289e-02,
475  3.91940313268087086e-15, 9.95054815464739996e-02,
476  // nsing = 23
477  -3.27482357793897640e-13, -1.64138890390689871e-05,
478  1.44278798346454523e-13, 8.96362542918265398e-05,
479  -5.00524303437266481e-14, -4.98699756861136127e-04,
480  1.28274026095767213e-14, 2.82359118537843949e-03,
481  -2.16009593993917109e-15, -1.62538825704327487e-02,
482  1.79368667683853708e-16, 9.50473084594884184e-02
483  };
484 
485  const double npi11[11] = { // precomputed values n * pi
486  0.00000000000000000e+00, 3.14159265358979324e+00, 6.28318530717958648e+00,
487  9.42477796076937972e+00, 1.25663706143591730e+01, 1.57079632679489662e+01,
488  1.88495559215387594e+01, 2.19911485751285527e+01, 2.51327412287183459e+01,
489  2.82743338823081391e+01, 3.14159265358979324e+01
490  };
491  const double a11[11] = { // precomputed Fourier coefficient prefactors
492  4.43113462726379007e-01, 3.79788034073635143e-01, 2.39122407410867584e-01,
493  1.10599187402169792e-01, 3.75782250080904725e-02, 9.37936104296856288e-03,
494  1.71974046186334976e-03, 2.31635559000523461e-04, 2.29192401420125452e-05,
495  1.66589592139340077e-06, 8.89504561311882155e-08
496  };
497  const double taylorarr11[11 * 6] = {
498  // real part imaginary part, low order coefficients last
499  // nsing = 0
500  -1.00000000000000000e+00, 0.00000000000000000e+00,
501  0.00000000000000000e-01, 1.12837916709551257e+00,
502  1.00000000000000000e+00, 0.00000000000000000e+00,
503  // nsing = 1
504  -5.92741768247463996e-01, -7.19914991991294310e-01,
505  -6.73156763521649944e-01, 8.14025039279059577e-01,
506  8.57089811121701143e-01, 4.00248106586639754e-01,
507  // nsing = 2
508  1.26114512111568737e-01, -7.46519337025968199e-01,
509  -8.47666863706379907e-01, 1.89347715957263646e-01,
510  5.39641485816297176e-01, 5.97805988669631615e-01,
511  // nsing = 3
512  4.43238482668529408e-01, -3.03563167310638372e-01,
513  -5.88095866853990048e-01, -2.32638360700858412e-01,
514  2.49595637924601714e-01, 5.77633779156009340e-01,
515  // nsing = 4
516  3.33690792296469441e-01, 3.97048587678703930e-02,
517  -2.66422678503135697e-01, -3.18469797424381480e-01,
518  8.48049724711137773e-02, 4.60546329221462864e-01,
519  // nsing = 5
520  1.42043544696751869e-01, 1.24094227867032671e-01,
521  -8.31224229982140323e-02, -2.40766729258442100e-01,
522  2.11669512031059302e-02, 3.48650139549945097e-01,
523  // nsing = 6
524  3.92113167048952835e-02, 9.03306084789976219e-02,
525  -1.82889636251263500e-02, -1.53816215444915245e-01,
526  3.88103861995563741e-03, 2.72090310854550347e-01,
527  // nsing = 7
528  7.37741897722738503e-03, 5.04625223970221539e-02,
529  -2.87394336989990770e-03, -9.96122819257496929e-02,
530  5.22745478269428248e-04, 2.23361039070072101e-01,
531  // nsing = 8
532  9.69251586187208358e-04, 2.83055679874589732e-02,
533  -3.24986363596307374e-04, -6.97056268370209313e-02,
534  5.17231862038123061e-05, 1.90681117197597520e-01,
535  // nsing = 9
536  9.01625563468897100e-05, 1.74961124275657019e-02,
537  -2.65745127697337342e-05, -5.22070356354932341e-02,
538  3.75952450449939411e-06, 1.67018782142871146e-01,
539  // nsing = 10
540  5.99057675687392260e-06, 1.17993805017130890e-02,
541  -1.57660578509526722e-06, -4.09165023743669707e-02,
542  2.00739683204152177e-07, 1.48879348585662670e-01
543  };
544 }
545 
546 std::complex<double> RooMath::faddeeva(std::complex<double> z)
547 {
548  return faddeeva_impl::faddeeva_smabmq_impl<double, 24, 6, 9>(
549  z.real(), z.imag(), 12., faddeeva_impl::a24,
551 }
552 
553 std::complex<double> RooMath::faddeeva_fast(std::complex<double> z)
554 {
555  return faddeeva_impl::faddeeva_smabmq_impl<double, 11, 3, 3>(
556  z.real(), z.imag(), 8., faddeeva_impl::a11,
558 }
559 
560 std::complex<double> RooMath::erfc(const std::complex<double> z)
561 {
562  double re = -z.real() * z.real() + z.imag() * z.imag();
563  double im = -2. * z.real() * z.imag();
564  faddeeva_impl::cexp(re, im);
565  return (z.real() >= 0.) ?
566  (std::complex<double>(re, im) *
567  faddeeva(std::complex<double>(-z.imag(), z.real()))) :
568  (2. - std::complex<double>(re, im) *
569  faddeeva(std::complex<double>(z.imag(), -z.real())));
570 }
571 
572 std::complex<double> RooMath::erfc_fast(const std::complex<double> z)
573 {
574  double re = -z.real() * z.real() + z.imag() * z.imag();
575  double im = -2. * z.real() * z.imag();
576  faddeeva_impl::cexp(re, im);
577  return (z.real() >= 0.) ?
578  (std::complex<double>(re, im) *
579  faddeeva_fast(std::complex<double>(-z.imag(), z.real()))) :
580  (2. - std::complex<double>(re, im) *
581  faddeeva_fast(std::complex<double>(z.imag(), -z.real())));
582 }
583 
584 std::complex<double> RooMath::erf(const std::complex<double> z)
585 {
586  double re = -z.real() * z.real() + z.imag() * z.imag();
587  double im = -2. * z.real() * z.imag();
588  faddeeva_impl::cexp(re, im);
589  return (z.real() >= 0.) ?
590  (1. - std::complex<double>(re, im) *
591  faddeeva(std::complex<double>(-z.imag(), z.real()))) :
592  (std::complex<double>(re, im) *
593  faddeeva(std::complex<double>(z.imag(), -z.real())) - 1.);
594 }
595 
596 std::complex<double> RooMath::erf_fast(const std::complex<double> z)
597 {
598  double re = -z.real() * z.real() + z.imag() * z.imag();
599  double im = -2. * z.real() * z.imag();
600  faddeeva_impl::cexp(re, im);
601  return (z.real() >= 0.) ?
602  (1. - std::complex<double>(re, im) *
603  faddeeva_fast(std::complex<double>(-z.imag(), z.real()))) :
604  (std::complex<double>(re, im) *
605  faddeeva_fast(std::complex<double>(z.imag(), -z.real())) - 1.);
606 }
607 
608 
610 {
611  // Interpolate array 'ya' with 'n' elements for 'x' (between 0 and 'n'-1)
612 
613  // Int to Double conversion is faster via array lookup than type conversion!
614  static Double_t itod[20] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
615  10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0} ;
616  int i,m,ns=1 ;
617  Double_t den,dif,dift/*,ho,hp,w*/,y,dy ;
618  Double_t c[20], d[20] ;
619 
620  dif = fabs(x) ;
621  for(i=1 ; i<=n ; i++) {
622  if ((dift=fabs(x-itod[i-1]))<dif) {
623  ns=i ;
624  dif=dift ;
625  }
626  c[i] = ya[i-1] ;
627  d[i] = ya[i-1] ;
628  }
629 
630  y=ya[--ns] ;
631  for(m=1 ; m<n; m++) {
632  for(i=1 ; i<=n-m ; i++) {
633  den=(c[i+1] - d[i])/itod[m] ;
634  d[i]=(x-itod[i+m-1])*den ;
635  c[i]=(x-itod[i-1])*den ;
636  }
637  dy = (2*ns)<(n-m) ? c[ns+1] : d[ns--] ;
638  y += dy ;
639  }
640  return y ;
641 }
642 
643 
644 
646 {
647  // Interpolate array 'ya' with 'n' elements for 'xa'
648 
649  // Int to Double conversion is faster via array lookup than type conversion!
650 // static Double_t itod[20] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0,
651 // 10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0} ;
652  int i,m,ns=1 ;
653  Double_t den,dif,dift,ho,hp,w,y,dy ;
654  Double_t c[20], d[20] ;
655 
656  dif = fabs(x-xa[0]) ;
657  for(i=1 ; i<=n ; i++) {
658  if ((dift=fabs(x-xa[i-1]))<dif) {
659  ns=i ;
660  dif=dift ;
661  }
662  c[i] = ya[i-1] ;
663  d[i] = ya[i-1] ;
664  }
665 
666  y=ya[--ns] ;
667  for(m=1 ; m<n; m++) {
668  for(i=1 ; i<=n-m ; i++) {
669  ho=xa[i-1]-x ;
670  hp=xa[i-1+m]-x ;
671  w=c[i+1]-d[i] ;
672  den=ho-hp ;
673  if (den==0.) {
674  oocoutE((TObject*)0,Eval) << "RooMath::interpolate ERROR: zero distance between points not allowed" << endl ;
675  return 0 ;
676  }
677  den = w/den ;
678  d[i]=hp*den ;
679  c[i]=ho*den;
680  }
681  dy = (2*ns)<(n-m) ? c[ns+1] : d[ns--] ;
682  y += dy ;
683  }
684  return y ;
685 }
686 
687 
688 #include <map>
689 void RooMath::warn(const char* oldfun, const char* newfun)
690 {
691  static std::map<const char*, int> nwarn;
692  if (nwarn[oldfun] < 1<<12) {
693  ++nwarn[oldfun];
694  if (newfun) {
695  std::cout << "[#0] WARN: RooMath::" << oldfun <<
696  " is deprecated, please use " <<
697  newfun << " instead." << std::endl;
698  } else {
699  std::cout << "[#0] WARN: RooMath::" << oldfun <<
700  " is deprecated, and no longer needed, "
701  "you can remove the call to " <<
702  oldfun << " entirely." << std::endl;
703  }
704  }
705 }
static Double_t interpolate(Double_t yArr[], Int_t nOrder, Double_t x)
Definition: RooMath.cxx:609
static const double a24[24]
Definition: RooMath.cxx:303
static std::complex< T > faddeeva_smabmq_impl(T zre, T zim, const T tm, const T(&a)[N], const T(&npi)[N], const T(&taylorarr)[N *NTAYLOR *2])
Definition: RooMath.cxx:123
const double a11[11]
Definition: RooMath.cxx:491
const double npi11[11]
Definition: RooMath.cxx:485
return c
double T(double x)
Definition: ChebyshevPol.h:34
static void warn(const char *oldfun, const char *newfun=0)
Definition: RooMath.cxx:689
#define N
int Int_t
Definition: RtypesCore.h:41
TArc * a
Definition: textangle.C:12
STL namespace.
double cos(double)
Double_t x[n]
Definition: legend1.C:17
#define oocoutE(o, a)
Definition: RooMsgService.h:48
#define NCF(TN, I, C)
Definition: cfortran.h:897
static std::complex< double > faddeeva(std::complex< double > z)
evaluate Faddeeva function for complex argument
Definition: RooMath.cxx:546
const double taylorarr11[11 *6]
Definition: RooMath.cxx:497
double sin(double)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static std::complex< double > faddeeva_fast(std::complex< double > z)
evaluate Faddeeva function for complex argument (fast version)
Definition: RooMath.cxx:553
static std::complex< double > erfc_fast(const std::complex< double > z)
complex erfc function (fast version)
Definition: RooMath.cxx:572
TMarker * m
Definition: textangle.C:8
static std::complex< double > erfc(const std::complex< double > z)
complex erfc function
Definition: RooMath.cxx:560
static void cexp(double &re, double &im)
Definition: RooMath.cxx:40
static const double npi24[24]
Definition: RooMath.cxx:293
#define ClassImp(name)
Definition: Rtypes.h:279
double f(double x)
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Mother of all ROOT objects.
Definition: TObject.h:37
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:584
double exp(double)
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
static std::complex< double > erf_fast(const std::complex< double > z)
complex erf function (fast version)
Definition: RooMath.cxx:596
const Int_t n
Definition: legend1.C:16
static const double taylorarr24[24 *12]
Definition: RooMath.cxx:313