/* specfunc/bessel_I0.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

/* Author:  G. Jungman */

#include <config.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "chebyshev.h"

#include "cheb_eval.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on SLATEC besi0 */

/* chebyshev expansions

 series for bi0        on the interval  0.          to  9.00000d+00

                                        with weighted error   2.46e-18

                                         log weighted error  17.61

                               significant figures required  17.90

                                    decimal places required  18.15

 series for ai0        on the interval  1.25000d-01 to  3.33333d-01

                                        with weighted error   7.87e-17

                                         log weighted error  16.10

                               significant figures required  14.69

                                    decimal places required  16.76

 series for ai02       on the interval  0.          to  1.25000d-01

                                        with weighted error   3.79e-17

                                         log weighted error  16.42

                               significant figures required  14.86

                                    decimal places required  17.09

*/

static double bi0_data[12] = {

  -.07660547252839144951,

  1.92733795399380827000,

   .22826445869203013390, 

   .01304891466707290428,

   .00043442709008164874,

   .00000942265768600193,

   .00000014340062895106,

   .00000000161384906966,

   .00000000001396650044,

   .00000000000009579451,

   .00000000000000053339,

   .00000000000000000245

};

static cheb_series bi0_cs = {

  bi0_data,

  11,

  -1, 1,

  11

};

static double ai0_data[21] = {

   .07575994494023796, 

   .00759138081082334,

   .00041531313389237,

   .00001070076463439,

  -.00000790117997921,

  -.00000078261435014,

   .00000027838499429,

   .00000000825247260,

  -.00000001204463945,

   .00000000155964859,

   .00000000022925563,

  -.00000000011916228,

   .00000000001757854,

   .00000000000112822,

  -.00000000000114684,

   .00000000000027155,

  -.00000000000002415,

  -.00000000000000608,

   .00000000000000314,

  -.00000000000000071,

   .00000000000000007

};

static cheb_series ai0_cs = {

  ai0_data,

  20,

  -1, 1,

  13

};

static double ai02_data[22] = {

   .05449041101410882,

   .00336911647825569,

   .00006889758346918,

   .00000289137052082,

   .00000020489185893,

   .00000002266668991,

   .00000000339623203,

   .00000000049406022,

   .00000000001188914,

  -.00000000003149915,

  -.00000000001321580,

  -.00000000000179419,

   .00000000000071801,

   .00000000000038529,

   .00000000000001539,

  -.00000000000004151,

  -.00000000000000954,

   .00000000000000382,

   .00000000000000176,

  -.00000000000000034,

  -.00000000000000027,

   .00000000000000003

};

static cheb_series ai02_cs = {

  ai02_data,

  21,

  -1, 1,

  11

};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_I0_scaled_e(const double x, gsl_sf_result * result)

{

  double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {

    result->val = 1.0 - y;

    result->err = 0.5*y*y;

    return GSL_SUCCESS;

  }

  else if(y <= 3.0) {

    const double ey = exp(-y);

    gsl_sf_result c;

    cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);

    result->val = ey * (2.75 + c.val);

    result->err = GSL_DBL_EPSILON * fabs(result->val) + ey * c.err;

    return GSL_SUCCESS;

  }

  else if(y <= 8.0) {

    const double sy = sqrt(y);

    gsl_sf_result c;

    cheb_eval_e(&ai0_cs, (48.0/y-11.0)/5.0, &c);

    result->val  = (0.375 + c.val) / sy;

    result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;

    result->err += c.err / sy;

    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_SUCCESS;

  }

  else {

    const double sy = sqrt(y);

    gsl_sf_result c;

    cheb_eval_e(&ai02_cs, 16.0/y-1.0, &c);

    result->val = (0.375 + c.val) / sy;

    result->err  = 2.0 * GSL_DBL_EPSILON * (0.375 + fabs(c.val)) / sy;

    result->err += c.err / sy;

    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_SUCCESS;

  }

}

int gsl_sf_bessel_I0_e(const double x, gsl_sf_result * result)

{

  double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y < 2.0 * GSL_SQRT_DBL_EPSILON) {

    result->val = 1.0;

    result->err = 0.5*y*y;

    return GSL_SUCCESS;

  }

  else if(y <= 3.0) {

    gsl_sf_result c;

    cheb_eval_e(&bi0_cs, y*y/4.5-1.0, &c);

    result->val  = 2.75 + c.val;

    result->err  = GSL_DBL_EPSILON * (2.75 + fabs(c.val));

    result->err += c.err;

    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_SUCCESS;

  }

  else if(y < GSL_LOG_DBL_MAX - 1.0) {

    const double ey = exp(y);

    gsl_sf_result b_scaled;

    gsl_sf_bessel_I0_scaled_e(x, &b_scaled);

    result->val  = ey * b_scaled.val;

    result->err  = ey * b_scaled.err + y*GSL_DBL_EPSILON*fabs(result->val);

    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

    return GSL_SUCCESS;

  }

  else {

    OVERFLOW_ERROR(result);

  }

}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_I0_scaled(const double x)

{

  EVAL_RESULT(gsl_sf_bessel_I0_scaled_e(x, &result); ) 

}

double gsl_sf_bessel_I0(const double x)

{

  EVAL_RESULT(gsl_sf_bessel_I0_e(x, &result); ) 

}