Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files 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 
25 typedef Double_t* pDouble_t;
26 
27 class RooMath {
28 public:
29 
30  virtual ~RooMath() {};
31 
32  /** @brief evaluate Faddeeva function for complex argument
33  *
34  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
35  * @date 2013-02-21
36  *
37  * Calculate the value of the Faddeeva function @f$w(z) = \exp(-z^2)
38  * \mathrm{erfc}(-i z)@f$.
39  *
40  * The method described in
41  *
42  * S.M. Abrarov, B.M. Quine: "Efficient algotithmic implementation of
43  * Voigt/complex error function based on exponential series approximation"
44  * published in Applied Mathematics and Computation 218 (2011) 1894-1902
45  * doi:10.1016/j.amc.2011.06.072
46  *
47  * is used. At the heart of the method (equation (14) of the paper) is the
48  * following Fourier series based approximation:
49  *
50  * @f[ w(z) \approx \frac{i}{2\sqrt{\pi}}\left(
51  * \sum^N_{n=0} a_n \tau_m\left(
52  * \frac{1-e^{i(n\pi+\tau_m z)}}{n\pi + \tau_m z} -
53  * \frac{1-e^{i(-n\pi+\tau_m z)}}{n\pi - \tau_m z}
54  * \right) - a_0 \frac{1-e^{i \tau_m z}}{z}
55  * \right) @f]
56  *
57  * The coefficients @f$a_b@f$ are given by:
58  *
59  * @f[ a_n=\frac{2\sqrt{\pi}}{\tau_m}
60  * \exp\left(-\frac{n^2\pi^2}{\tau_m^2}\right) @f]
61  *
62  * To achieve machine accuracy in double precision floating point arithmetic
63  * for most of the upper half of the complex plane, chose @f$t_m=12@f$ and
64  * @f$N=23@f$ as is done in the paper.
65  *
66  * There are two complications: For Im(z) negative, the exponent in the
67  * equation above becomes so large that the roundoff in the rest of the
68  * calculation is amplified enough that the result cannot be trusted.
69  * Therefore, for Im(z) < 0, the symmetry of the erfc function under the
70  * transformation z --> -z is used to avoid accuracy issues for Im(z) < 0 by
71  * formulating the problem such that the calculation can be done for Im(z) > 0
72  * where the accuracy of the method is fine, and some postprocessing then
73  * yields the desired final result.
74  *
75  * Second, the denominators in the equation above become singular at
76  * @f$z = n * pi / 12@f$ (for 0 <= n < 24). In a tiny disc around these
77  * points, Taylor expansions are used to overcome that difficulty.
78  *
79  * This routine precomputes everything it can, and tries to write out complex
80  * operations to minimise subroutine calls, e.g. for the multiplication of
81  * complex numbers.
82  *
83  * In the square -8 <= Re(z) <= 8, -8 <= Im(z) <= 8, the routine is accurate
84  * to better than 4e-13 relative, the average relative error is better than
85  * 7e-16. On a modern x86_64 machine, the routine is roughly three times as
86  * fast than the old CERNLIB implementation and offers better accuracy.
87  *
88  * For large @f$|z|@f$, the familiar continued fraction approximation
89  *
90  * @f[ w(z)=\frac{-iz/\sqrt{\pi}}{-z^2+\frac{1/2}{1+\frac{2/2}{-z^2 +
91  * \frac{3/2}{1+\frac{4/2}{-z^2+\frac{5/2}{1+\frac{6/2}{-z^2+\frac{7/2
92  * }{1+\frac{8/2}{-z^2+\frac{9/2}{1+\ldots}}}}}}}}}} @f]
93  *
94  * is used, truncated at the ellipsis ("...") in the formula; for @f$|z| >
95  * 12@f$, @f$Im(z)>0@f$ it will give full double precision at a smaller
96  * computational cost than the method described above. (For @f$|z|>12@f$,
97  * @f$Im(z)<0@f$, the symmetry property @f$w(x-iy)=2e^{-(x+iy)^2-w(x+iy)}@f$
98  * is used.
99  */
100  static std::complex<double> faddeeva(std::complex<double> z);
101  /** @brief evaluate Faddeeva function for complex argument (fast version)
102  *
103  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
104  * @date 2013-02-21
105  *
106  * Calculate the value of the Faddeeva function @f$w(z) = \exp(-z^2)
107  * \mathrm{erfc}(-i z)@f$.
108  *
109  * This is the "fast" version of the faddeeva routine above. Fast means that
110  * is takes roughly half the amount of CPU of the slow version of the
111  * routine, but is a little less accurate.
112  *
113  * To be fast, chose @f$t_m=8@f$ and @f$N=11@f$ which should give accuracies
114  * around 1e-7.
115  *
116  * In the square -8 <= Re(z) <= 8, -8 <= Im(z) <= 8, the routine is accurate
117  * to better than 4e-7 relative, the average relative error is better than
118  * 5e-9. On a modern x86_64 machine, the routine is roughly five times as
119  * fast than the old CERNLIB implementation, or about 30% faster than the
120  * interpolation/lookup table based fast method used previously in RooFit,
121  * and offers better accuracy than the latter (the relative error is roughly
122  * a factor 280 smaller than the old interpolation/table lookup routine).
123  *
124  * For large @f$|z|@f$, the familiar continued fraction approximation
125  *
126  * @f[ w(z)=\frac{-iz/\sqrt{\pi}}{-z^2+\frac{1/2}{1+\frac{2/2}{-z^2 +
127  * \frac{3/2}{1+\ldots}}}} @f]
128  *
129  * is used, truncated at the ellipsis ("...") in the formula; for @f$|z| >
130  * 8@f$, @f$Im(z)>0@f$ it will give full float precision at a smaller
131  * computational cost than the method described above. (For @f$|z|>8@f$,
132  * @f$Im(z)<0@f$, the symmetry property @f$w(x-iy)=2e^{-(x+iy)^2-w(x+iy)}@f$
133  * is used.
134  */
135  static std::complex<double> faddeeva_fast(std::complex<double> z);
136 
137  /** @brief complex erf function
138  *
139  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
140  * @date 2013-02-21
141  *
142  * Calculate erf(z) for complex z.
143  */
144  static std::complex<double> erf(const std::complex<double> z);
145 
146  /** @brief complex erf function (fast version)
147  *
148  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
149  * @date 2013-02-21
150  *
151  * Calculate erf(z) for complex z. Use the code in faddeeva_fast to save some time.
152  */
153  static std::complex<double> erf_fast(const std::complex<double> z);
154  /** @brief complex erfc function
155  *
156  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
157  * @date 2013-02-21
158  *
159  * Calculate erfc(z) for complex z.
160  */
161  static std::complex<double> erfc(const std::complex<double> z);
162  /** @brief complex erfc function (fast version)
163  *
164  * @author Manuel Schiller <manuel.schiller@nikhef.nl>
165  * @date 2013-02-21
166  *
167  * Calculate erfc(z) for complex z. Use the code in faddeeva_fast to save some time.
168  */
169  static std::complex<double> erfc_fast(const std::complex<double> z);
170 
171  // 1-D nth order polynomial interpolation routines
172  static Double_t interpolate(Double_t yArr[],Int_t nOrder, Double_t x) ;
173  static Double_t interpolate(Double_t xa[], Double_t ya[], Int_t n, Double_t x) ;
174 
175  static inline Double_t erf(Double_t x)
176  { return TMath::Erf(x); }
177 
178  static inline Double_t erfc(Double_t x)
179  { return TMath::Erfc(x); }
180 
181 };
182 
183 #endif