|
Packit |
67cb25 |
/* specfunc/debye.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
/* augmented to n=5 and 6 2005-11-08 by R. J. Mathar, http://www.strw.leidenuniv.nl/~mathar */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_debye.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
#include "check.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "chebyshev.h"
|
|
Packit |
67cb25 |
#include "cheb_eval.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double adeb1_data[17] = {
|
|
Packit |
67cb25 |
2.4006597190381410194,
|
|
Packit |
67cb25 |
0.1937213042189360089,
|
|
Packit |
67cb25 |
-0.62329124554895770e-02,
|
|
Packit |
67cb25 |
0.3511174770206480e-03,
|
|
Packit |
67cb25 |
-0.228222466701231e-04,
|
|
Packit |
67cb25 |
0.15805467875030e-05,
|
|
Packit |
67cb25 |
-0.1135378197072e-06,
|
|
Packit |
67cb25 |
0.83583361188e-08,
|
|
Packit |
67cb25 |
-0.6264424787e-09,
|
|
Packit |
67cb25 |
0.476033489e-10,
|
|
Packit |
67cb25 |
-0.36574154e-11,
|
|
Packit |
67cb25 |
0.2835431e-12,
|
|
Packit |
67cb25 |
-0.221473e-13,
|
|
Packit |
67cb25 |
0.17409e-14,
|
|
Packit |
67cb25 |
-0.1376e-15,
|
|
Packit |
67cb25 |
0.109e-16,
|
|
Packit |
67cb25 |
-0.9e-18
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series adeb1_cs = {
|
|
Packit |
67cb25 |
adeb1_data,
|
|
Packit |
67cb25 |
16,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
9
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double adeb2_data[18] = {
|
|
Packit |
67cb25 |
2.5943810232570770282,
|
|
Packit |
67cb25 |
0.2863357204530719834,
|
|
Packit |
67cb25 |
-0.102062656158046713e-01,
|
|
Packit |
67cb25 |
0.6049109775346844e-03,
|
|
Packit |
67cb25 |
-0.405257658950210e-04,
|
|
Packit |
67cb25 |
0.28633826328811e-05,
|
|
Packit |
67cb25 |
-0.2086394303065e-06,
|
|
Packit |
67cb25 |
0.155237875826e-07,
|
|
Packit |
67cb25 |
-0.11731280087e-08,
|
|
Packit |
67cb25 |
0.897358589e-10,
|
|
Packit |
67cb25 |
-0.69317614e-11,
|
|
Packit |
67cb25 |
0.5398057e-12,
|
|
Packit |
67cb25 |
-0.423241e-13,
|
|
Packit |
67cb25 |
0.33378e-14,
|
|
Packit |
67cb25 |
-0.2645e-15,
|
|
Packit |
67cb25 |
0.211e-16,
|
|
Packit |
67cb25 |
-0.17e-17,
|
|
Packit |
67cb25 |
0.1e-18
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series adeb2_cs = {
|
|
Packit |
67cb25 |
adeb2_data,
|
|
Packit |
67cb25 |
17,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double adeb3_data[17] = {
|
|
Packit |
67cb25 |
2.707737068327440945,
|
|
Packit |
67cb25 |
0.340068135211091751,
|
|
Packit |
67cb25 |
-0.12945150184440869e-01,
|
|
Packit |
67cb25 |
0.7963755380173816e-03,
|
|
Packit |
67cb25 |
-0.546360009590824e-04,
|
|
Packit |
67cb25 |
0.39243019598805e-05,
|
|
Packit |
67cb25 |
-0.2894032823539e-06,
|
|
Packit |
67cb25 |
0.217317613962e-07,
|
|
Packit |
67cb25 |
-0.16542099950e-08,
|
|
Packit |
67cb25 |
0.1272796189e-09,
|
|
Packit |
67cb25 |
-0.987963460e-11,
|
|
Packit |
67cb25 |
0.7725074e-12,
|
|
Packit |
67cb25 |
-0.607797e-13,
|
|
Packit |
67cb25 |
0.48076e-14,
|
|
Packit |
67cb25 |
-0.3820e-15,
|
|
Packit |
67cb25 |
0.305e-16,
|
|
Packit |
67cb25 |
-0.24e-17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series adeb3_cs = {
|
|
Packit |
67cb25 |
adeb3_data,
|
|
Packit |
67cb25 |
16,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double adeb4_data[17] = {
|
|
Packit |
67cb25 |
2.781869415020523460,
|
|
Packit |
67cb25 |
0.374976783526892863,
|
|
Packit |
67cb25 |
-0.14940907399031583e-01,
|
|
Packit |
67cb25 |
0.945679811437042e-03,
|
|
Packit |
67cb25 |
-0.66132916138933e-04,
|
|
Packit |
67cb25 |
0.4815632982144e-05,
|
|
Packit |
67cb25 |
-0.3588083958759e-06,
|
|
Packit |
67cb25 |
0.271601187416e-07,
|
|
Packit |
67cb25 |
-0.20807099122e-08,
|
|
Packit |
67cb25 |
0.1609383869e-09,
|
|
Packit |
67cb25 |
-0.125470979e-10,
|
|
Packit |
67cb25 |
0.9847265e-12,
|
|
Packit |
67cb25 |
-0.777237e-13,
|
|
Packit |
67cb25 |
0.61648e-14,
|
|
Packit |
67cb25 |
-0.4911e-15,
|
|
Packit |
67cb25 |
0.393e-16,
|
|
Packit |
67cb25 |
-0.32e-17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series adeb4_cs = {
|
|
Packit |
67cb25 |
adeb4_data,
|
|
Packit |
67cb25 |
16,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double adeb5_data[17] = {
|
|
Packit |
67cb25 |
2.8340269546834530149,
|
|
Packit |
67cb25 |
0.3994098857106266445,
|
|
Packit |
67cb25 |
-0.164566764773099646e-1,
|
|
Packit |
67cb25 |
0.10652138340664541e-2,
|
|
Packit |
67cb25 |
-0.756730374875418e-4,
|
|
Packit |
67cb25 |
0.55745985240273e-5,
|
|
Packit |
67cb25 |
-0.4190692330918e-6,
|
|
Packit |
67cb25 |
0.319456143678e-7,
|
|
Packit |
67cb25 |
-0.24613318171e-8,
|
|
Packit |
67cb25 |
0.1912801633e-9,
|
|
Packit |
67cb25 |
-0.149720049e-10,
|
|
Packit |
67cb25 |
0.11790312e-11,
|
|
Packit |
67cb25 |
-0.933329e-13,
|
|
Packit |
67cb25 |
0.74218e-14,
|
|
Packit |
67cb25 |
-0.5925e-15,
|
|
Packit |
67cb25 |
0.475e-16,
|
|
Packit |
67cb25 |
-0.39e-17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series adeb5_cs = {
|
|
Packit |
67cb25 |
adeb5_data,
|
|
Packit |
67cb25 |
16,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double adeb6_data[17] = {
|
|
Packit |
67cb25 |
2.8726727134130122113,
|
|
Packit |
67cb25 |
0.4174375352339027746,
|
|
Packit |
67cb25 |
-0.176453849354067873e-1,
|
|
Packit |
67cb25 |
0.11629852733494556e-2,
|
|
Packit |
67cb25 |
-0.837118027357117e-4,
|
|
Packit |
67cb25 |
0.62283611596189e-5,
|
|
Packit |
67cb25 |
-0.4718644465636e-6,
|
|
Packit |
67cb25 |
0.361950397806e-7,
|
|
Packit |
67cb25 |
-0.28030368010e-8,
|
|
Packit |
67cb25 |
0.2187681983e-9,
|
|
Packit |
67cb25 |
-0.171857387e-10,
|
|
Packit |
67cb25 |
0.13575809e-11,
|
|
Packit |
67cb25 |
-0.1077580e-12,
|
|
Packit |
67cb25 |
0.85893e-14,
|
|
Packit |
67cb25 |
-0.6872e-15,
|
|
Packit |
67cb25 |
0.552e-16,
|
|
Packit |
67cb25 |
-0.44e-17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
static cheb_series adeb6_cs = {
|
|
Packit |
67cb25 |
adeb6_data,
|
|
Packit |
67cb25 |
16,
|
|
Packit |
67cb25 |
-1.0, 1.0,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_debye_1_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double val_infinity = 1.64493406684822644;
|
|
Packit |
67cb25 |
const double xcut = -GSL_LOG_DBL_MIN;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = 1.0 - 0.25*x + x*x/36.0;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&adeb1_cs, t, &c);
|
|
Packit |
67cb25 |
result->val = c.val - 0.25 * x;
|
|
Packit |
67cb25 |
result->err = c.err + 0.25 * x * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
|
|
Packit |
67cb25 |
const int nexp = floor(xcut/x);
|
|
Packit |
67cb25 |
const double ex = exp(-x);
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
double xk = nexp * x;
|
|
Packit |
67cb25 |
double rk = nexp;
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
for(i=nexp; i>=1; i--) {
|
|
Packit |
67cb25 |
sum *= ex;
|
|
Packit |
67cb25 |
sum += (1.0 + 1.0/xk)/rk;
|
|
Packit |
67cb25 |
rk -= 1.0;
|
|
Packit |
67cb25 |
xk -= x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = val_infinity/x - sum*ex;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < xcut) {
|
|
Packit |
67cb25 |
result->val = (val_infinity - exp(-x)*(x+1.0)) / x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = val_infinity/x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_debye_2_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double val_infinity = 4.80822761263837714;
|
|
Packit |
67cb25 |
const double xcut = -GSL_LOG_DBL_MIN;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = 1.0 - x/3.0 + x*x/24.0;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&adeb2_cs, t, &c);
|
|
Packit |
67cb25 |
result->val = c.val - x/3.0;
|
|
Packit |
67cb25 |
result->err = c.err + GSL_DBL_EPSILON * x/3.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
|
|
Packit |
67cb25 |
const int nexp = floor(xcut/x);
|
|
Packit |
67cb25 |
const double ex = exp(-x);
|
|
Packit |
67cb25 |
double xk = nexp * x;
|
|
Packit |
67cb25 |
double rk = nexp;
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
for(i=nexp; i>=1; i--) {
|
|
Packit |
67cb25 |
sum *= ex;
|
|
Packit |
67cb25 |
sum += (1.0 + 2.0/xk + 2.0/(xk*xk)) / rk;
|
|
Packit |
67cb25 |
rk -= 1.0;
|
|
Packit |
67cb25 |
xk -= x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = val_infinity/(x*x) - 2.0 * sum * ex;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < xcut) {
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
const double sum = 2.0 + 2.0*x + x2;
|
|
Packit |
67cb25 |
result->val = (val_infinity - 2.0 * sum * exp(-x)) / x2;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = (val_infinity/x)/x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
CHECK_UNDERFLOW(result);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_debye_3_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double val_infinity = 19.4818182068004875;
|
|
Packit |
67cb25 |
const double xcut = -GSL_LOG_DBL_MIN;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = 1.0 - 3.0*x/8.0 + x*x/20.0;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&adeb3_cs, t, &c);
|
|
Packit |
67cb25 |
result->val = c.val - 0.375*x;
|
|
Packit |
67cb25 |
result->err = c.err + GSL_DBL_EPSILON * 0.375*x;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
|
|
Packit |
67cb25 |
const int nexp = floor(xcut/x);
|
|
Packit |
67cb25 |
const double ex = exp(-x);
|
|
Packit |
67cb25 |
double xk = nexp * x;
|
|
Packit |
67cb25 |
double rk = nexp;
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
for(i=nexp; i>=1; i--) {
|
|
Packit |
67cb25 |
double xk_inv = 1.0/xk;
|
|
Packit |
67cb25 |
sum *= ex;
|
|
Packit |
67cb25 |
sum += (((6.0*xk_inv + 6.0)*xk_inv + 3.0)*xk_inv + 1.0) / rk;
|
|
Packit |
67cb25 |
rk -= 1.0;
|
|
Packit |
67cb25 |
xk -= x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = val_infinity/(x*x*x) - 3.0 * sum * ex;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < xcut) {
|
|
Packit |
67cb25 |
const double x3 = x*x*x;
|
|
Packit |
67cb25 |
const double sum = 6.0 + 6.0*x + 3.0*x*x + x3;
|
|
Packit |
67cb25 |
result->val = (val_infinity - 3.0 * sum * exp(-x)) / x3;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = ((val_infinity/x)/x)/x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
CHECK_UNDERFLOW(result);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_debye_4_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double val_infinity = 99.5450644937635129;
|
|
Packit |
67cb25 |
const double xcut = -GSL_LOG_DBL_MIN;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = 1.0 - 2.0*x/5.0 + x*x/18.0;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&adeb4_cs, t, &c);
|
|
Packit |
67cb25 |
result->val = c.val - 2.0*x/5.0;
|
|
Packit |
67cb25 |
result->err = c.err + GSL_DBL_EPSILON * 2.0*x/5.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
|
|
Packit |
67cb25 |
const int nexp = floor(xcut/x);
|
|
Packit |
67cb25 |
const double ex = exp(-x);
|
|
Packit |
67cb25 |
double xk = nexp * x;
|
|
Packit |
67cb25 |
double rk = nexp;
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
for(i=nexp; i>=1; i--) {
|
|
Packit |
67cb25 |
double xk_inv = 1.0/xk;
|
|
Packit |
67cb25 |
sum *= ex;
|
|
Packit |
67cb25 |
sum += ((((24.0*xk_inv + 24.0)*xk_inv + 12.0)*xk_inv + 4.0)*xk_inv + 1.0) / rk;
|
|
Packit |
67cb25 |
rk -= 1.0;
|
|
Packit |
67cb25 |
xk -= x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = val_infinity/(x*x*x*x) - 4.0 * sum * ex;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < xcut) {
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
const double x4 = x2*x2;
|
|
Packit |
67cb25 |
const double sum = 24.0 + 24.0*x + 12.0*x2 + 4.0*x2*x + x4;
|
|
Packit |
67cb25 |
result->val = (val_infinity - 4.0 * sum * exp(-x)) / x4;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = (((val_infinity/x)/x)/x)/x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
CHECK_UNDERFLOW(result);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_debye_5_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double val_infinity = 610.405837190669483828710757875 ;
|
|
Packit |
67cb25 |
const double xcut = -GSL_LOG_DBL_MIN;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = 1.0 - 5.0*x/12.0 + 5.0*x*x/84.0;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&adeb5_cs, t, &c);
|
|
Packit |
67cb25 |
result->val = c.val - 5.0*x/12.0;
|
|
Packit |
67cb25 |
result->err = c.err + GSL_DBL_EPSILON * 5.0*x/12.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
|
|
Packit |
67cb25 |
const int nexp = floor(xcut/x);
|
|
Packit |
67cb25 |
const double ex = exp(-x);
|
|
Packit |
67cb25 |
double xk = nexp * x;
|
|
Packit |
67cb25 |
double rk = nexp;
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
for(i=nexp; i>=1; i--) {
|
|
Packit |
67cb25 |
double xk_inv = 1.0/xk;
|
|
Packit |
67cb25 |
sum *= ex;
|
|
Packit |
67cb25 |
sum += (((((120.0*xk_inv + 120.0)*xk_inv + 60.0)*xk_inv + 20.0)*xk_inv + 5.0)*xk_inv+ 1.0) / rk;
|
|
Packit |
67cb25 |
rk -= 1.0;
|
|
Packit |
67cb25 |
xk -= x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = val_infinity/(x*x*x*x*x) - 5.0 * sum * ex;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < xcut) {
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
const double x4 = x2*x2;
|
|
Packit |
67cb25 |
const double x5 = x4*x;
|
|
Packit |
67cb25 |
const double sum = 120.0 + 120.0*x + 60.0*x2 + 20.0*x2*x + 5.0*x4 + x5;
|
|
Packit |
67cb25 |
result->val = (val_infinity - 5.0 * sum * exp(-x)) / x5;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = ((((val_infinity/x)/x)/x)/x)/x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
CHECK_UNDERFLOW(result);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_debye_6_e(const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double val_infinity = 4356.06887828990661194792541535 ;
|
|
Packit |
67cb25 |
const double xcut = -GSL_LOG_DBL_MIN;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
result->val = 1.0 - 3.0*x/7.0 + x*x/16.0;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 4.0) {
|
|
Packit |
67cb25 |
const double t = x*x/8.0 - 1.0;
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&adeb6_cs, t, &c);
|
|
Packit |
67cb25 |
result->val = c.val - 3.0*x/7.0;
|
|
Packit |
67cb25 |
result->err = c.err + GSL_DBL_EPSILON * 3.0*x/7.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {
|
|
Packit |
67cb25 |
const int nexp = floor(xcut/x);
|
|
Packit |
67cb25 |
const double ex = exp(-x);
|
|
Packit |
67cb25 |
double xk = nexp * x;
|
|
Packit |
67cb25 |
double rk = nexp;
|
|
Packit |
67cb25 |
double sum = 0.0;
|
|
Packit |
67cb25 |
int i;
|
|
Packit |
67cb25 |
for(i=nexp; i>=1; i--) {
|
|
Packit |
67cb25 |
double xk_inv = 1.0/xk;
|
|
Packit |
67cb25 |
sum *= ex;
|
|
Packit |
67cb25 |
sum += ((((((720.0*xk_inv + 720.0)*xk_inv + 360.0)*xk_inv + 120.0)*xk_inv + 30.0)*xk_inv+ 6.0)*xk_inv+ 1.0) / rk;
|
|
Packit |
67cb25 |
rk -= 1.0;
|
|
Packit |
67cb25 |
xk -= x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = val_infinity/(x*x*x*x*x*x) - 6.0 * sum * ex;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < xcut) {
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
const double x4 = x2*x2;
|
|
Packit |
67cb25 |
const double x6 = x4*x2;
|
|
Packit |
67cb25 |
const double sum = 720.0 + 720.0*x + 360.0*x2 + 120.0*x2*x + 30.0*x4 + 6.0*x4*x +x6 ;
|
|
Packit |
67cb25 |
result->val = (val_infinity - 6.0 * sum * exp(-x)) / x6;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = (((((val_infinity/x)/x)/x)/x)/x)/x ;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * result->val;
|
|
Packit |
67cb25 |
CHECK_UNDERFLOW(result);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_debye_1(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_debye_1_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_debye_2(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_debye_2_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_debye_3(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_debye_3_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_debye_4(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_debye_4_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_debye_5(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_debye_5_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_debye_6(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_debye_6_e(x, &result));
|
|
Packit |
67cb25 |
}
|