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

   1 /* specfunc/bessel.c
   2  *
   3  * Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 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 /* Miscellaneous support functions for Bessel function evaluations.
  22  */
  23 #include <config.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_airy.h>
  27 #include <gsl/gsl_sf_elementary.h>
  28 #include <gsl/gsl_sf_exp.h>
  29 #include <gsl/gsl_sf_gamma.h>
  30 #include <gsl/gsl_sf_trig.h>
  31
  32 #include "error.h"
  33
  34 #include "bessel_amp_phase.h"
  35 #include "bessel_temme.h"
  36 #include "bessel.h"
  37
  38 #define CubeRoot2_  1.25992104989487316476721060728
  39
  40
  41
  42 /* Debye functions [Abramowitz+Stegun, 9.3.9-10] */
  43
  44 inline static double
  45 debye_u1(const double * tpow)
  46 {
  47   return (3.0*tpow[1] - 5.0*tpow[3])/24.0;
  48 }
  49
  50 inline static double
  51 debye_u2(const double * tpow)
  52 {
  53   return (81.0*tpow[2] - 462.0*tpow[4] + 385.0*tpow[6])/1152.0;
  54 }
  55
  56 inline
  57 static double debye_u3(const double * tpow)
  58 {
  59   return (30375.0*tpow[3] - 369603.0*tpow[5] + 765765.0*tpow[7] - 425425.0*tpow[9])/414720.0;
  60 }
  61
  62 inline
  63 static double debye_u4(const double * tpow)
  64 {
  65   return (4465125.0*tpow[4] - 94121676.0*tpow[6] + 349922430.0*tpow[8] -
  66           446185740.0*tpow[10] + 185910725.0*tpow[12])/39813120.0;
  67 }
  68
  69 inline
  70 static double debye_u5(const double * tpow)
  71 {
  72   return (1519035525.0*tpow[5]     - 49286948607.0*tpow[7] +
  73           284499769554.0*tpow[9]   - 614135872350.0*tpow[11] +
  74           566098157625.0*tpow[13]  - 188699385875.0*tpow[15])/6688604160.0;
  75 }
  76
  77 #if 0
  78 inline
  79 static double debye_u6(const double * tpow)
  80 {
  81   return (2757049477875.0*tpow[6] - 127577298354750.0*tpow[8] +
  82           1050760774457901.0*tpow[10] - 3369032068261860.0*tpow[12] +
  83           5104696716244125.0*tpow[14] - 3685299006138750.0*tpow[16] +
  84           1023694168371875.0*tpow[18])/4815794995200.0;
  85 }
  86 #endif
  87
  88
  89 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  90
  91 int
  92 gsl_sf_bessel_IJ_taylor_e(const double nu, const double x,
  93                              const int sign,
  94                              const int kmax,
  95                              const double threshold,
  96                              gsl_sf_result * result
  97                              )
  98 {
  99   /* CHECK_POINTER(result) */
 100
 101   if(nu < 0.0 || x < 0.0) {
 102     DOMAIN_ERROR(result);
 103   }
 104   else if(x == 0.0) {
 105     if(nu == 0.0) {
 106       result->val = 1.0;
 107       result->err = 0.0;
 108     }
 109     else {
 110       result->val = 0.0;
 111       result->err = 0.0;
 112     }
 113     return GSL_SUCCESS;
 114   }
 115   else {
 116     gsl_sf_result prefactor;   /* (x/2)^nu / Gamma(nu+1) */
 117     gsl_sf_result sum;
 118
 119     int stat_pre;
 120     int stat_sum;
 121     int stat_mul;
 122
 123     if(nu == 0.0) {
 124       prefactor.val = 1.0;
 125       prefactor.err = 0.0;
 126       stat_pre = GSL_SUCCESS;
 127     }
 128     else if(nu < INT_MAX-1) {
 129       /* Separate the integer part and use
 130        * y^nu / Gamma(nu+1) = y^N /N! y^f / (N+1)_f,
 131        * to control the error.
 132        */
 133       const int    N = (int)floor(nu + 0.5);
 134       const double f = nu - N;
 135       gsl_sf_result poch_factor;
 136       gsl_sf_result tc_factor;
 137       const int stat_poch = gsl_sf_poch_e(N+1.0, f, &poch_factor);
 138       const int stat_tc   = gsl_sf_taylorcoeff_e(N, 0.5*x, &tc_factor);
 139       const double p = pow(0.5*x,f);
 140       prefactor.val  = tc_factor.val * p / poch_factor.val;
 141       prefactor.err  = tc_factor.err * p / poch_factor.val;
 142       prefactor.err += fabs(prefactor.val) / poch_factor.val * poch_factor.err;
 143       prefactor.err += 2.0 * GSL_DBL_EPSILON * fabs(prefactor.val);
 144       stat_pre = GSL_ERROR_SELECT_2(stat_tc, stat_poch);
 145     }
 146     else {
 147       gsl_sf_result lg;
 148       const int stat_lg = gsl_sf_lngamma_e(nu+1.0, &lg);
 149       const double term1  = nu*log(0.5*x);
 150       const double term2  = lg.val;
 151       const double ln_pre = term1 - term2;
 152       const double ln_pre_err = GSL_DBL_EPSILON * (fabs(term1)+fabs(term2)) + lg.err;
 153       const int stat_ex = gsl_sf_exp_err_e(ln_pre, ln_pre_err, &prefactor);
 154       stat_pre = GSL_ERROR_SELECT_2(stat_ex, stat_lg);
 155     }
 156
 157     /* Evaluate the sum.
 158      * [Abramowitz+Stegun, 9.1.10]
 159      * [Abramowitz+Stegun, 9.6.7]
 160      */
 161     {
 162       const double y = sign * 0.25 * x*x;
 163       double sumk = 1.0;
 164       double term = 1.0;
 165       int k;
 166
 167       for(k=1; k<=kmax; k++) {
 168         term *= y/((nu+k)*k);
 169         sumk += term;
 170         if(fabs(term/sumk) < threshold) break;
 171       }
 172
 173       sum.val = sumk;
 174       sum.err = threshold * fabs(sumk);
 175
 176       stat_sum = ( k >= kmax ? GSL_EMAXITER : GSL_SUCCESS );
 177     }
 178
 179     stat_mul = gsl_sf_multiply_err_e(prefactor.val, prefactor.err,
 180                                         sum.val, sum.err,
 181                                         result);
 182
 183     return GSL_ERROR_SELECT_3(stat_mul, stat_pre, stat_sum);
 184   }
 185 }
 186
 187
 188 /* Hankel's Asymptotic Expansion - A&S 9.2.5
 189  *
 190  * x >> nu*nu+1
 191  * error ~ O( ((nu*nu+1)/x)^4 )
 192  *
 193  * empirical error analysis:
 194  *   choose  GSL_ROOT4_MACH_EPS * x > (nu*nu + 1)
 195  *
 196  * This is not especially useful. When the argument gets
 197  * large enough for this to apply, the cos() and sin()
 198  * start loosing digits. However, this seems inevitable
 199  * for this particular method.
 200  *
 201  * Wed Jun 25 14:39:38 MDT 2003 [GJ]
 202  * This function was inconsistent since the Q term did not
 203  * go to relative order eps^2. That's why the error estimate
 204  * originally given was screwy (it didn't make sense that the
 205  * "empirical" error was coming out O(eps^3)).
 206  * With Q to proper order, the error is O(eps^4).
 207  *
 208  * Sat Mar 15 05:16:18 GMT 2008 [BG]
 209  * Extended to use additional terms in the series to gain
 210  * higher accuracy.
 211  *
 212  */
 213
 214 int
 215 gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, gsl_sf_result * result)
 216 {
 217   double mu  = 4.0*nu*nu;
 218   double chi = x - (0.5*nu + 0.25)*M_PI;
 219
 220   double P   = 0.0;
 221   double Q   = 0.0;
 222
 223   double k = 0, t = 1;
 224   int convP, convQ;
 225
 226   do
 227     {
 228       t *= (k == 0) ? 1 : -(mu - (2*k-1)*(2*k-1)) / (k * (8 * x));
 229       convP = fabs(t) < GSL_DBL_EPSILON * fabs(P);
 230       P += t;
 231
 232       k++;
 233
 234       t *= (mu - (2*k-1)*(2*k-1)) / (k * (8 * x));
 235       convQ = fabs(t) < GSL_DBL_EPSILON * fabs(Q);
 236       Q += t;
 237
 238       /* To preserve the consistency of the series we need to exit
 239          when P and Q have the same number of terms */
 240
 241       if (convP && convQ && k > (nu / 2))
 242         break;
 243
 244       k++;
 245     }
 246   while (k < 1000);
 247
 248   {
 249     double pre = sqrt(2.0/(M_PI*x));
 250     double c   = cos(chi);
 251     double s   = sin(chi);
 252
 253     result->val  = pre * (c*P - s*Q);
 254     result->err  = pre * GSL_DBL_EPSILON * (fabs(c*P) + fabs(s*Q) + fabs(t)) * (1 + fabs(x));
 255     /* NB: final term accounts for phase error with large x */
 256   }
 257
 258   return GSL_SUCCESS;
 259 }
 260
 261
 262 /* x >> nu*nu+1
 263  */
 264 int
 265 gsl_sf_bessel_Ynu_asympx_e(const double nu, const double x, gsl_sf_result * result)
 266 {
 267   double ampl;
 268   double theta;
 269   double alpha = x;
 270   double beta  = -0.5*nu*M_PI;
 271   int stat_a = gsl_sf_bessel_asymp_Mnu_e(nu, x, &ampl);
 272   int stat_t = gsl_sf_bessel_asymp_thetanu_corr_e(nu, x, &theta);
 273   double sin_alpha = sin(alpha);
 274   double cos_alpha = cos(alpha);
 275   double sin_chi   = sin(beta + theta);
 276   double cos_chi   = cos(beta + theta);
 277   double sin_term     = sin_alpha * cos_chi + sin_chi * cos_alpha;
 278   double sin_term_mag = fabs(sin_alpha * cos_chi) + fabs(sin_chi * cos_alpha);
 279   result->val  = ampl * sin_term;
 280   result->err  = fabs(ampl) * GSL_DBL_EPSILON * sin_term_mag;
 281   result->err += fabs(result->val) * 2.0 * GSL_DBL_EPSILON;
 282
 283   if(fabs(alpha) > 1.0/GSL_DBL_EPSILON) {
 284     result->err *= 0.5 * fabs(alpha);
 285   }
 286   else if(fabs(alpha) > 1.0/GSL_SQRT_DBL_EPSILON) {
 287     result->err *= 256.0 * fabs(alpha) * GSL_SQRT_DBL_EPSILON;
 288   }
 289
 290   return GSL_ERROR_SELECT_2(stat_t, stat_a);
 291 }
 292
 293
 294 /* x >> nu*nu+1
 295  */
 296 int
 297 gsl_sf_bessel_Inu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)
 298 {
 299   double mu   = 4.0*nu*nu;
 300   double mum1 = mu-1.0;
 301   double mum9 = mu-9.0;
 302   double pre  = 1.0/sqrt(2.0*M_PI*x);
 303   double r    = mu/x;
 304   result->val = pre * (1.0 - mum1/(8.0*x) + mum1*mum9/(128.0*x*x));
 305   result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * fabs(0.1*r*r*r);
 306   return GSL_SUCCESS;
 307 }
 308
 309 /* x >> nu*nu+1
 310  */
 311 int
 312 gsl_sf_bessel_Knu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)
 313 {
 314   double mu   = 4.0*nu*nu;
 315   double mum1 = mu-1.0;
 316   double mum9 = mu-9.0;
 317   double pre  = sqrt(M_PI/(2.0*x));
 318   double r    = nu/x;
 319   result->val = pre * (1.0 + mum1/(8.0*x) + mum1*mum9/(128.0*x*x));
 320   result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val) + pre * fabs(0.1*r*r*r);
 321   return GSL_SUCCESS;
 322 }
 323
 324
 325 /* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.7]
 326  *
 327  * error:
 328  *   The error has the form u_N(t)/nu^N  where  0 <= t <= 1.
 329  *   It is not hard to show that |u_N(t)| is small for such t.
 330  *   We have N=6 here, and |u_6(t)| < 0.025, so the error is clearly
 331  *   bounded by 0.025/nu^6. This gives the asymptotic bound on nu
 332  *   seen below as nu ~ 100. For general MACH_EPS it will be
 333  *                     nu > 0.5 / MACH_EPS^(1/6)
 334  *   When t is small, the bound is even better because |u_N(t)| vanishes
 335  *   as t->0. In fact u_N(t) ~ C t^N as t->0, with C ~= 0.1.
 336  *   We write
 337  *                     err_N <= min(0.025, C(1/(1+(x/nu)^2))^3) / nu^6
 338  *   therefore
 339  *                     min(0.29/nu^2, 0.5/(nu^2+x^2)) < MACH_EPS^{1/3}
 340  *   and this is the general form.
 341  *
 342  * empirical error analysis, assuming 14 digit requirement:
 343  *   choose   x > 50.000 nu   ==>  nu >   3
 344  *   choose   x > 10.000 nu   ==>  nu >  15
 345  *   choose   x >  2.000 nu   ==>  nu >  50
 346  *   choose   x >  1.000 nu   ==>  nu >  75
 347  *   choose   x >  0.500 nu   ==>  nu >  80
 348  *   choose   x >  0.100 nu   ==>  nu >  83
 349  *
 350  * This makes sense. For x << nu, the error will be of the form u_N(1)/nu^N,
 351  * since the polynomial term will be evaluated near t=1, so the bound
 352  * on nu will become constant for small x. Furthermore, increasing x with
 353  * nu fixed will decrease the error.
 354  */
 355 int
 356 gsl_sf_bessel_Inu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)
 357 {
 358   int i;
 359   double z = x/nu;
 360   double root_term = hypot(1.0,z);
 361   double pre = 1.0/sqrt(2.0*M_PI*nu * root_term);
 362   double eta = root_term + log(z/(1.0+root_term));
 363   double ex_arg = ( z < 1.0/GSL_ROOT3_DBL_EPSILON ? nu*(-z + eta) : -0.5*nu/z*(1.0 - 1.0/(12.0*z*z)) );
 364   gsl_sf_result ex_result;
 365   int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);
 366   if(stat_ex == GSL_SUCCESS) {
 367     double t = 1.0/root_term;
 368     double sum;
 369     double tpow[16];
 370     tpow[0] = 1.0;
 371     for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];
 372     sum = 1.0 + debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) + debye_u3(tpow)/(nu*nu*nu)
 373           + debye_u4(tpow)/(nu*nu*nu*nu) + debye_u5(tpow)/(nu*nu*nu*nu*nu);
 374     result->val  = pre * ex_result.val * sum;
 375     result->err  = pre * ex_result.val / (nu*nu*nu*nu*nu*nu);
 376     result->err += pre * ex_result.err * fabs(sum);
 377     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 378     return GSL_SUCCESS;
 379   }
 380   else {
 381     result->val = 0.0;
 382     result->err = 0.0;
 383     return stat_ex;
 384   }
 385 }
 386
 387
 388 /* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.8]
 389  *
 390  * error:
 391  *   identical to that above for Inu_scaled
 392  */
 393 int
 394 gsl_sf_bessel_Knu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)
 395 {
 396   int i;
 397   double z = x/nu;
 398   double root_term = hypot(1.0,z);
 399   double pre = sqrt(M_PI/(2.0*nu*root_term));
 400   double eta = root_term + log(z/(1.0+root_term));
 401   double ex_arg = ( z < 1.0/GSL_ROOT3_DBL_EPSILON ? nu*(z - eta) : 0.5*nu/z*(1.0 + 1.0/(12.0*z*z)) );
 402   gsl_sf_result ex_result;
 403   int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);
 404   if(stat_ex == GSL_SUCCESS) {
 405     double t = 1.0/root_term;
 406     double sum;
 407     double tpow[16];
 408     tpow[0] = 1.0;
 409     for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];
 410     sum = 1.0 - debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) - debye_u3(tpow)/(nu*nu*nu)
 411           + debye_u4(tpow)/(nu*nu*nu*nu) - debye_u5(tpow)/(nu*nu*nu*nu*nu);
 412     result->val  = pre * ex_result.val * sum;
 413     result->err  = pre * ex_result.err * fabs(sum);
 414     result->err += pre * ex_result.val / (nu*nu*nu*nu*nu*nu);
 415     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 416     return GSL_SUCCESS;
 417   }
 418   else {
 419     result->val = 0.0;
 420     result->err = 0.0;
 421     return stat_ex;
 422   }
 423 }
 424
 425
 426 /* Evaluate J_mu(x),J_{mu+1}(x) and Y_mu(x),Y_{mu+1}(x)  for |mu| < 1/2
 427  */
 428 int
 429 gsl_sf_bessel_JY_mu_restricted(const double mu, const double x,
 430                                gsl_sf_result * Jmu, gsl_sf_result * Jmup1,
 431                                gsl_sf_result * Ymu, gsl_sf_result * Ymup1)
 432 {
 433   /* CHECK_POINTER(Jmu) */
 434   /* CHECK_POINTER(Jmup1) */
 435   /* CHECK_POINTER(Ymu) */
 436   /* CHECK_POINTER(Ymup1) */
 437
 438   if(x < 0.0 || fabs(mu) > 0.5) {
 439     Jmu->val   = 0.0;
 440     Jmu->err   = 0.0;
 441     Jmup1->val = 0.0;
 442     Jmup1->err = 0.0;
 443     Ymu->val   = 0.0;
 444     Ymu->err   = 0.0;
 445     Ymup1->val = 0.0;
 446     Ymup1->err = 0.0;
 447     GSL_ERROR ("error", GSL_EDOM);
 448   }
 449   else if(x == 0.0) {
 450     if(mu == 0.0) {
 451       Jmu->val   = 1.0;
 452       Jmu->err   = 0.0;
 453     }
 454     else {
 455       Jmu->val   = 0.0;
 456       Jmu->err   = 0.0;
 457     }
 458     Jmup1->val = 0.0;
 459     Jmup1->err = 0.0;
 460     Ymu->val   = 0.0;
 461     Ymu->err   = 0.0;
 462     Ymup1->val = 0.0;
 463     Ymup1->err = 0.0;
 464     GSL_ERROR ("error", GSL_EDOM);
 465   }
 466   else {
 467     int stat_Y;
 468     int stat_J;
 469
 470     if(x < 2.0) {
 471       /* Use Taylor series for J and the Temme series for Y.
 472        * The Taylor series for J requires nu > 0, so we shift
 473        * up one and use the recursion relation to get Jmu, in
 474        * case mu < 0.
 475        */
 476       gsl_sf_result Jmup2;
 477       int stat_J1 = gsl_sf_bessel_IJ_taylor_e(mu+1.0, x, -1, 100, GSL_DBL_EPSILON,  Jmup1);
 478       int stat_J2 = gsl_sf_bessel_IJ_taylor_e(mu+2.0, x, -1, 100, GSL_DBL_EPSILON, &Jmup2);
 479       double c = 2.0*(mu+1.0)/x;
 480       Jmu->val  = c * Jmup1->val - Jmup2.val;
 481       Jmu->err  = c * Jmup1->err + Jmup2.err;
 482       Jmu->err += 2.0 * GSL_DBL_EPSILON * fabs(Jmu->val);
 483       stat_J = GSL_ERROR_SELECT_2(stat_J1, stat_J2);
 484       stat_Y = gsl_sf_bessel_Y_temme(mu, x, Ymu, Ymup1);
 485       return GSL_ERROR_SELECT_2(stat_J, stat_Y);
 486     }
 487     else if(x < 1000.0) {
 488       double P, Q;
 489       double J_ratio;
 490       double J_sgn;
 491       const int stat_CF1 = gsl_sf_bessel_J_CF1(mu, x, &J_ratio, &J_sgn);
 492       const int stat_CF2 = gsl_sf_bessel_JY_steed_CF2(mu, x, &P, &Q);
 493       double Jprime_J_ratio = mu/x - J_ratio;
 494       double gamma = (P - Jprime_J_ratio)/Q;
 495       Jmu->val = J_sgn * sqrt(2.0/(M_PI*x) / (Q + gamma*(P-Jprime_J_ratio)));
 496       Jmu->err = 4.0 * GSL_DBL_EPSILON * fabs(Jmu->val);
 497       Jmup1->val = J_ratio * Jmu->val;
 498       Jmup1->err = fabs(J_ratio) * Jmu->err;
 499       Ymu->val = gamma * Jmu->val;
 500       Ymu->err = fabs(gamma) * Jmu->err;
 501       Ymup1->val = Ymu->val * (mu/x - P - Q/gamma);
 502       Ymup1->err = Ymu->err * fabs(mu/x - P - Q/gamma) + 4.0*GSL_DBL_EPSILON*fabs(Ymup1->val);
 503       return GSL_ERROR_SELECT_2(stat_CF1, stat_CF2);
 504     }
 505     else {
 506       /* Use asymptotics for large argument.
 507        */
 508       const int stat_J0 = gsl_sf_bessel_Jnu_asympx_e(mu,     x, Jmu);
 509       const int stat_J1 = gsl_sf_bessel_Jnu_asympx_e(mu+1.0, x, Jmup1);
 510       const int stat_Y0 = gsl_sf_bessel_Ynu_asympx_e(mu,     x, Ymu);
 511       const int stat_Y1 = gsl_sf_bessel_Ynu_asympx_e(mu+1.0, x, Ymup1);
 512       stat_J = GSL_ERROR_SELECT_2(stat_J0, stat_J1);
 513       stat_Y = GSL_ERROR_SELECT_2(stat_Y0, stat_Y1);
 514       return GSL_ERROR_SELECT_2(stat_J, stat_Y);
 515     }
 516   }
 517 }
 518
 519
 520 int
 521 gsl_sf_bessel_J_CF1(const double nu, const double x,
 522                     double * ratio, double * sgn)
 523 {
 524   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
 525   const double RECUR_SMALL = GSL_SQRT_DBL_MIN;
 526   const int maxiter = 10000;
 527   int n = 1;
 528   double Anm2 = 1.0;
 529   double Bnm2 = 0.0;
 530   double Anm1 = 0.0;
 531   double Bnm1 = 1.0;
 532   double a1 = x/(2.0*(nu+1.0));
 533   double An = Anm1 + a1*Anm2;
 534   double Bn = Bnm1 + a1*Bnm2;
 535   double an;
 536   double fn = An/Bn;
 537   double dn = a1;
 538   double s  = 1.0;
 539
 540   while(n < maxiter) {
 541     double old_fn;
 542     double del;
 543     n++;
 544     Anm2 = Anm1;
 545     Bnm2 = Bnm1;
 546     Anm1 = An;
 547     Bnm1 = Bn;
 548     an = -x*x/(4.0*(nu+n-1.0)*(nu+n));
 549     An = Anm1 + an*Anm2;
 550     Bn = Bnm1 + an*Bnm2;
 551
 552     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
 553       An /= RECUR_BIG;
 554       Bn /= RECUR_BIG;
 555       Anm1 /= RECUR_BIG;
 556       Bnm1 /= RECUR_BIG;
 557       Anm2 /= RECUR_BIG;
 558     } else if(fabs(An) < RECUR_SMALL || fabs(Bn) < RECUR_SMALL) {
 559       An /= RECUR_SMALL;
 560       Bn /= RECUR_SMALL;
 561       Anm1 /= RECUR_SMALL;
 562       Bnm1 /= RECUR_SMALL;
 563       Anm2 /= RECUR_SMALL;
 564       Bnm2 /= RECUR_SMALL;
 565     }
 566
 567     old_fn = fn;
 568     fn = An/Bn;
 569     del = old_fn/fn;
 570
 571     dn = 1.0 / (2.0*(nu+n)/x - dn);
 572     if(dn < 0.0) s = -s;
 573
 574     if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
 575   }
 576
 577   /* FIXME: we should return an error term here as well, because the
 578      error from this recurrence affects the overall error estimate. */
 579
 580   *ratio = fn;
 581   *sgn   = s;
 582
 583   if(n >= maxiter)
 584     GSL_ERROR ("error", GSL_EMAXITER);
 585   else
 586     return GSL_SUCCESS;
 587 }
 588
 589
 590
 591 /* Evaluate the continued fraction CF1 for J_{nu+1}/J_nu
 592  * using Gautschi (Euler) equivalent series.
 593  * This exhibits an annoying problem because the
 594  * a_k are not positive definite (in fact they are all negative).
 595  * There are cases when rho_k blows up. Example: nu=1,x=4.
 596  */
 597 #if 0
 598 int
 599 gsl_sf_bessel_J_CF1_ser(const double nu, const double x,
 600                         double * ratio, double * sgn)
 601 {
 602   const int maxk = 20000;
 603   double tk   = 1.0;
 604   double sum  = 1.0;
 605   double rhok = 0.0;
 606   double dk = 0.0;
 607   double s  = 1.0;
 608   int k;
 609
 610   for(k=1; k<maxk; k++) {
 611     double ak = -0.25 * (x/(nu+k)) * x/(nu+k+1.0);
 612     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
 613     tk  *= rhok;
 614     sum += tk;
 615
 616     dk = 1.0 / (2.0/x - (nu+k-1.0)/(nu+k) * dk);
 617     if(dk < 0.0) s = -s;
 618
 619     if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
 620   }
 621
 622   *ratio = x/(2.0*(nu+1.0)) * sum;
 623   *sgn   = s;
 624
 625   if(k == maxk)
 626     GSL_ERROR ("error", GSL_EMAXITER);
 627   else
 628     return GSL_SUCCESS;
 629 }
 630 #endif
 631
 632
 633 /* Evaluate the continued fraction CF1 for I_{nu+1}/I_nu
 634  * using Gautschi (Euler) equivalent series.
 635  */
 636 int
 637 gsl_sf_bessel_I_CF1_ser(const double nu, const double x, double * ratio)
 638 {
 639   const int maxk = 20000;
 640   double tk   = 1.0;
 641   double sum  = 1.0;
 642   double rhok = 0.0;
 643   int k;
 644
 645   for(k=1; k<maxk; k++) {
 646     double ak = 0.25 * (x/(nu+k)) * x/(nu+k+1.0);
 647     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
 648     tk  *= rhok;
 649     sum += tk;
 650     if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
 651   }
 652
 653   *ratio = x/(2.0*(nu+1.0)) * sum;
 654
 655   if(k == maxk)
 656     GSL_ERROR ("error", GSL_EMAXITER);
 657   else
 658     return GSL_SUCCESS;
 659 }
 660
 661
 662 int
 663 gsl_sf_bessel_JY_steed_CF2(const double nu, const double x,
 664                            double * P, double * Q)
 665 {
 666   const int max_iter = 10000;
 667   const double SMALL = 1.0e-100;
 668
 669   int i = 1;
 670
 671   double x_inv = 1.0/x;
 672   double a = 0.25 - nu*nu;
 673   double p = -0.5*x_inv;
 674   double q = 1.0;
 675   double br = 2.0*x;
 676   double bi = 2.0;
 677   double fact = a*x_inv/(p*p + q*q);
 678   double cr = br + q*fact;
 679   double ci = bi + p*fact;
 680   double den = br*br + bi*bi;
 681   double dr = br/den;
 682   double di = -bi/den;
 683   double dlr = cr*dr - ci*di;
 684   double dli = cr*di + ci*dr;
 685   double temp = p*dlr - q*dli;
 686   q = p*dli + q*dlr;
 687   p = temp;
 688   for (i=2; i<=max_iter; i++) {
 689     a  += 2*(i-1);
 690     bi += 2.0;
 691     dr = a*dr + br;
 692     di = a*di + bi;
 693     if(fabs(dr)+fabs(di) < SMALL) dr = SMALL;
 694     fact = a/(cr*cr+ci*ci);
 695     cr = br + cr*fact;
 696     ci = bi - ci*fact;
 697     if(fabs(cr)+fabs(ci) < SMALL) cr = SMALL;
 698     den = dr*dr + di*di;
 699     dr /= den;
 700     di /= -den;
 701     dlr = cr*dr - ci*di;
 702     dli = cr*di + ci*dr;
 703     temp = p*dlr - q*dli;
 704     q = p*dli + q*dlr;
 705     p = temp;
 706     if(fabs(dlr-1.0)+fabs(dli) < GSL_DBL_EPSILON) break;
 707   }
 708
 709   *P = p;
 710   *Q = q;
 711
 712   if(i == max_iter)
 713     GSL_ERROR ("error", GSL_EMAXITER);
 714   else
 715     return GSL_SUCCESS;
 716 }
 717
 718
 719 /* Evaluate continued fraction CF2, using Thompson-Barnett-Temme method,
 720  * to obtain values of exp(x)*K_nu and exp(x)*K_{nu+1}.
 721  *
 722  * This is unstable for small x; x > 2 is a good cutoff.
 723  * Also requires |nu| < 1/2.
 724  */
 725 int
 726 gsl_sf_bessel_K_scaled_steed_temme_CF2(const double nu, const double x,
 727                                        double * K_nu, double * K_nup1,
 728                                        double * Kp_nu)
 729 {
 730   const int maxiter = 10000;
 731
 732   int i = 1;
 733   double bi = 2.0*(1.0 + x);
 734   double di = 1.0/bi;
 735   double delhi = di;
 736   double hi    = di;
 737
 738   double qi   = 0.0;
 739   double qip1 = 1.0;
 740
 741   double ai = -(0.25 - nu*nu);
 742   double a1 = ai;
 743   double ci = -ai;
 744   double Qi = -ai;
 745
 746   double s = 1.0 + Qi*delhi;
 747
 748   for(i=2; i<=maxiter; i++) {
 749     double dels;
 750     double tmp;
 751     ai -= 2.0*(i-1);
 752     ci  = -ai*ci/i;
 753     tmp  = (qi - bi*qip1)/ai;
 754     qi   = qip1;
 755     qip1 = tmp;
 756     Qi += ci*qip1;
 757     bi += 2.0;
 758     di  = 1.0/(bi + ai*di);
 759     delhi = (bi*di - 1.0) * delhi;
 760     hi += delhi;
 761     dels = Qi*delhi;
 762     s += dels;
 763     if(fabs(dels/s) < GSL_DBL_EPSILON) break;
 764   }
 765
 766   hi *= -a1;
 767
 768   *K_nu   = sqrt(M_PI/(2.0*x)) / s;
 769   *K_nup1 = *K_nu * (nu + x + 0.5 - hi)/x;
 770   *Kp_nu  = - *K_nup1 + nu/x * *K_nu;
 771   if(i == maxiter)
 772     GSL_ERROR ("error", GSL_EMAXITER);
 773   else
 774     return GSL_SUCCESS;
 775 }
 776
 777
 778 int gsl_sf_bessel_cos_pi4_e(double y, double eps, gsl_sf_result * result)
 779 {
 780   const double sy = sin(y);
 781   const double cy = cos(y);
 782   const double s = sy + cy;
 783   const double d = sy - cy;
 784   const double abs_sum = fabs(cy) + fabs(sy);
 785   double seps;
 786   double ceps;
 787   if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
 788     const double e2 = eps*eps;
 789     seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));
 790     ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);
 791   }
 792   else {
 793     seps = sin(eps);
 794     ceps = cos(eps);
 795   }
 796   result->val = (ceps * s - seps * d)/ M_SQRT2;
 797   result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;
 798
 799   /* Try to account for error in evaluation of sin(y), cos(y).
 800    * This is a little sticky because we don't really know
 801    * how the library routines are doing their argument reduction.
 802    * However, we will make a reasonable guess.
 803    * FIXME ?
 804    */
 805   if(y > 1.0/GSL_DBL_EPSILON) {
 806     result->err *= 0.5 * y;
 807   }
 808   else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {
 809     result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;
 810   }
 811
 812   return GSL_SUCCESS;
 813 }
 814
 815
 816 int gsl_sf_bessel_sin_pi4_e(double y, double eps, gsl_sf_result * result)
 817 {
 818   const double sy = sin(y);
 819   const double cy = cos(y);
 820   const double s = sy + cy;
 821   const double d = sy - cy;
 822   const double abs_sum = fabs(cy) + fabs(sy);
 823   double seps;
 824   double ceps;
 825   if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
 826     const double e2 = eps*eps;
 827     seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));
 828     ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);
 829   }
 830   else {
 831     seps = sin(eps);
 832     ceps = cos(eps);
 833   }
 834   result->val = (ceps * d + seps * s)/ M_SQRT2;
 835   result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;
 836
 837   /* Try to account for error in evaluation of sin(y), cos(y).
 838    * See above.
 839    * FIXME ?
 840    */
 841   if(y > 1.0/GSL_DBL_EPSILON) {
 842     result->err *= 0.5 * y;
 843   }
 844   else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {
 845     result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;
 846   }
 847
 848   return GSL_SUCCESS;
 849 }
 850
 851
 852 /************************************************************************
 853  *                                                                      *
 854   Asymptotic approximations 8.11.5, 8.12.5, and 8.42.7 from
 855   G.N.Watson, A Treatise on the Theory of Bessel Functions,
 856   2nd Edition (Cambridge University Press, 1944).
 857   Higher terms in expansion for x near l given by
 858   Airey in Phil. Mag. 31, 520 (1916).
 859
 860   This approximation is accurate to near 0.1% at the boundaries
 861   between the asymptotic regions; well away from the boundaries
 862   the accuracy is better than 10^{-5}.
 863  *                                                                      *
 864  ************************************************************************/
 865 #if 0
 866 double besselJ_meissel(double nu, double x)
 867 {
 868   double beta = pow(nu, 0.325);
 869   double result;
 870
 871   /* Fitted matching points.   */
 872   double llimit = 1.1 * beta;
 873   double ulimit = 1.3 * beta;
 874
 875   double nu2 = nu * nu;
 876
 877   if (nu < 5. && x < 1.)
 878     {
 879       /* Small argument and order. Use a Taylor expansion. */
 880       int k;
 881       double xo2 = 0.5 * x;
 882       double gamfactor = pow(nu,nu) * exp(-nu) * sqrt(nu * 2. * M_PI)
 883         * (1. + 1./(12.*nu) + 1./(288.*nu*nu));
 884       double prefactor = pow(xo2, nu) / gamfactor;
 885       double C[5];
 886
 887       C[0] = 1.;
 888       C[1] = -C[0] / (nu+1.);
 889       C[2] = -C[1] / (2.*(nu+2.));
 890       C[3] = -C[2] / (3.*(nu+3.));
 891       C[4] = -C[3] / (4.*(nu+4.));
 892
 893       result = 0.;
 894       for(k=0; k<5; k++)
 895         result += C[k] * pow(xo2, 2.*k);
 896
 897       result *= prefactor;
 898     }
 899   else if(x < nu - llimit)
 900     {
 901       /* Small x region: x << l.    */
 902       double z = x / nu;
 903       double z2 = z*z;
 904       double rtomz2 = sqrt(1.-z2);
 905       double omz2_2 = (1.-z2)*(1.-z2);
 906
 907       /* Calculate Meissel exponent. */
 908       double term1 = 1./(24.*nu) * ((2.+3.*z2)/((1.-z2)*rtomz2) -2.);
 909       double term2 = - z2*(4. + z2)/(16.*nu2*(1.-z2)*omz2_2);
 910       double V_nu = term1 + term2;
 911
 912       /* Calculate the harmless prefactor. */
 913       double sterlingsum = 1. + 1./(12.*nu) + 1./(288*nu2);
 914       double harmless = 1. / (sqrt(rtomz2*2.*M_PI*nu) * sterlingsum);
 915
 916       /* Calculate the logarithm of the nu dependent prefactor. */
 917       double ln_nupre = rtomz2 + log(z) - log(1. + rtomz2);
 918
 919       result = harmless * exp(nu*ln_nupre - V_nu);
 920     }
 921   else if(x < nu + ulimit)
 922     {
 923       /* Intermediate region 1: x near nu. */
 924       double eps = 1.-nu/x;
 925       double eps_x = eps * x;
 926       double eps_x_2 = eps_x * eps_x;
 927       double xo6 = x/6.;
 928       double B[6];
 929       static double gam[6] = {2.67894, 1.35412, 1., 0.89298, 0.902745, 1.};
 930       static double sf[6] = {0.866025, 0.866025, 0., -0.866025, -0.866025, 0.};
 931
 932       /* Some terms are identically zero, because sf[] can be zero.
 933        * Some terms do not appear in the result.
 934        */
 935       B[0] = 1.;
 936       B[1] = eps_x;
 937       /* B[2] = 0.5 * eps_x_2 - 1./20.; */
 938       B[3] = eps_x * (eps_x_2/6. - 1./15.);
 939       B[4] = eps_x_2 * (eps_x_2 - 1.)/24. + 1./280.;
 940       /* B[5] = eps_x * (eps_x_2*(0.5*eps_x_2 - 1.)/60. + 43./8400.); */
 941
 942       result  = B[0] * gam[0] * sf[0] / pow(xo6, 1./3.);
 943       result += B[1] * gam[1] * sf[1] / pow(xo6, 2./3.);
 944       result += B[3] * gam[3] * sf[3] / pow(xo6, 4./3.);
 945       result += B[4] * gam[4] * sf[4] / pow(xo6, 5./3.);
 946
 947       result /= (3.*M_PI);
 948     }
 949   else
 950     {
 951       /* Region of very large argument. Use expansion
 952        * for x>>l, and we need not be very exacting.
 953        */
 954       double secb = x/nu;
 955       double sec2b= secb*secb;
 956
 957       double cotb = 1./sqrt(sec2b-1.);      /* cotb=cot(beta) */
 958
 959       double beta = acos(nu/x);
 960       double trigarg = nu/cotb - nu*beta - 0.25 * M_PI;
 961
 962       double cot3b = cotb * cotb * cotb;
 963       double cot6b = cot3b * cot3b;
 964
 965       double sum1, sum2, expterm, prefactor, trigcos;
 966
 967       sum1  = 2.0 + 3.0 * sec2b;
 968       trigarg -= sum1 * cot3b / (24.0 * nu);
 969
 970       trigcos = cos(trigarg);
 971
 972       sum2 = 4.0 + sec2b;
 973       expterm = sum2 * sec2b * cot6b / (16.0 * nu2);
 974
 975       expterm = exp(-expterm);
 976       prefactor = sqrt(2. * cotb / (nu * M_PI));
 977
 978       result = prefactor * expterm * trigcos;
 979     }
 980
 981   return  result;
 982 }
 983 #endif