ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
RooMath.h
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * File: $Id: RooMath.h,v 1.16 2007/05/11 09:11:30 verkerke Exp $
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 #ifndef ROO_MATH
17 #define ROO_MATH
18 
19 #include <cmath>
20 #include <complex>
21 
22 #include "Rtypes.h"
23 #include "TMath.h"
24 #include "RooComplex.h"
25 
26 #if defined(__my_func__)
27 #undef __my_func__
28 #endif
29 #if defined(WIN32)
30 #define __my_func__ __FUNCTION__
31 #else
32 #define __my_func__ __func__
33 #endif
34 
36 
37 class RooMath {
38 public:
39 
40  virtual ~RooMath() {};
41 
42  /** @brief evaluate Faddeeva function for complex argument
43  *
44  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
45  * @date 2013-02-21
46  *
47  * Calculate the value of the Faddeeva function @f$w(z) = \exp(-z^2)
48  * \mathrm{erfc}(-i z)@f$.
49  *
50  * The method described in
51  *
52  * S.M. Abrarov, B.M. Quine: "Efficient algotithmic implementation of
53  * Voigt/complex error function based on exponential series approximation"
54  * published in Applied Mathematics and Computation 218 (2011) 1894-1902
55  * doi:10.1016/j.amc.2011.06.072
56  *
57  * is used. At the heart of the method (equation (14) of the paper) is the
58  * following Fourier series based approximation:
59  *
60  * @f[ w(z) \approx \frac{i}{2\sqrt{\pi}}\left(
61  * \sum^N_{n=0} a_n \tau_m\left(
62  * \frac{1-e^{i(n\pi+\tau_m z)}}{n\pi + \tau_m z} -
63  * \frac{1-e^{i(-n\pi+\tau_m z)}}{n\pi - \tau_m z}
64  * \right) - a_0 \frac{1-e^{i \tau_m z}}{z}
65  * \right) @f]
66  *
67  * The coefficients @f$a_b@f$ are given by:
68  *
69  * @f[ a_n=\frac{2\sqrt{\pi}}{\tau_m}
70  * \exp\left(-\frac{n^2\pi^2}{\tau_m^2}\right) @f]
71  *
72  * To achieve machine accuracy in double precision floating point arithmetic
73  * for most of the upper half of the complex plane, chose @f$t_m=12@f$ and
74  * @f$N=23@f$ as is done in the paper.
75  *
76  * There are two complications: For Im(z) negative, the exponent in the
77  * equation above becomes so large that the roundoff in the rest of the
78  * calculation is amplified enough that the result cannot be trusted.
79  * Therefore, for Im(z) < 0, the symmetry of the erfc function under the
80  * transformation z --> -z is used to avoid accuracy issues for Im(z) < 0 by
81  * formulating the problem such that the calculation can be done for Im(z) > 0
82  * where the accuracy of the method is fine, and some postprocessing then
83  * yields the desired final result.
84  *
85  * Second, the denominators in the equation above become singular at
86  * @f$z = n * pi / 12@f$ (for 0 <= n < 24). In a tiny disc around these
87  * points, Taylor expansions are used to overcome that difficulty.
88  *
89  * This routine precomputes everything it can, and tries to write out complex
90  * operations to minimise subroutine calls, e.g. for the multiplication of
91  * complex numbers.
92  *
93  * In the square -8 <= Re(z) <= 8, -8 <= Im(z) <= 8, the routine is accurate
94  * to better than 4e-13 relative, the average relative error is better than
95  * 7e-16. On a modern x86_64 machine, the routine is roughly three times as
96  * fast than the old CERNLIB implementation and offers better accuracy.
97  *
98  * For large @f$|z|@f$, the familiar continued fraction approximation
99  *
100  * @f[ w(z)=\frac{-iz/\sqrt{\pi}}{-z^2+\frac{1/2}{1+\frac{2/2}{-z^2 +
101  * \frac{3/2}{1+\frac{4/2}{-z^2+\frac{5/2}{1+\frac{6/2}{-z^2+\frac{7/2
102  * }{1+\frac{8/2}{-z^2+\frac{9/2}{1+\ldots}}}}}}}}}} @f]
103  *
104  * is used, truncated at the ellipsis ("...") in the formula; for @f$|z| >
105  * 12@f$, @f$Im(z)>0@f$ it will give full double precision at a smaller
106  * computational cost than the method described above. (For @f$|z|>12@f$,
107  * @f$Im(z)<0@f$, the symmetry property @f$w(x-iy)=2e^{-(x+iy)^2-w(x+iy)}@f$
108  * is used.
109  */
110  static std::complex<double> faddeeva(std::complex<double> z);
111  /** @brief evaluate Faddeeva function for complex argument (fast version)
112  *
113  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
114  * @date 2013-02-21
115  *
116  * Calculate the value of the Faddeeva function @f$w(z) = \exp(-z^2)
117  * \mathrm{erfc}(-i z)@f$.
118  *
119  * This is the "fast" version of the faddeeva routine above. Fast means that
120  * is takes roughly half the amount of CPU of the slow version of the
121  * routine, but is a little less accurate.
122  *
123  * To be fast, chose @f$t_m=8@f$ and @f$N=11@f$ which should give accuracies
124  * around 1e-7.
125  *
126  * In the square -8 <= Re(z) <= 8, -8 <= Im(z) <= 8, the routine is accurate
127  * to better than 4e-7 relative, the average relative error is better than
128  * 5e-9. On a modern x86_64 machine, the routine is roughly five times as
129  * fast than the old CERNLIB implementation, or about 30% faster than the
130  * interpolation/lookup table based fast method used previously in RooFit,
131  * and offers better accuracy than the latter (the relative error is roughly
132  * a factor 280 smaller than the old interpolation/table lookup routine).
133  *
134  * For large @f$|z|@f$, the familiar continued fraction approximation
135  *
136  * @f[ w(z)=\frac{-iz/\sqrt{\pi}}{-z^2+\frac{1/2}{1+\frac{2/2}{-z^2 +
137  * \frac{3/2}{1+\ldots}}}} @f]
138  *
139  * is used, truncated at the ellipsis ("...") in the formula; for @f$|z| >
140  * 8@f$, @f$Im(z)>0@f$ it will give full float precision at a smaller
141  * computational cost than the method described above. (For @f$|z|>8@f$,
142  * @f$Im(z)<0@f$, the symmetry property @f$w(x-iy)=2e^{-(x+iy)^2-w(x+iy)}@f$
143  * is used.
144  */
145  static std::complex<double> faddeeva_fast(std::complex<double> z);
146 
147  /** @brief complex erf function
148  *
149  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
150  * @date 2013-02-21
151  *
152  * Calculate erf(z) for complex z.
153  */
154  static std::complex<double> erf(const std::complex<double> z);
155 
156  /** @brief complex erf function (fast version)
157  *
158  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
159  * @date 2013-02-21
160  *
161  * Calculate erf(z) for complex z. Use the code in faddeeva_fast to save some time.
162  */
163  static std::complex<double> erf_fast(const std::complex<double> z);
164  /** @brief complex erfc function
165  *
166  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
167  * @date 2013-02-21
168  *
169  * Calculate erfc(z) for complex z.
170  */
171  static std::complex<double> erfc(const std::complex<double> z);
172  /** @brief complex erfc function (fast version)
173  *
174  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
175  * @date 2013-02-21
176  *
177  * Calculate erfc(z) for complex z. Use the code in faddeeva_fast to save some time.
178  */
179  static std::complex<double> erfc_fast(const std::complex<double> z);
180 
181  // 1-D nth order polynomial interpolation routines
182  static Double_t interpolate(Double_t yArr[],Int_t nOrder, Double_t x) ;
183  static Double_t interpolate(Double_t xa[], Double_t ya[], Int_t n, Double_t x) ;
184 
185  static inline Double_t erf(Double_t x)
186  { return TMath::Erf(x); }
187 
188  static inline Double_t erfc(Double_t x)
189  { return TMath::Erfc(x); }
190 
191  /// deprecated function
193  { warn(__my_func__, "RooMath::faddeeva"); std::complex<Double_t> z = faddeeva(std::complex<Double_t>(re, im)); return RooComplex(z.real(), z.imag()); }
194  /// deprecated function
196  { warn(__my_func__, "RooMath::faddeeva"); std::complex<Double_t> z = faddeeva(std::complex<Double_t>(zz.re(), zz.im())); return RooComplex(z.real(), z.imag()); }
197  /// deprecated function
199  { warn(__my_func__, "RooMath::faddeeva_fast"); std::complex<Double_t> z = faddeeva_fast(std::complex<Double_t>(zz.re(), zz.im())); return RooComplex(z.real(), z.imag()); }
200  /// deprecated function
202  { warn(__my_func__, "RooMath::faddeeva_fast"); std::complex<Double_t> z = faddeeva_fast(std::complex<Double_t>(zz.re(), zz.im())); return z.real(); }
203  /// deprecated function
205  { warn(__my_func__, "RooMath::faddeeva_fast"); std::complex<Double_t> z = faddeeva_fast(std::complex<Double_t>(zz.re(), zz.im())); return z.imag(); }
206  /// deprecated function
208  { warn(__my_func__, "RooMath::faddeeva_fast"); std::complex<Double_t> z = faddeeva_fast(std::complex<Double_t>(zz.re(), zz.im())); return RooComplex(z.real(), z.imag()); }
209  /// deprecated function
211  { warn(__my_func__, "RooMath::faddeeva_fast"); std::complex<Double_t> z = faddeeva_fast(std::complex<Double_t>(zz.re(), zz.im())); return z.real(); }
212  /// deprecated function
214  { warn(__my_func__, "RooMath::faddeeva_fast"); std::complex<Double_t> z = faddeeva_fast(std::complex<Double_t>(zz.re(), zz.im())); return z.imag(); }
215  /// deprecated function
216  static void cacheCERF(Bool_t) { warn(__my_func__); }
217  /// deprecated function
218  static void cleanup() { warn(__my_func__); }
219  /// deprecated function
220  static void initFastCERF(Int_t /*reBins = 800*/, Double_t /*reMin = -4.0*/, Double_t /*reMax = 4.0*/,
221  Int_t /*imBins = 1000*/, Double_t /*imMin = -4.0*/, Double_t /*imMax = 6.0*/)
222  {
223  warn(__my_func__);
224  }
225 
226 private:
227  // deprecation warnings
228  static void warn(const char* oldfun, const char* newfun = 0);
229 
230  ClassDef(RooMath,0) // math utility routines
231 };
232 
233 #undef __my_func__
234 
235 #endif
static Double_t interpolate(Double_t yArr[], Int_t nOrder, Double_t x)
Definition: RooMath.cxx:609
#define __my_func__
Definition: RooMath.h:32
virtual ~RooMath()
Definition: RooMath.h:40
static void warn(const char *oldfun, const char *newfun=0)
Definition: RooMath.cxx:689
static Double_t ComplexErrFuncFastIm(const RooComplex &zz)
deprecated function
Definition: RooMath.h:204
Double_t re() const
Definition: RooComplex.h:79
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
Double_t x[n]
Definition: legend1.C:17
#define ClassDef(name, id)
Definition: Rtypes.h:254
static std::complex< double > faddeeva(std::complex< double > z)
evaluate Faddeeva function for complex argument
Definition: RooMath.cxx:546
Float_t z[5]
Definition: Ifit.C:16
static RooComplex ComplexErrFuncFast(const RooComplex &zz)
deprecated function
Definition: RooMath.h:198
Double_t Erfc(Double_t x)
Compute the complementary error function erfc(x).
Definition: TMath.cxx:197
static RooComplex ComplexErrFunc(Double_t re, Double_t im=0.)
deprecated function
Definition: RooMath.h:192
static std::complex< double > faddeeva_fast(std::complex< double > z)
evaluate Faddeeva function for complex argument (fast version)
Definition: RooMath.cxx:553
Double_t Erf(Double_t x)
Computation of the error function erf(x).
Definition: TMath.cxx:187
static void cleanup()
deprecated function
Definition: RooMath.h:218
static RooComplex ComplexErrFunc(const RooComplex &zz)
deprecated function
Definition: RooMath.h:195
static std::complex< double > erfc_fast(const std::complex< double > z)
complex erfc function (fast version)
Definition: RooMath.cxx:572
static void cacheCERF(Bool_t)
deprecated function
Definition: RooMath.h:216
static std::complex< double > erfc(const std::complex< double > z)
complex erfc function
Definition: RooMath.cxx:560
Double_t * pDouble_t
Definition: RooMath.h:35
double Double_t
Definition: RtypesCore.h:55
Double_t im() const
Definition: RooComplex.h:82
static void initFastCERF(Int_t, Double_t, Double_t, Int_t, Double_t, Double_t)
deprecated function
Definition: RooMath.h:220
Class RooComplex is a simple container class for complex values.
Definition: RooComplex.h:35
static Double_t ITPComplexErrFuncFastRe(const RooComplex &zz, Int_t)
deprecated function
Definition: RooMath.h:210
static Double_t erfc(Double_t x)
Definition: RooMath.h:188
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:584
static Double_t ITPComplexErrFuncFastIm(const RooComplex &zz, Int_t)
deprecated function
Definition: RooMath.h:213
static RooComplex ITPComplexErrFuncFast(const RooComplex &zz, Int_t)
deprecated function
Definition: RooMath.h:207
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 Double_t erf(Double_t x)
Definition: RooMath.h:185
static Double_t ComplexErrFuncFastRe(const RooComplex &zz)
deprecated function
Definition: RooMath.h:201