29 #ifndef CDPL_MATH_SPECIALFUNCTIONS_HPP
30 #define CDPL_MATH_SPECIALFUNCTIONS_HPP
68 template <
typename T1,
typename T2>
69 T1
sign(
const T1& a,
const T2& b);
115 template <
typename T>
116 T gammaPSer(
const T& a,
const T& x)
118 using namespace CDPL;
119 using namespace Math;
123 return std::numeric_limits<T>::quiet_NaN();
132 for (
unsigned int i = 0; i < 100; i++) {
138 if (TypeTraits<T>::abs(del) < TypeTraits<T>::abs(
sum) * std::numeric_limits<T>::epsilon())
142 return std::numeric_limits<T>::quiet_NaN();
146 template <
typename T>
147 T gammaQContFrac(
const T& a,
const T& x)
149 using namespace CDPL;
150 using namespace Math;
152 static const T EPS = std::numeric_limits<T>::epsilon();
153 static const T FP_MIN = std::numeric_limits<T>::min() / EPS;
160 for (
unsigned int i = 0; i < 100; i++) {
165 if (TypeTraits<T>::abs(d) < FP_MIN)
170 if (TypeTraits<T>::abs(c) < FP_MIN)
177 if (TypeTraits<T>::abs(del - 1) <= EPS)
181 return std::numeric_limits<T>::quiet_NaN();
188 template <
typename T>
196 for (
unsigned int i = 2; i <= n; i++)
202 template <
typename T>
205 T abs_a = TypeTraits<T>::abs(a);
206 T abs_b = TypeTraits<T>::abs(b);
209 T tmp = abs_b / abs_a;
211 return (abs_a * TypeTraits<T>::sqrt(
T(1) + tmp * tmp));
217 T tmp = abs_a / abs_b;
219 return (abs_b * TypeTraits<T>::sqrt(
T(1) + tmp * tmp));
222 template <
typename T1,
typename T2>
225 return (b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a));
228 template <
typename T>
231 static const long double cof[6] = {
236 0.1208650973866179e-2,
237 -0.5395239384953e-5};
239 long double ser = 1.000000000190015;
242 long double tmp = x + 5.5;
244 tmp -= (x + 0.5) * std::log(tmp);
246 for (
unsigned int i = 0; i < 6; i++)
249 return T(-tmp + std::log(2.5066282746310005 * ser / x));
252 template <
typename T>
255 if (x <
T(0) || a <=
T(0))
256 return std::numeric_limits<T>::quiet_NaN();
259 return (
T(1) - Detail::gammaPSer(a, x));
261 return Detail::gammaQContFrac(a, x);
264 template <
typename T>
267 return (
T(1) / (
T(1) + std::pow(std::abs((x - c) / a),
T(2) * b)));
272 #endif // CDPL_MATH_SPECIALFUNCTIONS_HPP