Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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