/* specfunc/bessel.c

 * 

 * Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 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 */

/* Miscellaneous support functions for Bessel function evaluations.

 */

#include <config.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_airy.h>

#include <gsl/gsl_sf_elementary.h>

#include <gsl/gsl_sf_exp.h>

#include <gsl/gsl_sf_gamma.h>

#include <gsl/gsl_sf_trig.h>

#include "error.h"

#include "bessel_amp_phase.h"

#include "bessel_temme.h"

#include "bessel.h"

#define CubeRoot2_  1.25992104989487316476721060728

/* Debye functions [Abramowitz+Stegun, 9.3.9-10] */

inline static double 

debye_u1(const double * tpow)

{

  return (3.0*tpow[1] - 5.0*tpow[3])/24.0;

}

inline static double 

debye_u2(const double * tpow)

{

  return (81.0*tpow[2] - 462.0*tpow[4] + 385.0*tpow[6])/1152.0;

}

inline

static double debye_u3(const double * tpow)

{

  return (30375.0*tpow[3] - 369603.0*tpow[5] + 765765.0*tpow[7] - 425425.0*tpow[9])/414720.0;

}

inline

static double debye_u4(const double * tpow)

{

  return (4465125.0*tpow[4] - 94121676.0*tpow[6] + 349922430.0*tpow[8] - 

          446185740.0*tpow[10] + 185910725.0*tpow[12])/39813120.0;

}

inline

static double debye_u5(const double * tpow)

{

  return (1519035525.0*tpow[5]     - 49286948607.0*tpow[7] + 

          284499769554.0*tpow[9]   - 614135872350.0*tpow[11] + 

          566098157625.0*tpow[13]  - 188699385875.0*tpow[15])/6688604160.0;

}

#if 0

inline

static double debye_u6(const double * tpow)

{

  return (2757049477875.0*tpow[6] - 127577298354750.0*tpow[8] + 

          1050760774457901.0*tpow[10] - 3369032068261860.0*tpow[12] + 

          5104696716244125.0*tpow[14] - 3685299006138750.0*tpow[16] + 

          1023694168371875.0*tpow[18])/4815794995200.0;

}

#endif

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

int

gsl_sf_bessel_IJ_taylor_e(const double nu, const double x,

                             const int sign,

                             const int kmax,

                             const double threshold,

                             gsl_sf_result * result

                             )

{

  /* CHECK_POINTER(result) */

  if(nu < 0.0 || x < 0.0) {

    DOMAIN_ERROR(result);

  }

  else if(x == 0.0) {

    if(nu == 0.0) {

      result->val = 1.0;

      result->err = 0.0;

    }

    else {

      result->val = 0.0;

      result->err = 0.0;

    }

    return GSL_SUCCESS;

  }

  else {

    gsl_sf_result prefactor;   /* (x/2)^nu / Gamma(nu+1) */

    gsl_sf_result sum;

    int stat_pre;

    int stat_sum;

    int stat_mul;

    if(nu == 0.0) {

      prefactor.val = 1.0;

      prefactor.err = 0.0;

      stat_pre = GSL_SUCCESS;

    }

    else if(nu < INT_MAX-1) {

      /* Separate the integer part and use

       * y^nu / Gamma(nu+1) = y^N /N! y^f / (N+1)_f,

       * to control the error.

       */

      const int    N = (int)floor(nu + 0.5);

      const double f = nu - N;

      gsl_sf_result poch_factor;

      gsl_sf_result tc_factor;

      const int stat_poch = gsl_sf_poch_e(N+1.0, f, &poch_factor);

      const int stat_tc   = gsl_sf_taylorcoeff_e(N, 0.5*x, &tc_factor);

      const double p = pow(0.5*x,f);

      prefactor.val  = tc_factor.val * p / poch_factor.val;

      prefactor.err  = tc_factor.err * p / poch_factor.val;

      prefactor.err += fabs(prefactor.val) / poch_factor.val * poch_factor.err;

      prefactor.err += 2.0 * GSL_DBL_EPSILON * fabs(prefactor.val);

      stat_pre = GSL_ERROR_SELECT_2(stat_tc, stat_poch);

    }

    else {

      gsl_sf_result lg;

      const int stat_lg = gsl_sf_lngamma_e(nu+1.0, &lg);

      const double term1  = nu*log(0.5*x);

      const double term2  = lg.val;

      const double ln_pre = term1 - term2;

      const double ln_pre_err = GSL_DBL_EPSILON * (fabs(term1)+fabs(term2)) + lg.err;

      const int stat_ex = gsl_sf_exp_err_e(ln_pre, ln_pre_err, &prefactor);

      stat_pre = GSL_ERROR_SELECT_2(stat_ex, stat_lg);

    }

    /* Evaluate the sum.

     * [Abramowitz+Stegun, 9.1.10]

     * [Abramowitz+Stegun, 9.6.7]

     */

    {

      const double y = sign * 0.25 * x*x;

      double sumk = 1.0;

      double term = 1.0;

      int k;

      for(k=1; k<=kmax; k++) {

        term *= y/((nu+k)*k);

        sumk += term;

        if(fabs(term/sumk) < threshold) break;

      }

      sum.val = sumk;

      sum.err = threshold * fabs(sumk);

      stat_sum = ( k >= kmax ? GSL_EMAXITER : GSL_SUCCESS );

    }

    stat_mul = gsl_sf_multiply_err_e(prefactor.val, prefactor.err,

                                        sum.val, sum.err,

                                        result);

    return GSL_ERROR_SELECT_3(stat_mul, stat_pre, stat_sum);

  }

}

