Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
gamma.icc
Go to the documentation of this file.
1
21#include "rheolef/compiler.h"
22#include <cmath>
23#include <cstdlib>
24namespace rheolef {
25template <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
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
T my_gamma(const T &x)
Definition gamma.icc:25