自学教程：C++ GSL_ERROR_SELECT_2函数代码示例

intgsl_sf_legendre_H3d_array(const int lmax, const double lambda, const double eta, double * result_array){  /* CHECK_POINTER(result_array) */ if(eta < 0.0 || lmax < 0) {    int ell;    for(ell=0; ell<=lmax; ell++) result_array[ell] = 0.0;    GSL_ERROR ("domain error", GSL_EDOM);  }  else if(eta > GSL_LOG_DBL_MAX) {    /* cosh(eta) is too big. */    int ell;    for(ell=0; ell<=lmax; ell++) result_array[ell] = 0.0;    GSL_ERROR ("overflow", GSL_EOVRFLW);  }  else if(lmax == 0) {    gsl_sf_result H0;    int stat = gsl_sf_legendre_H3d_e(0, lambda, eta, &H0);    result_array[0] = H0.val;    return stat;  }  else {    /* Not the most efficient method. But what the hell... it's simple.     */    gsl_sf_result r_Hlp1;    gsl_sf_result r_Hl;    int stat_lmax   = gsl_sf_legendre_H3d_e(lmax,   lambda, eta, &r_Hlp1);    int stat_lmaxm1 = gsl_sf_legendre_H3d_e(lmax-1, lambda, eta, &r_Hl);    int stat_max = GSL_ERROR_SELECT_2(stat_lmax, stat_lmaxm1);    const double coth_eta = 1.0/tanh(eta);    int stat_recursion = GSL_SUCCESS;    double Hlp1 = r_Hlp1.val;    double Hl   = r_Hl.val;    double Hlm1;    int ell;    result_array[lmax]   = Hlp1;    result_array[lmax-1] = Hl;    for(ell=lmax-1; ell>0; ell--) {      double root_term_0 = sqrt(lambda*lambda + (double)ell*ell);      double root_term_1 = sqrt(lambda*lambda + (ell+1.0)*(ell+1.0));      Hlm1 = ((2.0*ell + 1.0)*coth_eta*Hl - root_term_1 * Hlp1)/root_term_0;      result_array[ell-1] = Hlm1;      if(!(Hlm1 < GSL_DBL_MAX)) stat_recursion = GSL_EOVRFLW;      Hlp1 = Hl;      Hl   = Hlm1;    }    return GSL_ERROR_SELECT_2(stat_recursion, stat_max);  }}

intgsl_sf_hyperg_2F1_renorm_e(const double a, const double b, const double c,                              const double x,                              gsl_sf_result * result                              ){  const double rinta = floor(a + 0.5);  const double rintb = floor(b + 0.5);  const double rintc = floor(c + 0.5);  const int a_neg_integer = ( a < 0.0  &&  fabs(a - rinta) < locEPS );  const int b_neg_integer = ( b < 0.0  &&  fabs(b - rintb) < locEPS );  const int c_neg_integer = ( c < 0.0  &&  fabs(c - rintc) < locEPS );    if(c_neg_integer) {    if((a_neg_integer && a > c+0.1) || (b_neg_integer && b > c+0.1)) {      /* 2F1 terminates early */      result->val = 0.0;      result->err = 0.0;      return GSL_SUCCESS;    }    else {      /* 2F1 does not terminate early enough, so something survives */      /* [Abramowitz+Stegun, 15.1.2] */      gsl_sf_result g1, g2, g3, g4, g5;      double s1, s2, s3, s4, s5;      int stat = 0;      stat += gsl_sf_lngamma_sgn_e(a-c+1, &g1, &s1);      stat += gsl_sf_lngamma_sgn_e(b-c+1, &g2, &s2);      stat += gsl_sf_lngamma_sgn_e(a, &g3, &s3);      stat += gsl_sf_lngamma_sgn_e(b, &g4, &s4);      stat += gsl_sf_lngamma_sgn_e(-c+2, &g5, &s5);      if(stat != 0) {        DOMAIN_ERROR(result);      }      else {        gsl_sf_result F;        int stat_F = gsl_sf_hyperg_2F1_e(a-c+1, b-c+1, -c+2, x, &F);        double ln_pre_val = g1.val + g2.val - g3.val - g4.val - g5.val;        double ln_pre_err = g1.err + g2.err + g3.err + g4.err + g5.err;        double sg  = s1 * s2 * s3 * s4 * s5;        int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,                                              sg * F.val, F.err,                                              result);        return GSL_ERROR_SELECT_2(stat_e, stat_F);      }    }  }  else {    /* generic c */    gsl_sf_result F;    gsl_sf_result lng;    double sgn;    int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);    int stat_F = gsl_sf_hyperg_2F1_e(a, b, c, x, &F);    int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,                                          sgn*F.val, F.err,                                          result);    return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);  }}