/* Hankel's Asymptotic Expansion - A&S 9.2.5

 *  

 * x >> nu*nu+1

 * error ~ O( ((nu*nu+1)/x)^4 )

 *

 * empirical error analysis:

 *   choose  GSL_ROOT4_MACH_EPS * x > (nu*nu + 1)

 *

 * This is not especially useful. When the argument gets

 * large enough for this to apply, the cos() and sin()

 * start loosing digits. However, this seems inevitable

 * for this particular method.

 *

 * Wed Jun 25 14:39:38 MDT 2003 [GJ]

 * This function was inconsistent since the Q term did not

 * go to relative order eps^2. That's why the error estimate

 * originally given was screwy (it didn't make sense that the

 * "empirical" error was coming out O(eps^3)).

 * With Q to proper order, the error is O(eps^4).

 *

 * Sat Mar 15 05:16:18 GMT 2008 [BG]

 * Extended to use additional terms in the series to gain

 * higher accuracy.

 *

 */

int

gsl_sf_bessel_Jnu_asympx_e(const double nu, const double x, gsl_sf_result * result)

{

  double mu  = 4.0*nu*nu;

  double chi = x - (0.5*nu + 0.25)*M_PI;

  double P   = 0.0; 

  double Q   = 0.0;

  double k = 0, t = 1;

  int convP, convQ;

  do

    {

      t *= (k == 0) ? 1 : -(mu - (2*k-1)*(2*k-1)) / (k * (8 * x));

      convP = fabs(t) < GSL_DBL_EPSILON * fabs(P);

      P += t;

      k++;

      t *= (mu - (2*k-1)*(2*k-1)) / (k * (8 * x));

      convQ = fabs(t) < GSL_DBL_EPSILON * fabs(Q);

      Q += t;

      /* To preserve the consistency of the series we need to exit

         when P and Q have the same number of terms */

      if (convP && convQ && k > (nu / 2)) 

        break;

      k++;

    }

  while (k < 1000);

  {

    double pre = sqrt(2.0/(M_PI*x));

    double c   = cos(chi);

    double s   = sin(chi);

    result->val  = pre * (c*P - s*Q);

    result->err  = pre * GSL_DBL_EPSILON * (fabs(c*P) + fabs(s*Q) + fabs(t)) * (1 + fabs(x));

    /* NB: final term accounts for phase error with large x */

  }

  return GSL_SUCCESS;

}

/* x >> nu*nu+1

 */

int

gsl_sf_bessel_Ynu_asympx_e(const double nu, const double x, gsl_sf_result * result)

