Rheolef  7.1
an efficient C++ finite element environment
gamma.icc
Go to the documentation of this file.
1 #include "rheolef/compiler.h"
22 #include <cmath>
23 #include <cstdlib>
24 namespace rheolef {
25 template <class T> T my_gamma (const T& x) {
26  if (x == T(0)) return 1;
27  if (x == floor(x)) {
28  T value = 1;
29  for (size_t n = static_cast<int>(x)-1; n > 0; n--)
30  value *= T(1.*n);
31  // warning_macro ("gamma(" << x << ") = " << value);
32  return value;
33  }
34  static T sqrt_pi = sqrt(acos(T(-1.)));
35  if (x == -T(0.5)) return -2*sqrt_pi;
36  if (x-floor(x) == T(0.5)) {
37  T value = sqrt_pi;
38  for (size_t n = static_cast<int>(x); n != 0; n--)
39  value *= (n-0.5);
40  return value;
41  }
42  std::cerr << "gamma: " << x << " is not of integer of half order" << std::endl;
43  exit (1);
44  return 0;
45 }
46 } // namespace rheolef
rheolef::std value
Expr1::float_type T
Definition: field_expr.h:261
This file is part of Rheolef.
T my_gamma(const T &x)
Definition: gamma.icc:25