C++ GSL_ERROR_SELECT_2函数代码示例

int
gsl_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);
  }
}

int
gsl_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);
  }
}

static
int
bessel_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.
 */
static
int
hyperg_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;
  }
}

int
gsl_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
 */
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, &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)] */
int
gsl_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(). */
static
int
gamma_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);
  }
}

int
gsl_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);
}

int
gsl_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);
  }
}