staticintbessel_J_recur_asymp(const double nu, const double x,                     gsl_sf_result * Jnu, gsl_sf_result * Jnup1){  const double nu_cut = 25.0;  int n;  int steps = ceil(nu_cut - nu) + 1;  gsl_sf_result r_Jnp1;  gsl_sf_result r_Jn;  int stat_O1 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps + 1.0, x, &r_Jnp1);  int stat_O2 = gsl_sf_bessel_Jnu_asymp_Olver_e(nu + steps,       x, &r_Jn);  double r_fe = fabs(r_Jnp1.err/r_Jnp1.val) + fabs(r_Jn.err/r_Jn.val);  double Jnp1 = r_Jnp1.val;  double Jn   = r_Jn.val;  double Jnm1;  double Jnp1_save;  for(n=steps; n>0; n--) {    Jnm1 = 2.0*(nu+n)/x * Jn - Jnp1;    Jnp1 = Jn;    Jnp1_save = Jn;    Jn   = Jnm1;  }  Jnu->val = Jn;  Jnu->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jn);  Jnup1->val = Jnp1_save;  Jnup1->err = (r_fe + GSL_DBL_EPSILON * (steps + 1.0)) * fabs(Jnp1_save);  return GSL_ERROR_SELECT_2(stat_O1, stat_O2);}

int gsl_sf_bessel_yl_array(const int lmax, const double x, double * result_array){    /* CHECK_POINTER(result_array) */    if(lmax < 0 || x <= 0.0) {        GSL_ERROR ("error", GSL_EDOM);    } else if (lmax == 0) {        gsl_sf_result result;        int stat = gsl_sf_bessel_y0_e(x, &result);        result_array[0] = result.val;        return stat;    } else {        gsl_sf_result r_yell;        gsl_sf_result r_yellm1;        int stat_1 = gsl_sf_bessel_y1_e(x, &r_yell);        int stat_0 = gsl_sf_bessel_y0_e(x, &r_yellm1);        double yellp1;        double yell   = r_yell.val;        double yellm1 = r_yellm1.val;        int ell;        result_array[0] = yellm1;        result_array[1] = yell;        for(ell = 1; ell < lmax; ell++) {            yellp1 = (2*ell+1)/x * yell - yellm1;            result_array[ell+1] = yellp1;            yellm1 = yell;            yell   = yellp1;        }        return GSL_ERROR_SELECT_2(stat_0, stat_1);    }}

/* Evaluate bessel_I(nu, x), allowing nu < 0. * This is fine here because we do not not allow * nu to be a negative integer. * x > 0. */staticinthyperg_0F1_bessel_I(const double nu, const double x, gsl_sf_result * result){  if(x > GSL_LOG_DBL_MAX) {    OVERFLOW_ERROR(result);  }  if(nu < 0.0) {     const double anu = -nu;    const double s   = 2.0/M_PI * sin(anu*M_PI);    const double ex  = exp(x);    gsl_sf_result I;    gsl_sf_result K;    int stat_I = gsl_sf_bessel_Inu_scaled_e(anu, x, &I);    int stat_K = gsl_sf_bessel_Knu_scaled_e(anu, x, &K);    result->val  = ex * I.val + s * (K.val / ex);    result->err  = ex * I.err + fabs(s * K.err/ex);    result->err += fabs(s * (K.val/ex)) * GSL_DBL_EPSILON * anu * M_PI;    return GSL_ERROR_SELECT_2(stat_K, stat_I);  }  else {    const double ex  = exp(x);    gsl_sf_result I;    int stat_I = gsl_sf_bessel_Inu_scaled_e(nu, x, &I);    result->val = ex * I.val;    result->err = ex * I.err + GSL_DBL_EPSILON * fabs(result->val);    return stat_I;  }}

intgsl_sf_hydrogenicR_e(const int n, const int l,                        const double Z, const double r,                        gsl_sf_result * result){  if(n < 1 || l > n-1 || Z <= 0.0 || r < 0.0) {    DOMAIN_ERROR(result);  }  else {    double A = 2.0*Z/n;    gsl_sf_result norm;    int stat_norm = R_norm(n, l, Z, &norm);    double rho = A*r;    double ea = exp(-0.5*rho);    double pp = gsl_sf_pow_int(rho, l);    gsl_sf_result lag;    int stat_lag = gsl_sf_laguerre_n_e(n-l-1, 2*l+1, rho, &lag);    double W_val = norm.val * ea * pp;    double W_err = norm.err * ea * pp;    W_err += norm.val * ((0.5*rho + 1.0) * GSL_DBL_EPSILON) * ea * pp;    W_err += norm.val * ea * ((l+1.0) * GSL_DBL_EPSILON) * pp;    result->val  = W_val * lag.val;    result->err  = W_val * lag.err + W_err * fabs(lag.val);    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    if ((l == 0 || (r > 0 && l > 0)) && lag.val != 0.0         && stat_lag == GSL_SUCCESS && stat_norm == GSL_SUCCESS) {      CHECK_UNDERFLOW(result);    };    return GSL_ERROR_SELECT_2(stat_lag, stat_norm);  }}