{

  double ampl;

  double theta;

  double alpha = x;

  double beta  = -0.5*nu*M_PI;

  int stat_a = gsl_sf_bessel_asymp_Mnu_e(nu, x, &l;;

  int stat_t = gsl_sf_bessel_asymp_thetanu_corr_e(nu, x, &theta);

  double sin_alpha = sin(alpha);

  double cos_alpha = cos(alpha);

  double sin_chi   = sin(beta + theta);

  double cos_chi   = cos(beta + theta);

  double sin_term     = sin_alpha * cos_chi + sin_chi * cos_alpha;

  double sin_term_mag = fabs(sin_alpha * cos_chi) + fabs(sin_chi * cos_alpha);

  result->val  = ampl * sin_term;

  result->err  = fabs(ampl) * GSL_DBL_EPSILON * sin_term_mag;

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

  if(fabs(alpha) > 1.0/GSL_DBL_EPSILON) {

    result->err *= 0.5 * fabs(alpha);

  }

  else if(fabs(alpha) > 1.0/GSL_SQRT_DBL_EPSILON) {

    result->err *= 256.0 * fabs(alpha) * GSL_SQRT_DBL_EPSILON;

  }

  return GSL_ERROR_SELECT_2(stat_t, stat_a);

}

/* x >> nu*nu+1

 */

int

gsl_sf_bessel_Inu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)

{

  double mu   = 4.0*nu*nu;

  double mum1 = mu-1.0;

  double mum9 = mu-9.0;

  double pre  = 1.0/sqrt(2.0*M_PI*x);

  double r    = mu/x;

  result->val = pre * (1.0 - mum1/(8.0*x) + mum1*mum9/(128.0*x*x));

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

  return GSL_SUCCESS;

}

/* x >> nu*nu+1

 */

int

gsl_sf_bessel_Knu_scaled_asympx_e(const double nu, const double x, gsl_sf_result * result)

{

  double mu   = 4.0*nu*nu;

  double mum1 = mu-1.0;

  double mum9 = mu-9.0;

  double pre  = sqrt(M_PI/(2.0*x));

  double r    = nu/x;

  result->val = pre * (1.0 + mum1/(8.0*x) + mum1*mum9/(128.0*x*x));

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

  return GSL_SUCCESS;

}

/* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.7]

 *

 * error:

 *   The error has the form u_N(t)/nu^N  where  0 <= t <= 1.

 *   It is not hard to show that |u_N(t)| is small for such t.

 *   We have N=6 here, and |u_6(t)| < 0.025, so the error is clearly

 *   bounded by 0.025/nu^6. This gives the asymptotic bound on nu

 *   seen below as nu ~ 100. For general MACH_EPS it will be 

 *                     nu > 0.5 / MACH_EPS^(1/6)

 *   When t is small, the bound is even better because |u_N(t)| vanishes

 *   as t->0. In fact u_N(t) ~ C t^N as t->0, with C ~= 0.1.

 *   We write

 *                     err_N <= min(0.025, C(1/(1+(x/nu)^2))^3) / nu^6

 *   therefore

 *                     min(0.29/nu^2, 0.5/(nu^2+x^2)) < MACH_EPS^{1/3}

 *   and this is the general form.

 *

 * empirical error analysis, assuming 14 digit requirement:

 *   choose   x > 50.000 nu   ==>  nu >   3

 *   choose   x > 10.000 nu   ==>  nu >  15

 *   choose   x >  2.000 nu   ==>  nu >  50

 *   choose   x >  1.000 nu   ==>  nu >  75

 *   choose   x >  0.500 nu   ==>  nu >  80

 *   choose   x >  0.100 nu   ==>  nu >  83

 *

 * This makes sense. For x << nu, the error will be of the form u_N(1)/nu^N,

 * since the polynomial term will be evaluated near t=1, so the bound

 * on nu will become constant for small x. Furthermore, increasing x with

 * nu fixed will decrease the error.

 */

int

gsl_sf_bessel_Inu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)

{

  int i;

  double z = x/nu;

  double root_term = hypot(1.0,z);

  double pre = 1.0/sqrt(2.0*M_PI*nu * root_term);

  double eta = root_term + log(z/(1.0+root_term));

  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)) );

  gsl_sf_result ex_result;

  int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);

  if(stat_ex == GSL_SUCCESS) {

    double t = 1.0/root_term;

    double sum;

    double tpow[16];

    tpow[0] = 1.0;

    for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];

    sum = 1.0 + debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) + debye_u3(tpow)/(nu*nu*nu)

          + debye_u4(tpow)/(nu*nu*nu*nu) + debye_u5(tpow)/(nu*nu*nu*nu*nu);

    result->val  = pre * ex_result.val * sum;

    result->err  = pre * ex_result.val / (nu*nu*nu*nu*nu*nu);

    result->err += pre * ex_result.err * fabs(sum);

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

    return GSL_SUCCESS;

  }

  else {

    result->val = 0.0;

    result->err = 0.0;

    return stat_ex;

  }

}

