30#include "gsl/gsl_sf_gamma.h"
31#include "gsl/gsl_sf_pow_int.h"
41 std::vector< std::vector<int> > qt =
getQIndex(
n);
42 std::vector<double> g;
44 for (
size_t i = 0; i < qt.size(); ++i) {
46 for (
size_t j = 1; j < qt[i].size(); ++j) {
47 if (j > g.size()) g.push_back(
gN(x, j-1));
48 ei *= gsl_sf_pow_int(g[j-1], qt[i][j]);
61 for (
unsigned int i =
n; i <
_a.size(); i++)
62 hn +=
_a[i]/gsl_sf_pow_int(
_lambda, i)*gsl_sf_pow_int(x, i-
n)*
63 gsl_sf_fact(i)/gsl_sf_fact(i-
n);
70 if (
n == 0)
return 1+
exp(
hN(x, 0));
71 std::vector<double> hn(
n+1);
72 for (
int i = 0; i <=
n; i++) hn[i] =
hN(x, i);
73 double exp_h0 =
exp(hn[0]);
75 for (
size_t i = 0; i <
_h[
n].size(); ++i) {
76 double gnj =
_h[
n][i][0]*exp_h0;
77 for (
size_t j = 1; j <
_h[
n][i].size(); ++j)
78 gnj *= gsl_sf_pow_int(hn[j],
_h[
n][i][j]);
89std::vector< std::vector< std::vector<int> > >
Enge::_q;
90std::vector< std::vector< std::vector<int> > >
Enge::_h;
93 _q.push_back(std::vector< std::vector<int> >(1, std::vector<int>(3)) );
99 for (
size_t i =
_q.size(); i <
n+1; ++i) {
100 _q.push_back(std::vector< std::vector<int> >() );
101 for (
size_t j = 0; j <
_q[i-1].size(); ++j) {
102 size_t k_max =
_q[i-1][j].size();
103 std::vector<int> new_vec(
_q[i-1][j]);
105 new_vec[0] *= new_vec[1];
108 _q[i].push_back(new_vec);
109 for (
size_t k = 2; k < k_max; ++k) {
111 if (
_q[i-1][j][k] > 0) {
112 std::vector<int> new_vec(
_q[i-1][j]);
113 if ( k == k_max-1 ) new_vec.push_back(0);
114 new_vec[0] *= new_vec[k];
117 _q[i].push_back(new_vec);
123 if (
_h.size() == 0) {
125 _h.push_back(std::vector< std::vector<int> >());
127 _h.push_back(std::vector< std::vector<int> >());
128 _h[1].push_back(std::vector<int>(2, 1));
130 for (
size_t i =
_h.size(); i <
n+1; ++i) {
131 _h.push_back(std::vector< std::vector<int> >());
132 for (
size_t j = 0; j <
_h[i-1].size(); ++j) {
135 std::vector<int> new_vec(
_h[i-1][j]);
137 _h[i].push_back(new_vec);
138 for (
size_t k = 1; k <
_h[i-1][j].size(); ++k) {
139 if (
_h[i-1][j][k] > 0) {
140 std::vector<int> new_vec(
_h[i-1][j]);
141 if (k ==
_h[i-1][j].size()-1) new_vec.push_back(0);
142 new_vec[0] *= new_vec[k];
145 _h[i].push_back(new_vec);
154 : _a(
a), _lambda(lambda), _x0(x0) {
168 out <<
"Enge function l=" <<
_lambda <<
" x0=" <<
_x0 <<
" c=";
Tps< T > exp(const Tps< T > &x)
Exponential.
std::vector< std::vector< int > > CompactVector(std::vector< std::vector< int > > vec)
constexpr double e
The value of.
void rescale(double scaleFactor)
static std::vector< std::vector< std::vector< int > > > _h
static std::vector< std::vector< std::vector< int > > > _q
double gN(double x, int n) const
static void setEngeDiffIndices(size_t n)
static std::vector< std::vector< int > > getQIndex(int n)
double hN(double x, int n) const
double getEnge(double x, int n) const
std::ostream & print(std::ostream &out) const