int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(x <= 0.0) {    DOMAIN_ERROR(result);  }  else if(x < 2.0*GSL_DBL_MIN) {    OVERFLOW_ERROR(result);  }  else if(x <= 2.0) {    const double lx = log(x);    int stat_I1;    gsl_sf_result I1;    gsl_sf_result c;    cheb_eval_e(&bk1_cs, 0.5*x*x-1.0, &c);    stat_I1 = gsl_sf_bessel_I1_e(x, &I1);    result->val  = (lx-M_LN2)*I1.val + (0.75 + c.val)/x;    result->err  = c.err/x + fabs(lx)*I1.err;    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return stat_I1;  }  else {    gsl_sf_result K1_scaled;    int stat_K1 = gsl_sf_bessel_K1_scaled_e(x, &K1_scaled);    int stat_e  = gsl_sf_exp_mult_err_e(-x, 0.0,                                           K1_scaled.val, K1_scaled.err,					   result);    result->err = fabs(result->val) * (GSL_DBL_EPSILON*fabs(x) + K1_scaled.err/K1_scaled.val);    return GSL_ERROR_SELECT_2(stat_e, stat_K1);  }}

/* x >> nu*nu+1 */intgsl_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, &ampl);  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);}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.6)] */intgsl_sf_ellint_Ecomp_e(double k, gsl_mode_t mode, gsl_sf_result * result){  if(k*k >= 1.0) {    DOMAIN_ERROR(result);  }  else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {    /* [Abramowitz+Stegun, 17.3.36] */    const double y = 1.0 - k*k;    const double a[] = { 0.44325141463, 0.06260601220, 0.04757383546 };    const double b[] = { 0.24998368310, 0.09200180037, 0.04069697526 };    const double ta = 1.0 + y*(a[0] + y*(a[1] + a[2]*y));    const double tb = -y*log(y) * (b[0] + y*(b[1] + b[2]*y));    result->val = ta + tb;    result->err = 2.0 * GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else {    gsl_sf_result rf;    gsl_sf_result rd;    const double y = 1.0 - k*k;    const int rfstatus = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, &rf);    const int rdstatus = gsl_sf_ellint_RD_e(0.0, y, 1.0, mode, &rd);    result->val = rf.val - k*k/3.0 * rd.val;    result->err = rf.err + k*k/3.0 * rd.err;    return GSL_ERROR_SELECT_2(rfstatus, rdstatus);  }}

int gsl_sf_bessel_K0_e(const double x, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(x <= 0.0) {    DOMAIN_ERROR(result);  }  else if(x <= 2.0) {    const double lx = log(x);    int stat_I0;    gsl_sf_result I0;    gsl_sf_result c;    cheb_eval_e(&bk0_cs, 0.5*x*x-1.0, &c);    stat_I0 = gsl_sf_bessel_I0_e(x, &I0);    result->val  = (-lx+M_LN2)*I0.val - 0.25 + c.val;    result->err  = (fabs(lx) + M_LN2) * I0.err + c.err;    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return stat_I0;  }  else {    gsl_sf_result K0_scaled;    int stat_K0 = gsl_sf_bessel_K0_scaled_e(x, &K0_scaled);    int stat_e  = gsl_sf_exp_mult_err_e(-x, GSL_DBL_EPSILON*fabs(x),                                           K0_scaled.val, K0_scaled.err,					   result);    return GSL_ERROR_SELECT_2(stat_e, stat_K0);  }}

/* See note above for  gamma_inc_Q_CF(). */staticintgamma_inc_Q_CF_protected(const double a, const double x, gsl_sf_result * result){  if(fabs(1.0 - a + x) < 2.0*GSL_DBL_EPSILON) {    /*     * This is a problem region because when     * 1.0 - a + x = 0 the first term of the     * CF recursion is undefined.     *     * I missed this condition when I first     * implemented gamma_inc_Q_CF() function,     * so now I have to fix it by side-stepping     * this point, using the recursion relation     *   Q(a,x) = Q(a-1,x) + x^(a-1) e^(-z) / Gamma(a)     *          = Q(a-1,x) + D(a-1,x)     * to lower 'a' by one, giving an a=x point,     * which is fine.     */    gsl_sf_result D;    gsl_sf_result tmp_CF;    const int stat_tmp_CF = gamma_inc_Q_CF(a-1.0, x, &tmp_CF);    const int stat_D = gamma_inc_D(a-1.0, x, &D);    result->val  = tmp_CF.val + D.val;    result->err  = tmp_CF.err + D.err;    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_ERROR_SELECT_2(stat_tmp_CF, stat_D);  }  else {    return gamma_inc_Q_CF(a, x, result);  }}