/* nu -> Inf; uniform in x > 0  [Abramowitz+Stegun, 9.7.8]

 *

 * error:

 *   identical to that above for Inu_scaled

 */

int

gsl_sf_bessel_Knu_scaled_asymp_unif_e(const double nu, const double x, gsl_sf_result * result)

{

  int i;

  double z = x/nu;

  double root_term = hypot(1.0,z);

  double pre = sqrt(M_PI/(2.0*nu*root_term));

  double eta = root_term + log(z/(1.0+root_term));

  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)) );

  gsl_sf_result ex_result;

  int stat_ex = gsl_sf_exp_e(ex_arg, &ex_result);

  if(stat_ex == GSL_SUCCESS) {

    double t = 1.0/root_term;

    double sum;

    double tpow[16];

    tpow[0] = 1.0;

    for(i=1; i<16; i++) tpow[i] = t * tpow[i-1];

    sum = 1.0 - debye_u1(tpow)/nu + debye_u2(tpow)/(nu*nu) - debye_u3(tpow)/(nu*nu*nu)

          + debye_u4(tpow)/(nu*nu*nu*nu) - debye_u5(tpow)/(nu*nu*nu*nu*nu);

    result->val  = pre * ex_result.val * sum;

    result->err  = pre * ex_result.err * fabs(sum);

    result->err += pre * ex_result.val / (nu*nu*nu*nu*nu*nu);

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

    return GSL_SUCCESS;

  }

  else {

    result->val = 0.0;

    result->err = 0.0;

    return stat_ex;

  }

}

/* Evaluate J_mu(x),J_{mu+1}(x) and Y_mu(x),Y_{mu+1}(x)  for |mu| < 1/2

 */

int

gsl_sf_bessel_JY_mu_restricted(const double mu, const double x,

                               gsl_sf_result * Jmu, gsl_sf_result * Jmup1,

                               gsl_sf_result * Ymu, gsl_sf_result * Ymup1)

{

  /* CHECK_POINTER(Jmu) */

  /* CHECK_POINTER(Jmup1) */

  /* CHECK_POINTER(Ymu) */

  /* CHECK_POINTER(Ymup1) */

  if(x < 0.0 || fabs(mu) > 0.5) {

    Jmu->val   = 0.0;

    Jmu->err   = 0.0;

    Jmup1->val = 0.0;

    Jmup1->err = 0.0;

    Ymu->val   = 0.0;

    Ymu->err   = 0.0;

    Ymup1->val = 0.0;

    Ymup1->err = 0.0;

    GSL_ERROR ("error", GSL_EDOM);

  }

  else if(x == 0.0) {

    if(mu == 0.0) {

      Jmu->val   = 1.0;

      Jmu->err   = 0.0;

    }

    else {

      Jmu->val   = 0.0;

      Jmu->err   = 0.0;

    }

    Jmup1->val = 0.0;

    Jmup1->err = 0.0;

    Ymu->val   = 0.0;

    Ymu->err   = 0.0;

    Ymup1->val = 0.0;

    Ymup1->err = 0.0;

    GSL_ERROR ("error", GSL_EDOM);

  }

  else {

    int stat_Y;

    int stat_J;

    if(x < 2.0) {

      /* Use Taylor series for J and the Temme series for Y.

       * The Taylor series for J requires nu > 0, so we shift

       * up one and use the recursion relation to get Jmu, in

       * case mu < 0.

       */

      gsl_sf_result Jmup2;

      int stat_J1 = gsl_sf_bessel_IJ_taylor_e(mu+1.0, x, -1, 100, GSL_DBL_EPSILON,  Jmup1);

      int stat_J2 = gsl_sf_bessel_IJ_taylor_e(mu+2.0, x, -1, 100, GSL_DBL_EPSILON, &Jmup2);

      double c = 2.0*(mu+1.0)/x;

      Jmu->val  = c * Jmup1->val - Jmup2.val;

      Jmu->err  = c * Jmup1->err + Jmup2.err;

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

      stat_J = GSL_ERROR_SELECT_2(stat_J1, stat_J2);

      stat_Y = gsl_sf_bessel_Y_temme(mu, x, Ymu, Ymup1);

      return GSL_ERROR_SELECT_2(stat_J, stat_Y);

    }

    else if(x < 1000.0) {

      double P, Q;

      double J_ratio;

      double J_sgn;

      const int stat_CF1 = gsl_sf_bessel_J_CF1(mu, x, &J_ratio, &J_sgn);

      const int stat_CF2 = gsl_sf_bessel_JY_steed_CF2(mu, x, &P, &Q);

      double Jprime_J_ratio = mu/x - J_ratio;

      double gamma = (P - Jprime_J_ratio)/Q;

      Jmu->val = J_sgn * sqrt(2.0/(M_PI*x) / (Q + gamma*(P-Jprime_J_ratio)));

      Jmu->err = 4.0 * GSL_DBL_EPSILON * fabs(Jmu->val);

      Jmup1->val = J_ratio * Jmu->val;

      Jmup1->err = fabs(J_ratio) * Jmu->err;

      Ymu->val = gamma * Jmu->val;

      Ymu->err = fabs(gamma) * Jmu->err;

      Ymup1->val = Ymu->val * (mu/x - P - Q/gamma);

      Ymup1->err = Ymu->err * fabs(mu/x - P - Q/gamma) + 4.0*GSL_DBL_EPSILON*fabs(Ymup1->val);

      return GSL_ERROR_SELECT_2(stat_CF1, stat_CF2);

    }

    else {

      /* Use asymptotics for large argument.

       */

      const int stat_J0 = gsl_sf_bessel_Jnu_asympx_e(mu,     x, Jmu);

      const int stat_J1 = gsl_sf_bessel_Jnu_asympx_e(mu+1.0, x, Jmup1);

      const int stat_Y0 = gsl_sf_bessel_Ynu_asympx_e(mu,     x, Ymu);

      const int stat_Y1 = gsl_sf_bessel_Ynu_asympx_e(mu+1.0, x, Ymup1);

      stat_J = GSL_ERROR_SELECT_2(stat_J0, stat_J1);

      stat_Y = GSL_ERROR_SELECT_2(stat_Y0, stat_Y1);

      return GSL_ERROR_SELECT_2(stat_J, stat_Y);

    }

  }

}

