ROOT   6.14/05 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.
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)
#define h(i)
Definition: RSha256.hxx:106
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
#define c(i)
Definition: RSha256.hxx:101
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)