int gsl_sf_bessel_il_scaled_array(const int lmax, const double x, double * result_array){  if(x == 0.0) {    int ell;    result_array[0] = 1.0;    for (ell = lmax; ell >= 1; ell--) {      result_array[ell] = 0.0;    };    return GSL_SUCCESS;  } else {    int ell;    gsl_sf_result r_iellp1;    gsl_sf_result r_iell;    int stat_0 = gsl_sf_bessel_il_scaled_e(lmax+1, x, &r_iellp1);    int stat_1 = gsl_sf_bessel_il_scaled_e(lmax,   x, &r_iell);    double iellp1 = r_iellp1.val;    double iell   = r_iell.val;    double iellm1;    result_array[lmax] = iell;    for(ell = lmax; ell >= 1; ell--) {      iellm1 = iellp1 + (2*ell + 1)/x * iell;      iellp1 = iell;      iell   = iellm1;      result_array[ell-1] = iellm1;    }    return GSL_ERROR_SELECT_2(stat_0, stat_1);  }}

int gsl_sf_Ci_e(const double x, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(x <= 0.0) {    DOMAIN_ERROR(result);  }  else if(x <= 4.0) {    const double lx = log(x);    const double y  = (x*x-8.0)*0.125;    gsl_sf_result result_c;    cheb_eval_e(&ci_cs, y, &result_c);    result->val  = lx - 0.5 + result_c.val;    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(lx) + 0.5) + result_c.err;    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else {    gsl_sf_result sin_result;    gsl_sf_result cos_result;    int stat_sin = gsl_sf_sin_e(x, &sin_result);    int stat_cos = gsl_sf_cos_e(x, &cos_result);    gsl_sf_result f;    gsl_sf_result g;    fg_asymp(x, &f, &g);    result->val  = f.val*sin_result.val - g.val*cos_result.val;    result->err  = fabs(f.err*sin_result.val);    result->err += fabs(g.err*cos_result.val);    result->err += fabs(f.val*sin_result.err);    result->err += fabs(g.val*cos_result.err);    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_ERROR_SELECT_2(stat_sin, stat_cos);  }}

intgsl_sf_bessel_Knu_e(const double nu, const double x, gsl_sf_result * result){  gsl_sf_result b;  int stat_K = gsl_sf_bessel_Knu_scaled_e(nu, x, &b);  int stat_e = gsl_sf_exp_mult_err_e(-x, 0.0, b.val, b.err, result);  return GSL_ERROR_SELECT_2(stat_e, stat_K);}

intgsl_sf_hyperg_2F1_conj_renorm_e(const double aR, const double aI, const double c,                                   const double x,                                   gsl_sf_result * result                                   ){  const double rintc = floor(c  + 0.5);  const double rinta = floor(aR + 0.5);  const int a_neg_integer = ( aR < 0.0 && fabs(aR-rinta) < locEPS && aI == 0.0);  const int c_neg_integer = (  c < 0.0 && fabs(c - rintc) < locEPS );  if(c_neg_integer) {    if(a_neg_integer && aR > c+0.1) {      /* 2F1 terminates early */      result->val = 0.0;      result->err = 0.0;      return GSL_SUCCESS;    }    else {      /* 2F1 does not terminate early enough, so something survives */      /* [Abramowitz+Stegun, 15.1.2] */      gsl_sf_result g1, g2;      gsl_sf_result g3;      gsl_sf_result a1, a2;      int stat = 0;      stat += gsl_sf_lngamma_complex_e(aR-c+1, aI, &g1, &a1);      stat += gsl_sf_lngamma_complex_e(aR, aI, &g2, &a2);      stat += gsl_sf_lngamma_e(-c+2.0, &g3);      if(stat != 0) {        DOMAIN_ERROR(result);      }      else {        gsl_sf_result F;        int stat_F = gsl_sf_hyperg_2F1_conj_e(aR-c+1, aI, -c+2, x, &F);        double ln_pre_val = 2.0*(g1.val - g2.val) - g3.val;        double ln_pre_err = 2.0 * (g1.err + g2.err) + g3.err;        int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,                                              F.val, F.err,                                              result);        return GSL_ERROR_SELECT_2(stat_e, stat_F);      }    }  }  else {    /* generic c */    gsl_sf_result F;    gsl_sf_result lng;    double sgn;    int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);    int stat_F = gsl_sf_hyperg_2F1_conj_e(aR, aI, c, x, &F);    int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,                                          sgn*F.val, F.err,                                          result);    return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);  }}