int

gsl_sf_bessel_J_CF1(const double nu, const double x,

                    double * ratio, double * sgn)

{

  const double RECUR_BIG = GSL_SQRT_DBL_MAX;

  const double RECUR_SMALL = GSL_SQRT_DBL_MIN;

  const int maxiter = 10000;

  int n = 1;

  double Anm2 = 1.0;

  double Bnm2 = 0.0;

  double Anm1 = 0.0;

  double Bnm1 = 1.0;

  double a1 = x/(2.0*(nu+1.0));

  double An = Anm1 + a1*Anm2;

  double Bn = Bnm1 + a1*Bnm2;

  double an;

  double fn = An/Bn;

  double dn = a1;

  double s  = 1.0;

  while(n < maxiter) {

    double old_fn;

    double del;

    n++;

    Anm2 = Anm1;

    Bnm2 = Bnm1;

    Anm1 = An;

    Bnm1 = Bn;

    an = -x*x/(4.0*(nu+n-1.0)*(nu+n));

    An = Anm1 + an*Anm2;

    Bn = Bnm1 + an*Bnm2;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {

      An /= RECUR_BIG;

      Bn /= RECUR_BIG;

      Anm1 /= RECUR_BIG;

      Bnm1 /= RECUR_BIG;

      Anm2 /= RECUR_BIG;

    } else if(fabs(An) < RECUR_SMALL || fabs(Bn) < RECUR_SMALL) {

      An /= RECUR_SMALL;

      Bn /= RECUR_SMALL;

      Anm1 /= RECUR_SMALL;

      Bnm1 /= RECUR_SMALL;

      Anm2 /= RECUR_SMALL;

      Bnm2 /= RECUR_SMALL;

    }

    old_fn = fn;

    fn = An/Bn;

    del = old_fn/fn;

    dn = 1.0 / (2.0*(nu+n)/x - dn);

    if(dn < 0.0) s = -s;

    if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;

  }

  /* FIXME: we should return an error term here as well, because the

     error from this recurrence affects the overall error estimate. */

  *ratio = fn;

  *sgn   = s;

  if(n >= maxiter)

    GSL_ERROR ("error", GSL_EMAXITER);

  else

    return GSL_SUCCESS;

}

/* Evaluate the continued fraction CF1 for J_{nu+1}/J_nu

 * using Gautschi (Euler) equivalent series.

 * This exhibits an annoying problem because the

 * a_k are not positive definite (in fact they are all negative).

 * There are cases when rho_k blows up. Example: nu=1,x=4.

 */

