Blame specfunc/cheb_eval_mode.c
|
Packit |
67cb25 |
static inline int
|
|
Packit |
67cb25 |
cheb_eval_mode_e(const cheb_series * cs,
|
|
Packit |
67cb25 |
const double x,
|
|
Packit |
67cb25 |
gsl_mode_t mode,
|
|
Packit |
67cb25 |
gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
int j;
|
|
Packit |
67cb25 |
double d = 0.0;
|
|
Packit |
67cb25 |
double dd = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double y = (2.*x - cs->a - cs->b) / (cs->b - cs->a);
|
|
Packit |
67cb25 |
double y2 = 2.0 * y;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int eval_order;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(GSL_MODE_PREC(mode) == GSL_PREC_DOUBLE)
|
|
Packit |
67cb25 |
eval_order = cs->order;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
eval_order = cs->order_sp;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(j = eval_order; j>=1; j--) {
|
|
Packit |
67cb25 |
double temp = d;
|
|
Packit |
67cb25 |
d = y2*d - dd + cs->c[j];
|
|
Packit |
67cb25 |
dd = temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = y*d - dd + 0.5 * cs->c[0];
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(cs->c[eval_order]);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|