int gsl_sf_bessel_yl_e(int l, const double x, gsl_sf_result * result){    /* CHECK_POINTER(result) */    if(l < 0 || x <= 0.0) {        DOMAIN_ERROR(result);    }    else if(l == 0) {        return gsl_sf_bessel_y0_e(x, result);    }    else if(l == 1) {        return gsl_sf_bessel_y1_e(x, result);    }    else if(l == 2) {        return gsl_sf_bessel_y2_e(x, result);    }    else if(x < 3.0) {        return bessel_yl_small_x(l, x, result);    }    else if(GSL_ROOT3_DBL_EPSILON * x > (l*l + l + 1.0)) {        int status = gsl_sf_bessel_Ynu_asympx_e(l + 0.5, x, result);        double pre = sqrt((0.5*M_PI)/x);        result->val *= pre;        result->err *= pre;        return status;    }    else if(l > 40) {        int status = gsl_sf_bessel_Ynu_asymp_Olver_e(l + 0.5, x, result);        double pre = sqrt((0.5*M_PI)/x);        result->val *= pre;        result->err *= pre;        return status;    }    else {        /* recurse upward */        gsl_sf_result r_by;        gsl_sf_result r_bym;        int stat_1 = gsl_sf_bessel_y1_e(x, &r_by);        int stat_0 = gsl_sf_bessel_y0_e(x, &r_bym);        double bym = r_bym.val;        double by  = r_by.val;        double byp;        int j;        for(j=1; j<l; j++) {            byp = (2*j+1)/x*by - bym;            bym = by;            by  = byp;        }        result->val = by;        result->err = fabs(result->val) * (GSL_DBL_EPSILON + fabs(r_by.err/r_by.val) + fabs(r_bym.err/r_bym.val));        return GSL_ERROR_SELECT_2(stat_1, stat_0);    }}

int gsl_sf_bessel_kl_scaled_e(int l, const double x, gsl_sf_result * result){  if(l < 0 || x <= 0.0) {    DOMAIN_ERROR(result);  }  else if(l == 0) {    return gsl_sf_bessel_k0_scaled_e(x, result);  }  else if(l == 1) {    return gsl_sf_bessel_k1_scaled_e(x, result);  }  else if(l == 2) {    return gsl_sf_bessel_k2_scaled_e(x, result);  }  else if(x < 3.0) {    return bessel_kl_scaled_small_x(l, x, result);  }  else if(GSL_ROOT3_DBL_EPSILON * x > (l*l + l + 1)) {    int status = gsl_sf_bessel_Knu_scaled_asympx_e(l + 0.5, x, result);    double pre = sqrt((0.5*M_PI)/x);    result->val *= pre;    result->err *= pre;    return status;  }  else if(GSL_MIN(0.29/(l*l+1.0), 0.5/(l*l+1.0+x*x)) < GSL_ROOT3_DBL_EPSILON) {    int status = gsl_sf_bessel_Knu_scaled_asymp_unif_e(l + 0.5, x, result);    double pre = sqrt((0.5*M_PI)/x);    result->val *= pre;    result->err *= pre;    return status;  }  else {    /* recurse upward */    gsl_sf_result r_bk;    gsl_sf_result r_bkm;    int stat_1 = gsl_sf_bessel_k1_scaled_e(x, &r_bk);    int stat_0 = gsl_sf_bessel_k0_scaled_e(x, &r_bkm);    double bkp;    double bk  = r_bk.val;    double bkm = r_bkm.val;    int j;    for(j=1; j<l; j++) {       bkp = (2*j+1)/x*bk + bkm;      bkm = bk;      bk  = bkp;    }    result->val  = bk;    result->err  = fabs(bk) * (fabs(r_bk.err/r_bk.val) + fabs(r_bkm.err/r_bkm.val));    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_ERROR_SELECT_2(stat_1, stat_0);  }}

intgsl_sf_bessel_K_scaled_temme(const double nu, const double x,                             double * K_nu, double * K_nup1, double * Kp_nu){  const int max_iter = 15000;  const double half_x    = 0.5 * x;  const double ln_half_x = log(half_x);  const double half_x_nu = exp(nu*ln_half_x);  const double pi_nu   = M_PI * nu;  const double sigma   = -nu * ln_half_x;  const double sinrat  = (fabs(pi_nu) < GSL_DBL_EPSILON ? 1.0 : pi_nu/sin(pi_nu));  const double sinhrat = (fabs(sigma) < GSL_DBL_EPSILON ? 1.0 : sinh(sigma)/sigma);  const double ex = exp(x);  double sum0, sum1;  double fk, pk, qk, hk, ck;  int k = 0;  int stat_iter;  double g_1pnu, g_1mnu, g1, g2;  int stat_g = gsl_sf_temme_gamma(nu, &g_1pnu, &g_1mnu, &g1, &g2);  fk = sinrat * (cosh(sigma)*g1 - sinhrat*ln_half_x*g2);  pk = 0.5/half_x_nu * g_1pnu;  qk = 0.5*half_x_nu * g_1mnu;  hk = pk;  ck = 1.0;  sum0 = fk;  sum1 = hk;  while(k < max_iter) {    double del0;    double del1;    k++;    fk  = (k*fk + pk + qk)/(k*k-nu*nu);    ck *= half_x*half_x/k;    pk /= (k - nu);    qk /= (k + nu);    hk  = -k*fk + pk;    del0 = ck * fk;    del1 = ck * hk;    sum0 += del0;    sum1 += del1;    if(fabs(del0) < 0.5*fabs(sum0)*GSL_DBL_EPSILON) break;  }    *K_nu   = sum0 * ex;  *K_nup1 = sum1 * 2.0/x * ex;  *Kp_nu  = - *K_nup1 + nu/x * *K_nu;  stat_iter = ( k == max_iter ? GSL_EMAXITER : GSL_SUCCESS );  return GSL_ERROR_SELECT_2(stat_iter, stat_g);}