#if 0

int

gsl_sf_bessel_J_CF1_ser(const double nu, const double x,

                        double * ratio, double * sgn)

{

  const int maxk = 20000;

  double tk   = 1.0;

  double sum  = 1.0;

  double rhok = 0.0;

  double dk = 0.0;

  double s  = 1.0;

  int k;

  for(k=1; k

    double ak = -0.25 * (x/(nu+k)) * x/(nu+k+1.0);

    rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));

    tk  *= rhok;

    sum += tk;

    dk = 1.0 / (2.0/x - (nu+k-1.0)/(nu+k) * dk);

    if(dk < 0.0) s = -s;

    if(fabs(tk/sum) < GSL_DBL_EPSILON) break;

  }

  *ratio = x/(2.0*(nu+1.0)) * sum;

  *sgn   = s;

  if(k == maxk)

    GSL_ERROR ("error", GSL_EMAXITER);

  else

    return GSL_SUCCESS;

}

#endif

/* Evaluate the continued fraction CF1 for I_{nu+1}/I_nu

 * using Gautschi (Euler) equivalent series.

 */

int

gsl_sf_bessel_I_CF1_ser(const double nu, const double x, double * ratio)

{

  const int maxk = 20000;

  double tk   = 1.0;

  double sum  = 1.0;

  double rhok = 0.0;

  int k;

  for(k=1; k

    double ak = 0.25 * (x/(nu+k)) * x/(nu+k+1.0);

    rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));

    tk  *= rhok;

    sum += tk;

    if(fabs(tk/sum) < GSL_DBL_EPSILON) break;

  }

  *ratio = x/(2.0*(nu+1.0)) * sum;

  if(k == maxk)

    GSL_ERROR ("error", GSL_EMAXITER);

  else

    return GSL_SUCCESS;

}

int

gsl_sf_bessel_JY_steed_CF2(const double nu, const double x,

                           double * P, double * Q)

{

  const int max_iter = 10000;

  const double SMALL = 1.0e-100;

  int i = 1;

  double x_inv = 1.0/x;

  double a = 0.25 - nu*nu;

  double p = -0.5*x_inv;

  double q = 1.0;

  double br = 2.0*x;

  double bi = 2.0;

  double fact = a*x_inv/(p*p + q*q);

  double cr = br + q*fact;

  double ci = bi + p*fact;

  double den = br*br + bi*bi;

  double dr = br/den;

  double di = -bi/den;

  double dlr = cr*dr - ci*di;

  double dli = cr*di + ci*dr;

  double temp = p*dlr - q*dli;

  q = p*dli + q*dlr;

  p = temp;

  for (i=2; i<=max_iter; i++) {

    a  += 2*(i-1);

    bi += 2.0;

    dr = a*dr + br;

    di = a*di + bi;

    if(fabs(dr)+fabs(di) < SMALL) dr = SMALL;

    fact = a/(cr*cr+ci*ci);

    cr = br + cr*fact;

    ci = bi - ci*fact;

    if(fabs(cr)+fabs(ci) < SMALL) cr = SMALL;

    den = dr*dr + di*di;

    dr /= den;

    di /= -den;

    dlr = cr*dr - ci*di;

    dli = cr*di + ci*dr;

    temp = p*dlr - q*dli;

    q = p*dli + q*dlr;

    p = temp;

    if(fabs(dlr-1.0)+fabs(dli) < GSL_DBL_EPSILON) break;

  }

  *P = p;

  *Q = q;

  if(i == max_iter)

    GSL_ERROR ("error", GSL_EMAXITER);

  else

    return GSL_SUCCESS;

}

/* Evaluate continued fraction CF2, using Thompson-Barnett-Temme method,

 * to obtain values of exp(x)*K_nu and exp(x)*K_{nu+1}.

 *

 * This is unstable for small x; x > 2 is a good cutoff.

 * Also requires |nu| < 1/2.

 */

int

gsl_sf_bessel_K_scaled_steed_temme_CF2(const double nu, const double x,

                                       double * K_nu, double * K_nup1,

                                       double * Kp_nu)

