Logo ROOT   6.12/07
Reference Guide
SpecFuncMathCore.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: Andras Zsenei & Lorenzo Moneta 06/2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 , LCG ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
11 #if defined(__sun) || defined(__sgi) || defined(_WIN32) || defined(_AIX)
12 #define NOT_HAVE_TGAMMA
13 #endif
14 
15 
16 #include "SpecFuncCephes.h"
17 
18 
19 #include <cmath>
20 #include <limits>
21 
22 #ifndef PI
23 #define PI 3.14159265358979323846264338328 /* pi */
24 #endif
25 
26 // use cephes for functions which are also in C99
27 #define USE_CEPHES
28 
29 // platforms not implemening C99
30 // #if defined(__sun) || defined(__sgi) || defined(_WIN32) || defined(_AIX)
31 // #define USE_CEPHES
32 // #endif
33 
34 
35 namespace ROOT {
36 namespace Math {
37 
38 
39 
40 
41 
42 // (26.x.21.2) complementary error function
43 
44 double erfc(double x) {
45 
46 
47 #ifdef USE_CEPHES
48  // use cephes implementation
49  return ROOT::Math::Cephes::erfc(x);
50 #else
51  return ::erfc(x);
52 #endif
53 
54 }
55 
56 
57 // (26.x.21.1) error function
58 
59 double erf(double x) {
60 
61 
62 #ifdef USE_CEPHES
63  return ROOT::Math::Cephes::erf(x);
64 #else
65  return ::erf(x);
66 #endif
67 
68 
69 }
70 
71 
72 
73 
74 double lgamma(double z) {
75 
76 #ifdef USE_CEPHES
77  return ROOT::Math::Cephes::lgam(z);
78 #else
79  return ::lgamma(z);
80 #endif
81 
82 }
83 
84 
85 
86 
87 // (26.x.18) gamma function
88 
89 double tgamma(double x) {
90 
91 #ifdef USE_CEPHES
92  return ROOT::Math::Cephes::gamma(x);
93 #else
94  return ::tgamma(x);
95 #endif
96 
97 }
98 
99 double inc_gamma( double a, double x) {
100  return ROOT::Math::Cephes::igam(a,x);
101 }
102 
103 double inc_gamma_c( double a, double x) {
104  return ROOT::Math::Cephes::igamc(a,x);
105 }
106 
107 
108 // [5.2.1.3] beta function
109 // (26.x.19)
110 
111 double beta(double x, double y) {
112  return std::exp(lgamma(x)+lgamma(y)-lgamma(x+y));
113 }
114 
115 double inc_beta( double x, double a, double b) {
116  return ROOT::Math::Cephes::incbet(a,b,x);
117 }
118 
119 // Sine integral
120 // Translated from CERNLIB SININT (C336) by B. List 29.4.2010
121 
122 double sinint(double x) {
123 
124  static const double z1 = 1, r8 = z1/8;
125 
126  static const double pih = PI/2;
127 
128  static const double s[16] = {
129  +1.95222097595307108, -0.68840423212571544,
130  +0.45518551322558484, -0.18045712368387785,
131  +0.04104221337585924, -0.00595861695558885,
132  +0.00060014274141443, -0.00004447083291075,
133  +0.00000253007823075, -0.00000011413075930,
134  +0.00000000418578394, -0.00000000012734706,
135  +0.00000000000326736, -0.00000000000007168,
136  +0.00000000000000136, -0.00000000000000002};
137 
138  static const double p[29] = {
139  +0.96074783975203596, -0.03711389621239806,
140  +0.00194143988899190, -0.00017165988425147,
141  +0.00002112637753231, -0.00000327163256712,
142  +0.00000060069211615, -0.00000012586794403,
143  +0.00000002932563458, -0.00000000745695921,
144  +0.00000000204105478, -0.00000000059502230,
145  +0.00000000018322967, -0.00000000005920506,
146  +0.00000000001996517, -0.00000000000699511,
147  +0.00000000000253686, -0.00000000000094929,
148  +0.00000000000036552, -0.00000000000014449,
149  +0.00000000000005851, -0.00000000000002423,
150  +0.00000000000001025, -0.00000000000000442,
151  +0.00000000000000194, -0.00000000000000087,
152  +0.00000000000000039, -0.00000000000000018,
153  +0.00000000000000008};
154 
155  static const double q[25] = {
156  +0.98604065696238260, -0.01347173820829521,
157  +0.00045329284116523, -0.00003067288651655,
158  +0.00000313199197601, -0.00000042110196496,
159  +0.00000006907244830, -0.00000001318321290,
160  +0.00000000283697433, -0.00000000067329234,
161  +0.00000000017339687, -0.00000000004786939,
162  +0.00000000001403235, -0.00000000000433496,
163  +0.00000000000140273, -0.00000000000047306,
164  +0.00000000000016558, -0.00000000000005994,
165  +0.00000000000002237, -0.00000000000000859,
166  +0.00000000000000338, -0.00000000000000136,
167  +0.00000000000000056, -0.00000000000000024,
168  +0.00000000000000010};
169 
170  double h;
171  if (std::abs(x) <= 8) {
172  double y = r8*x;
173  h = 2*y*y-1;
174  double alfa = h+h;
175  double b0 = 0;
176  double b1 = 0;
177  double b2 = 0;
178  for (int i = 15; i >= 0; --i) {
179  b0 = s[i]+alfa*b1-b2;
180  b2 = b1;
181  b1 = b0;
182  }
183  h = y*(b0-b2);
184  } else {
185  double r = 1/x;
186  h = 128*r*r-1;
187  double alfa = h+h;
188  double b0 = 0;
189  double b1 = 0;
190  double b2 = 0;
191  for (int i = 28; i >= 0; --i) {
192  b0 = p[i]+alfa*b1-b2;
193  b2 = b1;
194  b1 = b0;
195  }
196  double pp = b0-h*b2;
197  b1 = 0;
198  b2 = 0;
199  for (int i = 24; i >= 0; --i) {
200  b0 = q[i]+alfa*b1-b2;
201  b2 = b1;
202  b1 = b0;
203  }
204  h = (x > 0 ? pih : -pih)-r*(r*pp*std::sin(x)+(b0-h*b2)*std::cos(x));
205  }
206  return h;
207 }
208 
209 // Real part of the cosine integral
210 // Translated from CERNLIB COSINT (C336) by B. List 29.4.2010
211 
212 double cosint(double x) {
213 
214  static const double z1 = 1, r32 = z1/32;
215 
216  static const double ce = 0.57721566490153286;
217 
218  static const double c[16] = {
219  +1.94054914648355493, +0.94134091328652134,
220  -0.57984503429299276, +0.30915720111592713,
221  -0.09161017922077134, +0.01644374075154625,
222  -0.00197130919521641, +0.00016925388508350,
223  -0.00001093932957311, +0.00000055223857484,
224  -0.00000002239949331, +0.00000000074653325,
225  -0.00000000002081833, +0.00000000000049312,
226  -0.00000000000001005, +0.00000000000000018};
227 
228  static const double p[29] = {
229  +0.96074783975203596, -0.03711389621239806,
230  +0.00194143988899190, -0.00017165988425147,
231  +0.00002112637753231, -0.00000327163256712,
232  +0.00000060069211615, -0.00000012586794403,
233  +0.00000002932563458, -0.00000000745695921,
234  +0.00000000204105478, -0.00000000059502230,
235  +0.00000000018322967, -0.00000000005920506,
236  +0.00000000001996517, -0.00000000000699511,
237  +0.00000000000253686, -0.00000000000094929,
238  +0.00000000000036552, -0.00000000000014449,
239  +0.00000000000005851, -0.00000000000002423,
240  +0.00000000000001025, -0.00000000000000442,
241  +0.00000000000000194, -0.00000000000000087,
242  +0.00000000000000039, -0.00000000000000018,
243  +0.00000000000000008};
244 
245  static const double q[25] = {
246  +0.98604065696238260, -0.01347173820829521,
247  +0.00045329284116523, -0.00003067288651655,
248  +0.00000313199197601, -0.00000042110196496,
249  +0.00000006907244830, -0.00000001318321290,
250  +0.00000000283697433, -0.00000000067329234,
251  +0.00000000017339687, -0.00000000004786939,
252  +0.00000000001403235, -0.00000000000433496,
253  +0.00000000000140273, -0.00000000000047306,
254  +0.00000000000016558, -0.00000000000005994,
255  +0.00000000000002237, -0.00000000000000859,
256  +0.00000000000000338, -0.00000000000000136,
257  +0.00000000000000056, -0.00000000000000024,
258  +0.00000000000000010};
259 
260  double h = 0;
261  if(x == 0) {
262  h = - std::numeric_limits<double>::infinity();
263  } else if (std::abs(x) <= 8) {
264  h = r32*x*x-1;
265  double alfa = h+h;
266  double b0 = 0;
267  double b1 = 0;
268  double b2 = 0;
269  for (int i = 15; i >= 0; --i) {
270  b0 = c[i]+alfa*b1-b2;
271  b2 = b1;
272  b1 = b0;
273  }
274  h = ce+std::log(std::abs(x))-b0+h*b2;
275  } else {
276  double r = 1/x;
277  h = 128*r*r-1;
278  double alfa = h+h;
279  double b0 = 0;
280  double b1 = 0;
281  double b2 = 0;
282  for (int i = 28; i >= 0; --i) {
283  b0 = p[i]+alfa*b1-b2;
284  b2 = b1;
285  b1 = b0;
286  }
287  double pp = b0-h*b2;
288  b1 = 0;
289  b2 = 0;
290  for (int i = 24; i >= 0; --i) {
291  b0 = q[i]+alfa*b1-b2;
292  b2 = b1;
293  b1 = b0;
294  }
295  h = r*((b0-h*b2)*std::sin(x)-r*pp*std::cos(x));
296  }
297  return h;
298 }
299 
300 
301 
302 
303 } // namespace Math
304 } // namespace ROOT
305 
306 
307 
308 
309 
double erf(double x)
Error function encountered in integrating the normal distribution.
double igam(double a, double x)
double cosint(double x)
Calculates the real part of the cosine integral Re(Ci).
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
double sinint(double x)
Calculates the sine integral.
TH1 * h
Definition: legend2.C:5
double lgam(double x)
double inc_beta(double x, double a, double b)
Calculates the normalized (regularized) incomplete beta function.
#define PI
double tgamma(double x)
The gamma function is defined to be the extension of the factorial to real numbers.
double cos(double)
double beta(double x, double y)
Calculates the beta function.
double erfc(double x)
Complementary error function.
double inc_gamma_c(double a, double x)
Calculates the normalized (regularized) upper incomplete gamma function (upper integral) ...
Double_t x[n]
Definition: legend1.C:17
double sin(double)
double gamma(double x)
ROOT::R::TRInterface & r
Definition: Object.C:4
auto * a
Definition: textangle.C:12
double incbet(double aa, double bb, double xx)
DESCRIPTION:
double erf(double x)
double igamc(double a, double x)
incomplete complementary gamma function igamc(a, x) = 1 - igam(a, x)
Double_t y[n]
Definition: legend1.C:17
static constexpr double s
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
double erfc(double a)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
double lgamma(double x)
Calculates the logarithm of the gamma function.
double exp(double)
float * q
Definition: THbookFile.cxx:87
double inc_gamma(double a, double x)
Calculates the normalized (regularized) lower incomplete gamma function (lower integral) ...
double log(double)