intgsl_sf_bessel_In_scaled_array(const int nmin, const int nmax, const double x, double * result_array){  /* CHECK_POINTER(result_array) */  if(nmax < nmin || nmin < 0) {    int j;    for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0;    GSL_ERROR ("domain error", GSL_EDOM);  }  else if(x == 0.0) {    int j;    for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0;    if(nmin == 0) result_array[0] = 1.0;    return GSL_SUCCESS;  }  else if(nmax == 0) {    gsl_sf_result I0_scaled;    int stat = gsl_sf_bessel_I0_scaled_e(x, &I0_scaled);    result_array[0] = I0_scaled.val;    return stat;  }  else {    const double ax = fabs(x);    const double two_over_x = 2.0/ax;    /* starting values */    gsl_sf_result r_Inp1;    gsl_sf_result r_In;    int stat_0 = gsl_sf_bessel_In_scaled_e(nmax+1, ax, &r_Inp1);    int stat_1 = gsl_sf_bessel_In_scaled_e(nmax,   ax, &r_In);    double Inp1 = r_Inp1.val;    double In   = r_In.val;    double Inm1;    int n;    for(n=nmax; n>=nmin; n--) {      result_array[n-nmin] = In;      Inm1 = Inp1 + n * two_over_x * In;      Inp1 = In;      In   = Inm1;    }    /* deal with signs */    if(x < 0.0) {      for(n=nmin; n<=nmax; n++) {        if(GSL_IS_ODD(n)) result_array[n-nmin] = -result_array[n-nmin];      }    }    return GSL_ERROR_SELECT_2(stat_0, stat_1);  }}

/* series for small a and x, but not defined for a == 0 */static intgamma_inc_series(double a, double x, gsl_sf_result * result){  gsl_sf_result Q;  gsl_sf_result G;  const int stat_Q = gamma_inc_Q_series(a, x, &Q);  const int stat_G = gsl_sf_gamma_e(a, &G);  result->val = Q.val * G.val;  result->err = fabs(Q.val * G.err) + fabs(Q.err * G.val);  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);  return GSL_ERROR_SELECT_2(stat_Q, stat_G);}

int gsl_sf_bessel_Kn_scaled_array(const int nmin, const int nmax, const double x, double * result_array){  /* CHECK_POINTER(result_array) */  if(nmin < 0 || nmax < nmin || x <= 0.0) {    int j;    for(j=0; j<=nmax-nmin; j++) result_array[j] = 0.0;    GSL_ERROR ("domain error", GSL_EDOM);  }  else if(nmax == 0) {    gsl_sf_result b;    int stat = gsl_sf_bessel_K0_scaled_e(x, &b);    result_array[0] = b.val;    return stat;  }  else {    double two_over_x = 2.0/x;    gsl_sf_result r_Knm1;    gsl_sf_result r_Kn;    int stat_0 = gsl_sf_bessel_Kn_scaled_e(nmin,   x, &r_Knm1);    int stat_1 = gsl_sf_bessel_Kn_scaled_e(nmin+1, x, &r_Kn);    int stat = GSL_ERROR_SELECT_2(stat_0, stat_1);    double Knp1;    double Kn   = r_Kn.val;    double Knm1 = r_Knm1.val;    int n;    for(n=nmin+1; n<=nmax+1; n++) {      if(Knm1 < GSL_DBL_MAX) {        result_array[n-1-nmin] = Knm1;        Knp1 = Knm1 + n * two_over_x * Kn;        Knm1 = Kn;        Kn   = Knp1;      }      else {        /* Overflow. Set the rest of the elements to         * zero and bug out.         * FIXME: Note: this relies on the convention         * that the test x < DBL_MIN fails for x not         * a number. This may be only an IEEE convention,         * so the portability is unclear.         */        int j;        for(j=n; j<=nmax+1; j++) result_array[j-1-nmin] = 0.0;        GSL_ERROR ("overflow", GSL_EOVRFLW);      }    }    return stat;  }}