{

  const int maxiter = 10000;

  int i = 1;

  double bi = 2.0*(1.0 + x);

  double di = 1.0/bi;

  double delhi = di;

  double hi    = di;

  double qi   = 0.0;

  double qip1 = 1.0;

  double ai = -(0.25 - nu*nu);

  double a1 = ai;

  double ci = -ai;

  double Qi = -ai;

  double s = 1.0 + Qi*delhi;

  for(i=2; i<=maxiter; i++) {

    double dels;

    double tmp;

    ai -= 2.0*(i-1);

    ci  = -ai*ci/i;

    tmp  = (qi - bi*qip1)/ai;

    qi   = qip1;

    qip1 = tmp;

    Qi += ci*qip1;

    bi += 2.0;

    di  = 1.0/(bi + ai*di);

    delhi = (bi*di - 1.0) * delhi;

    hi += delhi;

    dels = Qi*delhi;

    s += dels;

    if(fabs(dels/s) < GSL_DBL_EPSILON) break;

  }

  hi *= -a1;

  *K_nu   = sqrt(M_PI/(2.0*x)) / s;

  *K_nup1 = *K_nu * (nu + x + 0.5 - hi)/x;

  *Kp_nu  = - *K_nup1 + nu/x * *K_nu;

  if(i == maxiter)

    GSL_ERROR ("error", GSL_EMAXITER);

  else

    return GSL_SUCCESS;

}

int gsl_sf_bessel_cos_pi4_e(double y, double eps, gsl_sf_result * result)

{

  const double sy = sin(y);

  const double cy = cos(y);

  const double s = sy + cy;

  const double d = sy - cy;

  const double abs_sum = fabs(cy) + fabs(sy);

  double seps;

  double ceps;

  if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {

    const double e2 = eps*eps;

    seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));

    ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);

  }

  else {

    seps = sin(eps);

    ceps = cos(eps);

  }

  result->val = (ceps * s - seps * d)/ M_SQRT2;

  result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;

  /* Try to account for error in evaluation of sin(y), cos(y).

   * This is a little sticky because we don't really know

   * how the library routines are doing their argument reduction.

   * However, we will make a reasonable guess.

   * FIXME ?

   */

  if(y > 1.0/GSL_DBL_EPSILON) {

    result->err *= 0.5 * y;

  }

  else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {

    result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;

  }

  return GSL_SUCCESS;

}

int gsl_sf_bessel_sin_pi4_e(double y, double eps, gsl_sf_result * result)

{

  const double sy = sin(y);

  const double cy = cos(y);

  const double s = sy + cy;

  const double d = sy - cy;

  const double abs_sum = fabs(cy) + fabs(sy);

  double seps;

  double ceps;

  if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {

    const double e2 = eps*eps;

    seps = eps * (1.0 - e2/6.0 * (1.0 - e2/20.0));

    ceps = 1.0 - e2/2.0 * (1.0 - e2/12.0);

  }

  else {

    seps = sin(eps);

    ceps = cos(eps);

  }

  result->val = (ceps * d + seps * s)/ M_SQRT2;

  result->err = 2.0 * GSL_DBL_EPSILON * (fabs(ceps) + fabs(seps)) * abs_sum / M_SQRT2;

  /* Try to account for error in evaluation of sin(y), cos(y).

   * See above.

   * FIXME ?

   */

  if(y > 1.0/GSL_DBL_EPSILON) {

    result->err *= 0.5 * y;

  }

  else if(y > 1.0/GSL_SQRT_DBL_EPSILON) {

    result->err *= 256.0 * y * GSL_SQRT_DBL_EPSILON;

  }

  return GSL_SUCCESS;

}

/************************************************************************

 *                                                                      *

  Asymptotic approximations 8.11.5, 8.12.5, and 8.42.7 from

  G.N.Watson, A Treatise on the Theory of Bessel Functions,

  2nd Edition (Cambridge University Press, 1944).

  Higher terms in expansion for x near l given by

  Airey in Phil. Mag. 31, 520 (1916).

  This approximation is accurate to near 0.1% at the boundaries

  between the asymptotic regions; well away from the boundaries

  the accuracy is better than 10^{-5}.

 *                                                                      *

 ************************************************************************/

#if 0

double besselJ_meissel(double nu, double x)

{

  double beta = pow(nu, 0.325);

  double result;

  /* Fitted matching points.   */

  double llimit = 1.1 * beta;

  double ulimit = 1.3 * beta;

  double nu2 = nu * nu;

  if (nu < 5. && x < 1.)

    {

      /* Small argument and order. Use a Taylor expansion. */

      int k;