htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/bessel_K0.c

   1 /* specfunc/bessel_K0.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_exp.h>
  26 #include <gsl/gsl_sf_bessel.h>
  27
  28 #include "error.h"
  29
  30 #include "chebyshev.h"
  31 #include "cheb_eval.c"
  32
  33 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  34
  35 /* based on SLATEC bk0(), bk0e() */
  36
  37 /* chebyshev expansions
  38
  39  series for bk0        on the interval  0.          to  4.00000d+00
  40                                         with weighted error   3.57e-19
  41                                          log weighted error  18.45
  42                                significant figures required  17.99
  43                                     decimal places required  18.97
  44
  45  series for ak0        on the interval  1.25000d-01 to  5.00000d-01
  46                                         with weighted error   5.34e-17
  47                                          log weighted error  16.27
  48                                significant figures required  14.92
  49                                     decimal places required  16.89
  50
  51  series for ak02       on the interval  0.          to  1.25000d-01
  52                                         with weighted error   2.34e-17
  53                                          log weighted error  16.63
  54                                significant figures required  14.67
  55                                     decimal places required  17.20
  56 */
  57
  58 static double bk0_data[11] = {
  59   -0.03532739323390276872,
  60    0.3442898999246284869,
  61    0.03597993651536150163,
  62    0.00126461541144692592,
  63    0.00002286212103119451,
  64    0.00000025347910790261,
  65    0.00000000190451637722,
  66    0.00000000001034969525,
  67    0.00000000000004259816,
  68    0.00000000000000013744,
  69    0.00000000000000000035
  70 };
  71 static cheb_series bk0_cs = {
  72   bk0_data,
  73   10,
  74   -1, 1,
  75   10
  76 };
  77
  78 static double ak0_data[17] = {
  79   -0.07643947903327941,
  80   -0.02235652605699819,
  81    0.00077341811546938,
  82   -0.00004281006688886,
  83    0.00000308170017386,
  84   -0.00000026393672220,
  85    0.00000002563713036,
  86   -0.00000000274270554,
  87    0.00000000031694296,
  88   -0.00000000003902353,
  89    0.00000000000506804,
  90   -0.00000000000068895,
  91    0.00000000000009744,
  92   -0.00000000000001427,
  93    0.00000000000000215,
  94   -0.00000000000000033,
  95    0.00000000000000005
  96 };
  97 static cheb_series ak0_cs = {
  98   ak0_data,
  99   16,
 100   -1, 1,
 101   10
 102 };
 103
 104 static double ak02_data[14] = {
 105   -0.01201869826307592,
 106   -0.00917485269102569,
 107    0.00014445509317750,
 108   -0.00000401361417543,
 109    0.00000015678318108,
 110   -0.00000000777011043,
 111    0.00000000046111825,
 112   -0.00000000003158592,
 113    0.00000000000243501,
 114   -0.00000000000020743,
 115    0.00000000000001925,
 116   -0.00000000000000192,
 117    0.00000000000000020,
 118   -0.00000000000000002
 119 };
 120 static cheb_series ak02_cs = {
 121   ak02_data,
 122   13,
 123   -1, 1,
 124   8
 125 };
 126
 127
 128 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 129
 130 int gsl_sf_bessel_K0_scaled_e(const double x, gsl_sf_result * result)
 131 {
 132   /* CHECK_POINTER(result) */
 133
 134   if(x <= 0.0) {
 135     DOMAIN_ERROR(result);
 136   }
 137   else if(x <= 2.0) {
 138     const double lx = log(x);
 139     const double ex = exp(x);
 140     int stat_I0;
 141     gsl_sf_result I0;
 142     gsl_sf_result c;
 143     cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c);
 144     stat_I0 = gsl_sf_bessel_I0_e(x, &I0);
 145     result->val  = ex * ((-lx+M_LN2)*I0.val - 0.25 + c.val);
 146     result->err  = ex * ((M_LN2+fabs(lx))*I0.err + c.err);
 147     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 148     return stat_I0;
 149   }
 150   else if(x <= 8.0) {
 151     const double sx = sqrt(x);
 152     gsl_sf_result c;
 153     cheb_eval_e(&ak0_cs, (16.0/x-5.0)/3.0, &c);
 154     result->val  = (1.25 + c.val) / sx;
 155     result->err  = c.err / sx;
 156     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 157     return GSL_SUCCESS;
 158   }
 159   else {
 160     const double sx = sqrt(x);
 161     gsl_sf_result c;
 162     cheb_eval_e(&ak02_cs, 16.0/x-1.0, &c);
 163     result->val  = (1.25 + c.val) / sx;
 164     result->err  = (c.err + GSL_DBL_EPSILON) / sx;
 165     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 166     return GSL_SUCCESS;
 167   }
 168 }
 169
 170
 171 int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result)
 172 {
 173   /* CHECK_POINTER(result) */
 174
 175   if(x <= 0.0) {
 176     DOMAIN_ERROR(result);
 177   }
 178   else if(x <= 2.0) {
 179     const double lx = log(x);
 180     int stat_I0;
 181     gsl_sf_result I0;
 182     gsl_sf_result c;
 183     cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c);
 184     stat_I0 = gsl_sf_bessel_I0_e(x, &I0);
 185     result->val  = (-lx+M_LN2)*I0.val - 0.25 + c.val;
 186     result->err  = (fabs(lx) + M_LN2) * I0.err + c.err;
 187     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 188     return stat_I0;
 189   }
 190   else {
 191     gsl_sf_result K0_scaled;
 192     int stat_K0 = gsl_sf_bessel_K0_scaled_e(x, &K0_scaled);
 193     int stat_e  = gsl_sf_exp_mult_err_e(-x, GSL_DBL_EPSILON*fabs(x),
 194                                            K0_scaled.val, K0_scaled.err,
 195                                            result);
 196     return GSL_ERROR_SELECT_2(stat_e, stat_K0);
 197   }
 198 }
 199
 200
 201 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 202
 203 #include "eval.h"
 204
 205 double gsl_sf_bessel_K0_scaled(const double x)
 206 {
 207   EVAL_RESULT(gsl_sf_bessel_K0_scaled_e(x, &result));
 208 }
 209
 210 double gsl_sf_bessel_K0(const double x)
 211 {
 212   EVAL_RESULT(gsl_sf_bessel_K0_e(x, &result));
 213 }
 214