intgsl_sf_airy_Bi_deriv_e(const double x, gsl_mode_t mode, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(x < -1.0) {    gsl_sf_result a;    gsl_sf_result p;    int status_ap = airy_deriv_mod_phase(x, mode, &a, &p);    double s    = sin(p.val);    result->val  = a.val * s;    result->err  = fabs(result->val * p.err) + fabs(s * a.err);    result->err += GSL_DBL_EPSILON * fabs(result->val);    return status_ap;  }  else if(x < 1.0) {    const double x3 = x*x*x;    const double x2 = x*x;    gsl_sf_result result_c1;    gsl_sf_result result_c2;    cheb_eval_mode_e(&bif_cs, x3, mode, &result_c1);    cheb_eval_mode_e(&big_cs, x3, mode, &result_c2);    result->val  = x2 * (result_c1.val + 0.25) + result_c2.val + 0.5;    result->err  = x2 * result_c1.err + result_c2.err;    result->err += GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else if(x < 2.0) {    const double z = (2.0*x*x*x - 9.0) / 7.0;    gsl_sf_result result_c1;    gsl_sf_result result_c2;    cheb_eval_mode_e(&bif2_cs, z, mode, &result_c1);    cheb_eval_mode_e(&big2_cs, z, mode, &result_c2);    result->val  = x*x * (result_c1.val + 0.25) + result_c2.val + 0.5;    result->err  = x*x * result_c1.err + result_c2.err;    result->err += GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else if(x < GSL_ROOT3_DBL_MAX*GSL_ROOT3_DBL_MAX) {    gsl_sf_result result_bps;    const double arg = 2.0*(x*sqrt(x)/3.0);    int stat_b = gsl_sf_airy_Bi_deriv_scaled_e(x, mode, &result_bps);    int stat_e = gsl_sf_exp_mult_err_e(arg, 1.5*fabs(arg*GSL_DBL_EPSILON),                                          result_bps.val, result_bps.err,                                          result);    return GSL_ERROR_SELECT_2(stat_e, stat_b);  }  else {    OVERFLOW_ERROR(result);  }}

/* Continued fraction for Q. * * Q(a,x) = D(a,x) a/x F(a,x) * * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no): * * Since the Gautschi equivalent series method for CF evaluation may lead * to singularities, I have replaced it with the modified Lentz algorithm * given in * * I J Thompson and A R Barnett * Coulomb and Bessel Functions of Complex Arguments and Order * J Computational Physics 64:490-509 (1986) * * In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been * removed. * * Identification of terms between the above equation for F(a, x) and * the first equation in the appendix of Thompson&Barnett is as follows: * *    b_0 = 0, b_n = 1 for all n > 0 * *    a_1 = 1 *    a_n = (n/2-a)/x    for n even *    a_n = (n-1)/(2x)   for n odd * */staticintgamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result){  gsl_sf_result D;  gsl_sf_result F;  const int stat_D = gamma_inc_D(a, x, &D);  const int stat_F = gamma_inc_F_CF(a, x, &F);  result->val  = D.val * (a/x) * F.val;  result->err  = D.err * fabs((a/x) * F.val) + fabs(D.val * a/x * F.err);  return GSL_ERROR_SELECT_2(stat_F, stat_D);}

static intgamma_inc_CF(double a, double x, gsl_sf_result * result){  gsl_sf_result F;  gsl_sf_result pre;  const int stat_F = gamma_inc_F_CF(a, x, &F);  const int stat_E = gsl_sf_exp_e((a-1.0)*log(x) - x, &pre);  result->val = F.val * pre.val;  result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err);  result->err += (2.0 + fabs(a)) * GSL_DBL_EPSILON * fabs(result->val);  return GSL_ERROR_SELECT_2(stat_F, stat_E);}

int gsl_sf_bessel_Kn_scaled_e(int n, const double x, gsl_sf_result * result){  n = abs(n); /* K(-n, z) = K(n, z) */  /* CHECK_POINTER(result) */  if(x <= 0.0) {    DOMAIN_ERROR(result);  }  else if(n == 0) {    return gsl_sf_bessel_K0_scaled_e(x, result);  }  else if(n == 1) {    return gsl_sf_bessel_K1_scaled_e(x, result);  }  else if(x <= 5.0) {    return bessel_Kn_scaled_small_x(n, x, result);  }  else if(GSL_ROOT3_DBL_EPSILON * x > 0.25 * (n*n + 1)) {    return gsl_sf_bessel_Knu_scaled_asympx_e((double)n, x, result);  }  else if(GSL_MIN(0.29/(n*n), 0.5/(n*n + x*x)) < GSL_ROOT3_DBL_EPSILON) {    return gsl_sf_bessel_Knu_scaled_asymp_unif_e((double)n, x, result);  }  else {    /* Upward recurrence. [Gradshteyn + Ryzhik, 8.471.1] */    double two_over_x = 2.0/x;    gsl_sf_result r_b_jm1;    gsl_sf_result r_b_j;    int stat_0 = gsl_sf_bessel_K0_scaled_e(x, &r_b_jm1);    int stat_1 = gsl_sf_bessel_K1_scaled_e(x, &r_b_j);    double b_jm1 = r_b_jm1.val;    double b_j   = r_b_j.val;    double b_jp1;    int j;    for(j=1; j<n; j++) {      b_jp1 = b_jm1 + j * two_over_x * b_j;      b_jm1 = b_j;      b_j   = b_jp1;     }         result->val  = b_j;    result->err  = n * (fabs(b_j) * (fabs(r_b_jm1.err/r_b_jm1.val) + fabs(r_b_j.err/r_b_j.val)));    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_ERROR_SELECT_2(stat_0, stat_1);  }}

