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