Jpp  15.0.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
18 namespace JMATH {}
19 namespace JPP { using namespace JMATH; }
20 
21 namespace 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) < 10.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  * \param a a
89  * \param x x
90  * \return function value
91  */
92  inline double Gamma(const double a, const double x)
93  {
94  using namespace std;
95 
96  const int max = 100;
97 
98  if (x < 0.0) { THROW(JValueOutOfRange, "x < 0 " << x); }
99  if (a <= 0.0) { THROW(JValueOutOfRange, "a <= 0 " << a); }
100 
101  const double gln = lgamma(a);
102 
103  if (x < a + 1.0) {
104 
105  if (x <= 0.0) {
106  THROW(JValueOutOfRange, "x <= 0 " << x);
107  }
108 
109  double ap = a;
110  double sum = 1.0 /a;
111  double del = sum;
112 
113  for (int i = 1; i != max; ++i) {
114 
115  ap += 1.0;
116  del *= x/ap;
117  sum += del;
118 
119  if (fabs(del) < fabs(sum)*numeric_limits<double>::epsilon()) {
120  return sum*exp(-x + a*log(x) - gln);
121  }
122  }
123 
124  } else {
125 
126  double b = x + 1.0 - a;
127  double c = numeric_limits<double>::epsilon() / numeric_limits<double>::min();
128  double d = 1.0 / b;
129  double h = d;
130 
131  for (int i = 1; i != max; ++i) {
132 
133  const double an = -i * (i-a);
134 
135  b += 2.0;
136  d = an*d + b;
137 
138  if (fabs(d) < numeric_limits<double>::min()) {
139  d = numeric_limits<double>::min();
140  }
141 
142  c = b + an/c;
143 
144  if (fabs(c) < numeric_limits<double>::min()) {
145  c = numeric_limits<double>::min();
146  }
147 
148  d = 1.0/d;
149 
150  const double del = d*c;
151 
152  h *= del;
153 
154  if (fabs(del - 1.0) < numeric_limits<double>::epsilon()) {
155  return 1.0 - exp(-x + a*log(x) - gln) * h;
156  }
157  }
158 
159  THROW(JValueOutOfRange, "a " << a);
160  }
161 
162  return 0.0;
163  }
164 
165 
166  /**
167  * Legendre polynome.
168  *
169  * \param n degree
170  * \param x x
171  * \return function value
172  */
173  inline double legendre(const unsigned int n, const double x)
174  {
175  switch (n) {
176 
177  case 0:
178  return 1.0;
179 
180  case 1:
181  return x;
182 
183  default:
184  {
185  double p0;
186  double p1 = 1.0;
187  double p2 = x;
188 
189  for (unsigned int i = 2; i <= n; ++i) {
190  p0 = p1;
191  p1 = p2;
192  p2 = ((2*i-1) * x*p1 - (i-1) * p0) / i;
193  }
194 
195  return p2;
196  }
197  }
198  }
199 
200 
201  /**
202  * Binomial function.
203  *
204  * \param n n
205  * \param k k
206  * \return function value
207  */
208  inline double binomial(const int n, const int k)
209  {
210  if (k == 0 || n == k) {
211  return 1.0;
212  }
213 
214  if (n <= 0 || k < 0 || n < k) {
215  return 0.0;
216  }
217 
218  const int k1 = std::min(k, n - k);
219  const int k2 = n - k1;
220 
221  double value = k2 + 1;
222 
223  for (int i = k1; i != 1; --i) {
224  value *= (double) (k2 + i) / (double) i;
225  }
226 
227  return value;
228  }
229 }
230 
231 #endif
Exceptions.
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
TPaveText * p1
double binomial(const int n, const int k)
Binomial function.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
double Gamma(const double a, const double x)
Incomplete gamma function.
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` typeset -Z 4 STRING JOpera1D -f hydrophone.root
const int n
Definition: JPolint.hh:660
Mathematical constants.
static const double PI
Mathematical constants.
double legendre(const unsigned int n, const double x)
Legendre polynome.
p2
Definition: module-Z:fit.sh:74
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
then JCalibrateToT a
Definition: JTuneHV.sh:116
double gauss(const double x, const double sigma)
Gauss function (normalised to 1 at x = 0).
double Gauss(const double x, const double sigma)
Normalised Gauss function.
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
double u[N+1]
Definition: JPolint.hh:739