intgsl_sf_bessel_Jn_array(int nmin, int nmax, double x, double * result_array){  /* CHECK_POINTER(result_array) */  if(nmin < 0 || nmax < nmin) {    int n;    for(n=nmax; n>=nmin; n--) {      result_array[n-nmin] = 0.0;    }    GSL_ERROR ("domain error", GSL_EDOM);  }  else if(x == 0.0) {    int n;    for(n=nmax; n>=nmin; n--) {      result_array[n-nmin] = 0.0;    }    if(nmin == 0) result_array[0] = 1.0;    return GSL_SUCCESS;  }  else {    gsl_sf_result r_Jnp1;    gsl_sf_result r_Jn;    int stat_np1 = gsl_sf_bessel_Jn_e(nmax+1, x, &r_Jnp1);    int stat_n   = gsl_sf_bessel_Jn_e(nmax,   x, &r_Jn);    int stat = GSL_ERROR_SELECT_2(stat_np1, stat_n);    double Jnp1 = r_Jnp1.val;    double Jn   = r_Jn.val;    double Jnm1;    int n;    if(stat == GSL_SUCCESS) {      for(n=nmax; n>=nmin; n--) {        result_array[n-nmin] = Jn;        Jnm1 = -Jnp1 + 2.0*n/x * Jn;        Jnp1 = Jn;        Jn   = Jnm1;      }    }    else {      for(n=nmax; n>=nmin; n--) {        result_array[n-nmin] = 0.0;      }    }    return stat;  }}

static intgamma_inc_a_gt_0(double a, double x, gsl_sf_result * result){  /* x > 0 and a > 0; use result for Q */  gsl_sf_result Q;  gsl_sf_result G;  const int stat_Q = gsl_sf_gamma_inc_Q_e(a, x, &Q);  const int stat_G = gsl_sf_gamma_e(a, &G);  result->val = G.val * Q.val;  result->err = fabs(G.val * Q.err) + fabs(G.err * Q.val);  result->err += 2.0*GSL_DBL_EPSILON * fabs(result->val);  return GSL_ERROR_SELECT_2(stat_G, stat_Q);}

intgsl_sf_dilog_e(const double x, gsl_sf_result * result){  if(x >= 0.0) {    return dilog_xge0(x, result);  }  else {    gsl_sf_result d1, d2;    int stat_d1 = dilog_xge0( -x, &d1);    int stat_d2 = dilog_xge0(x*x, &d2);    result->val  = -d1.val + 0.5 * d2.val;    result->err  =  d1.err + 0.5 * d2.err;    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_ERROR_SELECT_2(stat_d1, stat_d2);  }}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.3)] */intgsl_sf_ellint_P_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result){  const double sin_phi  = sin(phi);  const double sin2_phi = sin_phi  * sin_phi;  const double sin3_phi = sin2_phi * sin_phi;  const double x = 1.0 - sin2_phi;  const double y = 1.0 - k*k*sin2_phi;  gsl_sf_result rf;  gsl_sf_result rj;  const int rfstatus = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);  const int rjstatus = gsl_sf_ellint_RJ_e(x, y, 1.0, 1.0 + n*sin2_phi, mode, &rj);  result->val  = sin_phi * rf.val - n/3.0*sin3_phi * rj.val;  result->err  = GSL_DBL_EPSILON * fabs(sin_phi * rf.val);  result->err += n/3.0 * fabs(sin3_phi*rj.err);  return GSL_ERROR_SELECT_2(rfstatus, rjstatus);}

int gsl_sf_bessel_j2_e(const double x, gsl_sf_result * result){  double ax = fabs(x);  /* CHECK_POINTER(result) */    if(x == 0.0) {    result->val = 0.0;    result->err = 0.0;    return GSL_SUCCESS;  }  else if(ax < 4.0*GSL_SQRT_DBL_MIN) {    UNDERFLOW_ERROR(result);  }  else if(ax < 1.3) {    const double y  = x*x;    const double c1 = -1.0/14.0;    const double c2 =  1.0/504.0;    const double c3 = -1.0/33264.0;    const double c4 =  1.0/3459456.0;    const double c5 = -1.0/518918400;    const double c6 =  1.0/105859353600.0;    const double c7 = -1.0/28158588057600.0;    const double c8 =  1.0/9461285587353600.0;    const double c9 = -1.0/3916972233164390400.0;    const double sum = 1.0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*(c7+y*(c8+y*c9))))))));    result->val = y/15.0 * sum;    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else {    gsl_sf_result cos_result;    gsl_sf_result sin_result;    const int stat_cos = gsl_sf_cos_e(x, &cos_result);    const int stat_sin = gsl_sf_sin_e(x, &sin_result);    const double cos_x = cos_result.val;    const double sin_x = sin_result.val;    const double f = (3.0/(x*x) - 1.0);    result->val  = (f * sin_x - 3.0*cos_x/x)/x;    result->err  = fabs(f * sin_result.err/x) + fabs((3.0*cos_result.err/x)/x);    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(f*sin_x/x) + 3.0*fabs(cos_x/(x*x)));    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_ERROR_SELECT_2(stat_cos, stat_sin);  }}