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

   1 /* specfunc/bessel_temme.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 /* Calculate series for Y_nu and K_nu for small x and nu.
  23  * This is applicable for x < 2 and |nu|<=1/2.
  24  * These functions assume x > 0.
  25  */
  26 #include <config.h>
  27 #include <gsl/gsl_math.h>
  28 #include <gsl/gsl_errno.h>
  29 #include <gsl/gsl_mode.h>
  30 #include "bessel_temme.h"
  31
  32 #include "chebyshev.h"
  33 #include "cheb_eval.c"
  34
  35 /* nu = (x+1)/4, -1<x<1, 1/(2nu)(1/Gamma[1-nu]-1/Gamma[1+nu]) */
  36 static double g1_dat[14] = {
  37   -1.14516408366268311786898152867,
  38    0.00636085311347084238122955495,
  39    0.00186245193007206848934643657,
  40    0.000152833085873453507081227824,
  41    0.000017017464011802038795324732,
  42   -6.4597502923347254354668326451e-07,
  43   -5.1819848432519380894104312968e-08,
  44    4.5189092894858183051123180797e-10,
  45    3.2433227371020873043666259180e-11,
  46    6.8309434024947522875432400828e-13,
  47    2.8353502755172101513119628130e-14,
  48   -7.9883905769323592875638087541e-16,
  49   -3.3726677300771949833341213457e-17,
  50   -3.6586334809210520744054437104e-20
  51 };
  52 static cheb_series g1_cs = {
  53   g1_dat,
  54   13,
  55   -1, 1,
  56   7
  57 };
  58
  59 /* nu = (x+1)/4, -1<x<1,  1/2 (1/Gamma[1-nu]+1/Gamma[1+nu]) */
  60 static double g2_dat[15] =
  61 {
  62   1.882645524949671835019616975350,
  63  -0.077490658396167518329547945212,
  64  -0.018256714847324929419579340950,
  65   0.0006338030209074895795923971731,
  66   0.0000762290543508729021194461175,
  67  -9.5501647561720443519853993526e-07,
  68  -8.8927268107886351912431512955e-08,
  69  -1.9521334772319613740511880132e-09,
  70  -9.4003052735885162111769579771e-11,
  71   4.6875133849532393179290879101e-12,
  72   2.2658535746925759582447545145e-13,
  73  -1.1725509698488015111878735251e-15,
  74  -7.0441338200245222530843155877e-17,
  75  -2.4377878310107693650659740228e-18,
  76  -7.5225243218253901727164675011e-20
  77 };
  78 static cheb_series g2_cs = {
  79   g2_dat,
  80   14,
  81   -1, 1,
  82   8
  83 };
  84
  85
  86 static
  87 int
  88 gsl_sf_temme_gamma(const double nu, double * g_1pnu, double * g_1mnu, double * g1, double * g2)
  89 {
  90   const double anu = fabs(nu);    /* functions are even */
  91   const double x = 4.0*anu - 1.0;
  92   gsl_sf_result r_g1;
  93   gsl_sf_result r_g2;
  94   cheb_eval_e(&g1_cs, x, &r_g1);
  95   cheb_eval_e(&g2_cs, x, &r_g2);
  96   *g1 = r_g1.val;
  97   *g2 = r_g2.val;
  98   *g_1mnu = 1.0/(r_g2.val + nu * r_g1.val);
  99   *g_1pnu = 1.0/(r_g2.val - nu * r_g1.val);
 100   return GSL_SUCCESS;
 101 }
 102
 103
 104 int
 105 gsl_sf_bessel_Y_temme(const double nu, const double x,
 106                       gsl_sf_result * Ynu,
 107                       gsl_sf_result * Ynup1)
 108 {
 109   const int max_iter = 15000;
 110
 111   const double half_x = 0.5 * x;
 112   const double ln_half_x = log(half_x);
 113   const double half_x_nu = exp(nu*ln_half_x);
 114   const double pi_nu   = M_PI * nu;
 115   const double alpha   = pi_nu / 2.0;
 116   const double sigma   = -nu * ln_half_x;
 117   const double sinrat  = (fabs(pi_nu) < GSL_DBL_EPSILON ? 1.0 : pi_nu/sin(pi_nu));
 118   const double sinhrat = (fabs(sigma) < GSL_DBL_EPSILON ? 1.0 : sinh(sigma)/sigma);
 119   const double sinhalf = (fabs(alpha) < GSL_DBL_EPSILON ? 1.0 : sin(alpha)/alpha);
 120   const double sin_sqr = nu*M_PI*M_PI*0.5 * sinhalf*sinhalf;
 121
 122   double sum0, sum1;
 123   double fk, pk, qk, hk, ck;
 124   int k = 0;
 125   int stat_iter;
 126
 127   double g_1pnu, g_1mnu, g1, g2;
 128   int stat_g = gsl_sf_temme_gamma(nu, &g_1pnu, &g_1mnu, &g1, &g2);
 129
 130   fk = 2.0/M_PI * sinrat * (cosh(sigma)*g1 - sinhrat*ln_half_x*g2);
 131   pk = 1.0/M_PI /half_x_nu * g_1pnu;
 132   qk = 1.0/M_PI *half_x_nu * g_1mnu;
 133   hk = pk;
 134   ck = 1.0;
 135
 136   sum0 = fk + sin_sqr * qk;
 137   sum1 = pk;
 138
 139   while(k < max_iter) {
 140     double del0;
 141     double del1;
 142     double gk;
 143     k++;
 144     fk  = (k*fk + pk + qk)/(k*k-nu*nu);
 145     ck *= -half_x*half_x/k;
 146     pk /= (k - nu);
 147     qk /= (k + nu);
 148     gk  = fk + sin_sqr * qk;
 149     hk  = -k*gk + pk;
 150     del0 = ck * gk;
 151     del1 = ck * hk;
 152     sum0 += del0;
 153     sum1 += del1;
 154     if(fabs(del0) < 0.5*(1.0 + fabs(sum0))*GSL_DBL_EPSILON) break;
 155   }
 156
 157   Ynu->val   = -sum0;
 158   Ynu->err   = (2.0 + 0.5*k) * GSL_DBL_EPSILON * fabs(Ynu->val);
 159   Ynup1->val = -sum1 * 2.0/x;
 160   Ynup1->err = (2.0 + 0.5*k) * GSL_DBL_EPSILON * fabs(Ynup1->val);
 161
 162   stat_iter = ( k >= max_iter ? GSL_EMAXITER : GSL_SUCCESS );
 163   return GSL_ERROR_SELECT_2(stat_iter, stat_g);
 164 }
 165
 166
 167 int
 168 gsl_sf_bessel_K_scaled_temme(const double nu, const double x,
 169                              double * K_nu, double * K_nup1, double * Kp_nu)
 170 {
 171   const int max_iter = 15000;
 172
 173   const double half_x    = 0.5 * x;
 174   const double ln_half_x = log(half_x);
 175   const double half_x_nu = exp(nu*ln_half_x);
 176   const double pi_nu   = M_PI * nu;
 177   const double sigma   = -nu * ln_half_x;
 178   const double sinrat  = (fabs(pi_nu) < GSL_DBL_EPSILON ? 1.0 : pi_nu/sin(pi_nu));
 179   const double sinhrat = (fabs(sigma) < GSL_DBL_EPSILON ? 1.0 : sinh(sigma)/sigma);
 180   const double ex = exp(x);
 181
 182   double sum0, sum1;
 183   double fk, pk, qk, hk, ck;
 184   int k = 0;
 185   int stat_iter;
 186
 187   double g_1pnu, g_1mnu, g1, g2;
 188   int stat_g = gsl_sf_temme_gamma(nu, &g_1pnu, &g_1mnu, &g1, &g2);
 189
 190   fk = sinrat * (cosh(sigma)*g1 - sinhrat*ln_half_x*g2);
 191   pk = 0.5/half_x_nu * g_1pnu;
 192   qk = 0.5*half_x_nu * g_1mnu;
 193   hk = pk;
 194   ck = 1.0;
 195   sum0 = fk;
 196   sum1 = hk;
 197   while(k < max_iter) {
 198     double del0;
 199     double del1;
 200     k++;
 201     fk  = (k*fk + pk + qk)/(k*k-nu*nu);
 202     ck *= half_x*half_x/k;
 203     pk /= (k - nu);
 204     qk /= (k + nu);
 205     hk  = -k*fk + pk;
 206     del0 = ck * fk;
 207     del1 = ck * hk;
 208     sum0 += del0;
 209     sum1 += del1;
 210     if(fabs(del0) < 0.5*fabs(sum0)*GSL_DBL_EPSILON) break;
 211   }
 212
 213   *K_nu   = sum0 * ex;
 214   *K_nup1 = sum1 * 2.0/x * ex;
 215   *Kp_nu  = - *K_nup1 + nu/x * *K_nu;
 216
 217   stat_iter = ( k == max_iter ? GSL_EMAXITER : GSL_SUCCESS );
 218   return GSL_ERROR_SELECT_2(stat_iter, stat_g);
 219 }