|
Packit |
67cb25 |
/* specfunc/bessel_K0.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
* Copyright (C) 2016 Pavel Holoborodko, Patrick Alken
|
|
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 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_poly.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_exp.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_bessel.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "chebyshev.h"
|
|
Packit |
67cb25 |
#include "cheb_eval.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
Minimax rational approximation for [0,1), peak relative error = 2.04*GSL_DBL_EPSILON.
|
|
Packit |
67cb25 |
Source: http://www.advanpix.com/?p=3812
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double k0_poly[8] = {
|
|
Packit |
67cb25 |
1.1593151565841244842077226e-01,
|
|
Packit |
67cb25 |
2.7898287891460317300886539e-01,
|
|
Packit |
67cb25 |
2.5248929932161220559969776e-02,
|
|
Packit |
67cb25 |
8.4603509072136578707676406e-04,
|
|
Packit |
67cb25 |
1.4914719243067801775856150e-05,
|
|
Packit |
67cb25 |
1.6271068931224552553548933e-07,
|
|
Packit |
67cb25 |
1.2082660336282566759313543e-09,
|
|
Packit |
67cb25 |
6.6117104672254184399933971e-12
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static double i0_poly[7] = {
|
|
Packit |
67cb25 |
1.0000000000000000044974165e+00,
|
|
Packit |
67cb25 |
2.4999999999999822316775454e-01,
|
|
Packit |
67cb25 |
2.7777777777892149148858521e-02,
|
|
Packit |
67cb25 |
1.7361111083544590676709592e-03,
|
|
Packit |
67cb25 |
6.9444476047072424198677755e-05,
|
|
Packit |
67cb25 |
1.9288265756466775034067979e-06,
|
|
Packit |
67cb25 |
3.9908220583262192851839992e-08
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
Chebyshev expansion for [1,8], peak relative error = 1.28*GSL_DBL_EPSILON.
|
|
Packit |
67cb25 |
Source: Pavel Holoborodko.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double ak0_data[24] = {
|
|
Packit |
67cb25 |
-3.28737867094650101e-02,
|
|
Packit |
67cb25 |
-4.49369057710236880e-02,
|
|
Packit |
67cb25 |
+2.98149992004308095e-03,
|
|
Packit |
67cb25 |
-3.03693649396187920e-04,
|
|
Packit |
67cb25 |
+3.91085569307646836e-05,
|
|
Packit |
67cb25 |
-5.86872422399215952e-06,
|
|
Packit |
67cb25 |
+9.82873709937322009e-07,
|
|
Packit |
67cb25 |
-1.78978645055651171e-07,
|
|
Packit |
67cb25 |
+3.48332306845240957e-08,
|
|
Packit |
67cb25 |
-7.15909210462546599e-09,
|
|
Packit |
67cb25 |
+1.54019930048919494e-09,
|
|
Packit |
67cb25 |
-3.44555485579194210e-10,
|
|
Packit |
67cb25 |
+7.97356101783753023e-11,
|
|
Packit |
67cb25 |
-1.90090968913069735e-11,
|
|
Packit |
67cb25 |
+4.65295609304114621e-12,
|
|
Packit |
67cb25 |
-1.16614287433470780e-12,
|
|
Packit |
67cb25 |
+2.98554375218596891e-13,
|
|
Packit |
67cb25 |
-7.79276979512292169e-14,
|
|
Packit |
67cb25 |
+2.07027467168948402e-14,
|
|
Packit |
67cb25 |
-5.58987860393825313e-15,
|
|
Packit |
67cb25 |
+1.53202965950646914e-15,
|
|
Packit |
67cb25 |
-4.25737536712188186e-16,
|
|
Packit |
67cb25 |
+1.19840238501357389e-16,
|
|
Packit |
67cb25 |
-3.41407346762502397e-17
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static cheb_series ak0_cs = {
|
|
Packit |
67cb25 |
ak0_data,
|
|
Packit |
67cb25 |
23,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
10
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*
|
|
Packit |
67cb25 |
Chebyshev expansion for [8,inf), peak relative error = 1.25*GSL_DBL_EPSILON.
|
|
Packit |
67cb25 |
Source: SLATEC/dbsk0e.f
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static double ak02_data[14] = {
|
|
Packit |
67cb25 |
-.1201869826307592240E-1,
|
|
Packit |
67cb25 |
-.9174852691025695311E-2,
|
|
Packit |
67cb25 |
+.1444550931775005821E-3,
|
|
Packit |
67cb25 |
-.4013614175435709729E-5,
|
|
Packit |
67cb25 |
+.1567831810852310673E-6,
|
|
Packit |
67cb25 |
-.7770110438521737710E-8,
|
|
Packit |
67cb25 |
+.4611182576179717883E-9,
|
|
Packit |
67cb25 |
-.3158592997860565771E-10,
|
|
Packit |
67cb25 |
+.2435018039365041128E-11,
|
|
Packit |
67cb25 |
-.2074331387398347898E-12,
|
|
Packit |
67cb25 |
+.1925787280589917085E-13,
|
|
Packit |
67cb25 |
-.1927554805838956104E-14,
|
|
Packit |
67cb25 |
+.2062198029197818278E-15,
|
|
Packit |
67cb25 |
-.2341685117579242403E-16
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static cheb_series ak02_cs = {
|
|
Packit |
67cb25 |
ak02_data,
|
|
Packit |
67cb25 |
13,
|
|
Packit |
67cb25 |
-1, 1,
|
|
Packit |
67cb25 |
8
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result)
|
|
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 < 1.0) {
|
|
Packit |
67cb25 |
const double lx = log(x);
|
|
Packit |
67cb25 |
const double ex = exp(x);
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
result->val = ex * (gsl_poly_eval(k0_poly,8,x2)-lx*(1.0+0.25*x2*gsl_poly_eval(i0_poly,7,0.25*x2)));
|
|
Packit |
67cb25 |
result->err = ex * (1.6+fabs(lx)*0.6) * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 8.0) {
|
|
Packit |
67cb25 |
const double sx = sqrt(x);
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&ak0_cs, (16.0/x-9.0)/7.0, &c);
|
|
Packit |
67cb25 |
result->val = (1.203125 + c.val) / sx; /* 1.203125 = 77/64 */
|
|
Packit |
67cb25 |
result->err = c.err / sx;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double sx = sqrt(x);
|
|
Packit |
67cb25 |
gsl_sf_result c;
|
|
Packit |
67cb25 |
cheb_eval_e(&ak02_cs, 16.0/x-1.0, &c);
|
|
Packit |
67cb25 |
result->val = (1.25 + c.val) / sx;
|
|
Packit |
67cb25 |
result->err = (c.err + GSL_DBL_EPSILON) / sx;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result)
|
|
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 < 1.0) {
|
|
Packit |
67cb25 |
const double lx = log(x);
|
|
Packit |
67cb25 |
const double x2 = x*x;
|
|
Packit |
67cb25 |
result->val = gsl_poly_eval(k0_poly,8,x2)-lx*(1.0+0.25*x2*gsl_poly_eval(i0_poly,7,0.25*x2));
|
|
Packit |
67cb25 |
result->err = (1.6+fabs(lx)*0.6) * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result K0_scaled;
|
|
Packit |
67cb25 |
int stat_K0 = gsl_sf_bessel_K0_scaled_e(x, &K0_scaled);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(-x, GSL_DBL_EPSILON*fabs(x),
|
|
Packit |
67cb25 |
K0_scaled.val, K0_scaled.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K0);
|
|
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_bessel_K0_scaled(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_bessel_K0_scaled_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_bessel_K0(const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_bessel_K0_e(x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|