Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
35namespace ROOT {
36namespace Math {
37
38
39
40
41
42// (26.x.21.2) complementary error function
43
44double erfc(double x) {
45
46
47#ifdef USE_CEPHES
48 // use cephes implementation
50#else
51 return ::erfc(x);
52#endif
53
54}
55
56
57// (26.x.21.1) error function
58
59double erf(double x) {
60
61
62#ifdef USE_CEPHES
64#else
65 return ::erf(x);
66#endif
67
68
69}
70
71
72
73
74double lgamma(double z) {
75
76#ifdef USE_CEPHES
78#else
79 return ::lgamma(z);
80#endif
81
82}
83
84
85
86
87// (26.x.18) gamma function
88
89double tgamma(double x) {
90
91#ifdef USE_CEPHES
93#else
94 return ::tgamma(x);
95#endif
96
97}
98
99double inc_gamma( double a, double x) {
101}
102
103double inc_gamma_c( double a, double x) {
105}
106
107
108// [5.2.1.3] beta function
109// (26.x.19)
110
111double beta(double x, double y) {
112 return std::exp(lgamma(x)+lgamma(y)-lgamma(x+y));
113}
114
115double inc_beta( double x, double a, double b) {
117}
118
119// Sine integral
120// Translated from CERNLIB SININT (C336) by B. List 29.4.2010
121
122double 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
212double 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
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
#define PI
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
float * q
double inc_gamma_c(double a, double x)
Calculates the normalized (regularized) upper incomplete gamma function (upper integral)
double beta(double x, double y)
Calculates the beta function.
double inc_beta(double x, double a, double b)
Calculates the normalized (regularized) incomplete beta function.
double erfc(double x)
Complementary error function.
double sinint(double x)
Calculates the sine integral.
double tgamma(double x)
The gamma function is defined to be the extension of the factorial to real numbers.
double lgamma(double x)
Calculates the logarithm of the gamma function.
double cosint(double x)
Calculates the real part of the cosine integral Re(Ci).
double inc_gamma(double a, double x)
Calculates the normalized (regularized) lower incomplete gamma function (lower integral)
double erf(double x)
Error function encountered in integrating the normal distribution.
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
Namespace for new Math classes and functions.
double erfc(double a)
double erf(double x)
double incbet(double aa, double bb, double xx)
DESCRIPTION:
double igam(double a, double x)
double lgam(double x)
double igamc(double a, double x)
incomplete complementary gamma function igamc(a, x) = 1 - igam(a, x)
double gamma(double x)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.