Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JMath/JMathSupportkit.hh
Go to the documentation of this file.
1#ifndef __JMATHSUPPORTKIT__
2#define __JMATHSUPPORTKIT__
3
4#include <limits>
5#include <cmath>
6
7#include "JMath/JConstants.hh"
8#include "JLang/JException.hh"
9
10
11/**
12 * \file
13 *
14 * Auxiliary methods for mathematics.
15 * \author mdejong
16 */
17
18namespace JMATH {}
19namespace JPP { using namespace JMATH; }
20
21namespace JMATH {
22
24
25
26 /**
27 * Gauss function (normalised to 1 at x = 0).
28 *
29 * \param x x
30 * \param sigma sigma
31 * \return function value
32 */
33 inline double gauss(const double x, const double sigma)
34 {
35 const double u = x / sigma;
36
37 if (fabs(u) < 20.0)
38 return exp(-0.5*u*u);
39 else
40 return 0.0;
41 }
42
43
44 /**
45 * Gauss function (normalised to 1 at x = x0).
46 *
47 * \param x x
48 * \param x0 central value
49 * \param sigma sigma
50 * \return function value
51 */
52 inline double gauss(const double x, const double x0, const double sigma)
53 {
54 return gauss(x - x0, sigma);
55 }
56
57
58 /**
59 * Normalised Gauss function.
60 *
61 * \param x x
62 * \param sigma sigma
63 * \return function value
64 */
65 inline double Gauss(const double x, const double sigma)
66 {
67 return gauss(x, sigma) / sqrt(2.0*PI) / sigma;
68 }
69
70
71 /**
72 * Normalised Gauss function.
73 *
74 * \param x x
75 * \param x0 central value
76 * \param sigma sigma
77 * \return function value
78 */
79 inline double Gauss(const double x, const double x0, const double sigma)
80 {
81 return Gauss(x - x0, sigma);
82 }
83
84
85 /**
86 * Incomplete gamma function.
87 *
88 * Source code is taken from reference:
89 * Numerical Recipes in C++, W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery,
90 * Cambridge University Press.
91 *
92 * \param a a
93 * \param x x
94 * \return function value
95 */
96 inline double Gamma(const double a, const double x)
97 {
98 using namespace std;
99
100 const int max = 1000000;
101
102 if (x < 0.0) { THROW(JValueOutOfRange, "x < 0 " << x); }
103 if (a <= 0.0) { THROW(JValueOutOfRange, "a <= 0 " << a); }
104
105 const double gln = lgamma(a);
106
107 if (x < a + 1.0) {
108
109 if (x < 0.0) {
110 THROW(JValueOutOfRange, "x <= 0 " << x);
111 }
112
113 if (x == 0.0) {
114 return 0.0;
115 }
116
117 double ap = a;
118 double sum = 1.0 / a;
119 double del = sum;
120
121 for (int i = 0; i != max; ++i) {
122
123 ap += 1.0;
124 del *= x/ap;
125 sum += del;
126
127 if (fabs(del) < fabs(sum)*numeric_limits<double>::epsilon()) {
128 return sum*exp(-x + a*log(x) - gln);
129 }
130 }
131
132 THROW(JValueOutOfRange, "i == " << max);
133
134 } else {
135
136 const double FPMIN = numeric_limits<double>::min() / numeric_limits<double>::epsilon();
137
138 double b = x + 1.0 - a;
139 double c = 1.0 / FPMIN;
140 double d = 1.0 / b;
141 double h = d;
142
143 for (int i = 1; i != max; ++i) {
144
145 const double an = -i * (i-a);
146
147 b += 2.0;
148 d = an*d + b;
149
150 if (fabs(d) < FPMIN) {
151 d = FPMIN;
152 }
153
154 c = b + an/c;
155
156 if (fabs(c) < FPMIN) {
157 c = FPMIN;
158 }
159
160 d = 1.0/d;
161
162 const double del = d*c;
163
164 h *= del;
165
166 if (fabs(del - 1.0) < numeric_limits<double>::epsilon()) {
167 return 1.0 - exp(-x + a*log(x) - gln) * h;
168 }
169 }
170
171 THROW(JValueOutOfRange, "i == " << max);
172 }
173 }
174
175
176 /**
177 * Legendre polynome.
178 *
179 * \param n degree
180 * \param x x
181 * \return function value
182 */
183 inline double legendre(const size_t n, const double x)
184 {
185 switch (n) {
186
187 case 0:
188 return 1.0;
189
190 case 1:
191 return x;
192
193 default:
194 {
195 double p0;
196 double p1 = 1.0;
197 double p2 = x;
198
199 for (size_t i = 2; i <= n; ++i) {
200 p0 = p1;
201 p1 = p2;
202 p2 = ((2*i-1) * x*p1 - (i-1) * p0) / i;
203 }
204
205 return p2;
206 }
207 }
208 }
209
210
211 /**
212 * Binomial function.
213 *
214 * \param n n
215 * \param k k
216 * \return function value
217 */
218 inline double binomial(const size_t n, const size_t k)
219 {
220 if (n == 0 || n < k) {
221 return 0.0;
222 }
223
224 if (k == 0 || n == k) {
225 return 1.0;
226 }
227
228 const int k1 = std::min(k, n - k);
229 const int k2 = n - k1;
230
231 double value = k2 + 1;
232
233 for (int i = k1; i != 1; --i) {
234 value *= (double) (k2 + i) / (double) i;
235 }
236
237 return value;
238 }
239
240
241 /**
242 * Poisson probability density distribition.
243 *
244 * \param n number of occurences
245 * \param mu expectation value
246 * \return probability
247 */
248 inline double poisson(const size_t n, const double mu)
249 {
250 using namespace std;
251
252 if (mu > 0.0) {
253
254 if (n > 0)
255 return exp(n*log(mu) - lgamma(n+1) - mu);
256 else
257 return exp(-mu);
258 } else if (mu == 0.0) {
259
260 return (n == 0 ? 1.0 : 0.0);
261 }
262
263 THROW(JValueOutOfRange, "mu <= 0 " << mu);
264 }
265
266
267 /**
268 * Poisson cumulative density distribition.
269 *
270 * \param n number of occurences
271 * \param mu expectation value
272 * \return probability
273 */
274 inline double Poisson(const size_t n, const double mu)
275 {
276 return 1.0 - Gamma(n + 1, mu);
277 }
278}
279
280#endif
TPaveText * p1
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Mathematical constants.
Exception for accessing a value in a collection that is outside of its range.
Auxiliary classes and methods for mathematical operations.
Definition JEigen3D.hh:88
double Gamma(const double a, const double x)
Incomplete gamma function.
double Gauss(const double x, const double sigma)
Normalised Gauss function.
double poisson(const size_t n, const double mu)
Poisson probability density distribition.
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
double Poisson(const size_t n, const double mu)
Poisson cumulative density distribition.
double legendre(const size_t n, const double x)
Legendre polynome.
double binomial(const size_t n, const size_t k)
Binomial function.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).