Chemical Data Processing Library C++ API - Version 1.1.1
SpecialFunctions.hpp
Go to the documentation of this file.
1 /*
2  * SpecialFunctions.hpp
3  *
4  * This file is part of the Chemical Data Processing Toolkit
5  *
6  * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
7  *
8  * This library is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU Lesser General Public
10  * License as published by the Free Software Foundation; either
11  * version 2 of the License, or (at your option) any later version.
12  *
13  * This library is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  * Lesser General Public License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public License
19  * along with this library; see the file COPYING. If not, write to
20  * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
21  * Boston, MA 02111-1307, USA.
22  */
23 
29 #ifndef CDPL_MATH_SPECIALFUNCTIONS_HPP
30 #define CDPL_MATH_SPECIALFUNCTIONS_HPP
31 
32 #include <cstddef>
33 #include <cmath>
34 #include <limits>
35 
36 #include "CDPL/Math/TypeTraits.hpp"
37 
38 
39 namespace CDPL
40 {
41 
42  namespace Math
43  {
44 
50  template <typename T>
51  T factorial(unsigned int n);
52 
59  template <typename T>
60  T pythag(const T& a, const T& b);
61 
68  template <typename T1, typename T2>
69  T1 sign(const T1& a, const T2& b);
70 
76  template <typename T>
77  T lnGamma(const T& z);
78 
85  template <typename T>
86  T gammaQ(const T& a, const T& x);
87 
96  template <typename T>
97  T generalizedBell(const T& x, const T& a, const T& b, const T& c);
98  } // namespace Math
99 } // namespace CDPL
100 
101 
102 // Implementation
103 // \cond DOC_IMPL_DETAILS
104 
105 namespace CDPL
106 {
107 
108  namespace Math
109  {
110 
111  namespace Detail
112  {
113 
114  // Computes the incomplete gamma function P(a, x) evaluated by its series representation.
115  template <typename T>
116  T gammaPSer(const T& a, const T& x)
117  {
118  using namespace CDPL;
119  using namespace Math;
120 
121  if (x <= 0) {
122  if (x < 0)
123  return std::numeric_limits<T>::quiet_NaN();
124 
125  return 0;
126  }
127 
128  T ap = a;
129  T del = 1 / a;
130  T sum = del;
131 
132  for (unsigned int i = 0; i < 100; i++) {
133  ++ap;
134 
135  del *= x / ap;
136  sum += del;
137 
138  if (TypeTraits<T>::abs(del) < TypeTraits<T>::abs(sum) * std::numeric_limits<T>::epsilon())
139  return (sum * std::exp(-x + a * log(x) - CDPL::Math::lnGamma(a)));
140  }
141 
142  return std::numeric_limits<T>::quiet_NaN();
143  }
144 
145  // Computes the incomplete gamma function Q(a, x) evaluated by its continued fraction representation.
146  template <typename T>
147  T gammaQContFrac(const T& a, const T& x)
148  {
149  using namespace CDPL;
150  using namespace Math;
151 
152  static const T EPS = std::numeric_limits<T>::epsilon();
153  static const T FP_MIN = std::numeric_limits<T>::min() / EPS;
154 
155  T b = x + 1 - a;
156  T c = 1 / FP_MIN;
157  T d = 1 / b;
158  T h = d;
159 
160  for (unsigned int i = 0; i < 100; i++) {
161  T an = -i * (i - a);
162  b += 2;
163  d = an * d + b;
164 
165  if (TypeTraits<T>::abs(d) < FP_MIN)
166  d = FP_MIN;
167 
168  c = b + an / c;
169 
170  if (TypeTraits<T>::abs(c) < FP_MIN)
171  c = FP_MIN;
172 
173  d = 1 / d;
174  T del = d * c;
175  h *= del;
176 
177  if (TypeTraits<T>::abs(del - 1) <= EPS)
178  return (std::exp(-x + a * log(x) - CDPL::Math::lnGamma(a)) * h);
179  }
180 
181  return std::numeric_limits<T>::quiet_NaN();
182  }
183  } // namespace Detail
184  } // namespace Math
185 } // namespace CDPL
186 
187 
188 template <typename T>
189 T CDPL::Math::factorial(unsigned int n)
190 {
191  T f(1);
192 
193  if (n == 0)
194  return f;
195 
196  for (unsigned int i = 2; i <= n; i++)
197  f *= T(i);
198 
199  return f;
200 }
201 
202 template <typename T>
203 T CDPL::Math::pythag(const T& a, const T& b)
204 {
205  T abs_a = TypeTraits<T>::abs(a);
206  T abs_b = TypeTraits<T>::abs(b);
207 
208  if (abs_a > abs_b) {
209  T tmp = abs_b / abs_a;
210 
211  return (abs_a * TypeTraits<T>::sqrt(T(1) + tmp * tmp));
212  }
213 
214  if (abs_b == 0)
215  return 0;
216 
217  T tmp = abs_a / abs_b;
218 
219  return (abs_b * TypeTraits<T>::sqrt(T(1) + tmp * tmp));
220 }
221 
222 template <typename T1, typename T2>
223 T1 CDPL::Math::sign(const T1& a, const T2& b)
224 {
225  return (b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));
226 }
227 
228 template <typename T>
229 T CDPL::Math::lnGamma(const T& xx)
230 {
231  static const long double cof[6] = {
232  76.18009172947146,
233  -86.50532032941677,
234  24.01409824083091,
235  -1.231739572450155,
236  0.1208650973866179e-2,
237  -0.5395239384953e-5};
238 
239  long double ser = 1.000000000190015;
240  long double x = xx;
241  long double y = xx;
242  long double tmp = x + 5.5;
243 
244  tmp -= (x + 0.5) * std::log(tmp);
245 
246  for (unsigned int i = 0; i < 6; i++)
247  ser += cof[i] / ++y;
248 
249  return T(-tmp + std::log(2.5066282746310005 * ser / x));
250 }
251 
252 template <typename T>
253 T CDPL::Math::gammaQ(const T& a, const T& x)
254 {
255  if (x < T(0) || a <= T(0))
256  return std::numeric_limits<T>::quiet_NaN();
257 
258  if (x < (a + T(1)))
259  return (T(1) - Detail::gammaPSer(a, x));
260 
261  return Detail::gammaQContFrac(a, x);
262 }
263 
264 template <typename T>
265 T CDPL::Math::generalizedBell(const T& x, const T& a, const T& b, const T& c)
266 {
267  return (T(1) / (T(1) + std::pow(std::abs((x - c) / a), T(2) * b)));
268 }
269 
270 // \endcond
271 
272 #endif // CDPL_MATH_SPECIALFUNCTIONS_HPP
CDPL::Math::pythag
T pythag(const T &a, const T &b)
Computes without destructive underflow or overflow.
CDPL::Math::generalizedBell
T generalizedBell(const T &x, const T &a, const T &b, const T &c)
Computes the generalized bell function at x.
CDPL::Math::gammaQ
T gammaQ(const T &a, const T &x)
Computes the incomplete gamma function (see [NRIC] for details).
CDPL::Math::sum
GridElementSum< E >::ResultType sum(const GridExpression< E > &e)
Definition: GridExpression.hpp:416
TypeTraits.hpp
Definition of type traits.
CDPL::Chem::AtomType::T
const unsigned int T
Specifies Hydrogen (Tritium).
Definition: AtomType.hpp:67
CDPL
The namespace of the Chemical Data Processing Library.
CDPL::Math::factorial
T factorial(unsigned int n)
Computes the factorial of the non-negative integer n.
CDPL::Math::sign
T1 sign(const T1 &a, const T2 &b)
Returns the magnitude of parameter a times the sign of parameter b.
CDPL::Math::lnGamma
T lnGamma(const T &z)
Computes for .