gamma.icc
Go to the documentation of this file.
#include "rheolef/compiler.h"
#include <cmath>
#include <cstdlib>
namespace rheolef {
template <class T> T my_gamma (const T& x) {
if (x == T(0)) return 1;
if (x == floor(x)) {
T value = 1;
for (size_t n = static_cast<int>(x)-1; n > 0; n--)
value *= T(1.*n);
// warning_macro ("gamma(" << x << ") = " << value);
return value;
}
static T sqrt_pi = sqrt(acos(T(-1.)));
if (x == -T(0.5)) return -2*sqrt_pi;
if (x-floor(x) == T(0.5)) {
T value = sqrt_pi;
for (size_t n = static_cast<int>(x); n != 0; n--)
value *= (n-0.5);
return value;
}
std::cerr << "gamma: " << x << " is not of integer of half order" << std::endl;
exit (1);
return 0;
}
} // namespace rheolef
rheolef::my_gamma
T my_gamma(const T &x)
Definition: gamma.icc:25
rheolef::value
rheolef::std value
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
T
Expr1::float_type T
Definition: field_expr.h:218