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

   1 /* specfunc/bessel_K1.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 besk1(), besk1e() */
  36
  37 /* chebyshev expansions
  38
  39  series for bk1        on the interval  0.          to  4.00000d+00
  40                                         with weighted error   7.02e-18
  41                                          log weighted error  17.15
  42                                significant figures required  16.73
  43                                     decimal places required  17.67
  44
  45  series for ak1        on the interval  1.25000d-01 to  5.00000d-01
  46                                         with weighted error   6.06e-17
  47                                          log weighted error  16.22
  48                                significant figures required  15.41
  49                                     decimal places required  16.83
  50
  51  series for ak12       on the interval  0.          to  1.25000d-01
  52                                         with weighted error   2.58e-17
  53                                          log weighted error  16.59
  54                                significant figures required  15.22
  55                                     decimal places required  17.16
  56 */
  57
  58 static double bk1_data[11] = {
  59    0.0253002273389477705,
  60   -0.3531559607765448760,
  61   -0.1226111808226571480,
  62   -0.0069757238596398643,
  63   -0.0001730288957513052,
  64   -0.0000024334061415659,
  65   -0.0000000221338763073,
  66   -0.0000000001411488392,
  67   -0.0000000000006666901,
  68   -0.0000000000000024274,
  69   -0.0000000000000000070
  70 };
  71
  72 static cheb_series bk1_cs = {
  73   bk1_data,
  74   10,
  75   -1, 1,
  76   8
  77 };
  78
  79 static double ak1_data[17] = {
  80    0.27443134069738830,
  81    0.07571989953199368,
  82   -0.00144105155647540,
  83    0.00006650116955125,
  84   -0.00000436998470952,
  85    0.00000035402774997,
  86   -0.00000003311163779,
  87    0.00000000344597758,
  88   -0.00000000038989323,
  89    0.00000000004720819,
  90   -0.00000000000604783,
  91    0.00000000000081284,
  92   -0.00000000000011386,
  93    0.00000000000001654,
  94   -0.00000000000000248,
  95    0.00000000000000038,
  96   -0.00000000000000006
  97 };
  98 static cheb_series ak1_cs = {
  99   ak1_data,
 100   16,
 101   -1, 1,
 102   9
 103 };
 104
 105 static double ak12_data[14] = {
 106    0.06379308343739001,
 107    0.02832887813049721,
 108   -0.00024753706739052,
 109    0.00000577197245160,
 110   -0.00000020689392195,
 111    0.00000000973998344,
 112   -0.00000000055853361,
 113    0.00000000003732996,
 114   -0.00000000000282505,
 115    0.00000000000023720,
 116   -0.00000000000002176,
 117    0.00000000000000215,
 118   -0.00000000000000022,
 119    0.00000000000000002
 120 };
 121 static cheb_series ak12_cs = {
 122   ak12_data,
 123   13,
 124   -1, 1,
 125   7
 126 };
 127
 128
 129 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 130
 131 int gsl_sf_bessel_K1_scaled_e(const double x, gsl_sf_result * result)
 132 {
 133   /* CHECK_POINTER(result) */
 134
 135   if(x <= 0.0) {
 136     DOMAIN_ERROR(result);
 137   }
 138   else if(x < 2.0*GSL_DBL_MIN) {
 139     OVERFLOW_ERROR(result);
 140   }
 141   else if(x <= 2.0) {
 142     const double lx = log(x);
 143     const double ex = exp(x);
 144     int stat_I1;
 145     gsl_sf_result I1;
 146     gsl_sf_result c;
 147     cheb_eval_e(&bk1_cs, 0.5*x*x-1.0, &c);
 148     stat_I1 = gsl_sf_bessel_I1_e(x, &I1);
 149     result->val  = ex * ((lx-M_LN2)*I1.val + (0.75 + c.val)/x);
 150     result->err  = ex * (c.err/x + fabs(lx)*I1.err);
 151     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 152     return stat_I1;
 153   }
 154   else if(x <= 8.0) {
 155     const double sx = sqrt(x);
 156     gsl_sf_result c;
 157     cheb_eval_e(&ak1_cs, (16.0/x-5.0)/3.0, &c);
 158     result->val  = (1.25 + c.val) / sx;
 159     result->err  = c.err / sx;
 160     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 161     return GSL_SUCCESS;
 162   }
 163   else {
 164     const double sx = sqrt(x);
 165     gsl_sf_result c;
 166     cheb_eval_e(&ak12_cs, 16.0/x-1.0, &c);
 167     result->val  = (1.25 + c.val) / sx;
 168     result->err  = c.err / sx;
 169     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 170     return GSL_SUCCESS;
 171   }
 172 }
 173
 174
 175 int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result)
 176 {
 177   /* CHECK_POINTER(result) */
 178
 179   if(x <= 0.0) {
 180     DOMAIN_ERROR(result);
 181   }
 182   else if(x < 2.0*GSL_DBL_MIN) {
 183     OVERFLOW_ERROR(result);
 184   }
 185   else if(x <= 2.0) {
 186     const double lx = log(x);
 187     int stat_I1;
 188     gsl_sf_result I1;
 189     gsl_sf_result c;
 190     cheb_eval_e(&bk1_cs, 0.5*x*x-1.0, &c);
 191     stat_I1 = gsl_sf_bessel_I1_e(x, &I1);
 192     result->val  = (lx-M_LN2)*I1.val + (0.75 + c.val)/x;
 193     result->err  = c.err/x + fabs(lx)*I1.err;
 194     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 195     return stat_I1;
 196   }
 197   else {
 198     gsl_sf_result K1_scaled;
 199     int stat_K1 = gsl_sf_bessel_K1_scaled_e(x, &K1_scaled);
 200     int stat_e  = gsl_sf_exp_mult_err_e(-x, 0.0,
 201                                            K1_scaled.val, K1_scaled.err,
 202                                            result);
 203     result->err = fabs(result->val) * (GSL_DBL_EPSILON*fabs(x) + K1_scaled.err/K1_scaled.val);
 204     return GSL_ERROR_SELECT_2(stat_e, stat_K1);
 205   }
 206 }
 207
 208 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 209
 210 #include "eval.h"
 211
 212 double gsl_sf_bessel_K1_scaled(const double x)
 213 {
 214   EVAL_RESULT(gsl_sf_bessel_K1_scaled_e(x, &result));
 215 }
 216
 217 double gsl_sf_bessel_K1(const double x)
 218 {
 219   EVAL_RESULT(gsl_sf_bessel_K1_e(x, &result));
 220 }