C++ cheb_eval_e函数代码示例

static double olver_A3(double z, double abs_zeta)
{
  if(z < 0.9) {
    const double x = 20.0*z/9.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&A3_lt1_cs, x, &c);
    return c.val;
  }
  else if(z < 1.1) {
    double a = 1.0-z;
    const double c0 = -0.000354211971457743840771125759200;
    const double c1 = -0.000312322527890318832782774881353;
    const double c2 =  0.000277947465383133980329617631915;
    const double c3 =  0.000919803044747966977054155192400;
    const double c4 =  0.001147600388275977640983696906320;
    const double c5 =  0.000869239326123625742931772044544;
    const double c6 =  0.000287392257282507334785281718027;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*c6)))));
  }
  else {
    const double x   = 11.0/(5.0*z) - 1.0;
    const double zi2 = 1.0/(z*z);
    gsl_sf_result c;
    cheb_eval_e(&A3_gt1_cs, x, &c);
    return  c.val * zi2*zi2*zi2;
  }
}

static double olver_B2(double z, double abs_zeta)
{
  if(z < 0.8) {
    const double x = 5.0*z/2.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&B2_lt1_cs, x, &c);
    return  c.val / z;
  }
  else if(z <= 1.2) {
    const double a = 1.0-z;
    const double c0 = 0.00055221307672129279005986982501;
    const double c1 = 0.00089586516310476929281129228969;
    const double c2 = 0.00067015003441569770883539158863;
    const double c3 = 0.00010166263361949045682945811828;
    const double c4 = -0.00044086345133806887291336488582;
    const double c5 = -0.00073963081508788743392883072523;
    const double c6 = -0.00076745494377839561259903887331;
    const double c7 = -0.00060829038106040362291568012663;
    const double c8 = -0.00037128707528893496121336168683;
    const double c9 = -0.00014116325105702609866850307176;
    return c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*c9))))))));
  }
  else {
    const double zi = 1.0/z;
    const double x  = 12.0/5.0 * zi - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&B2_gt1_cs, x, &c);
    return c.val * zi*zi*zi;
  }
}

static double olver_B3(double z, double abs_zeta)
{
  if(z < 0.8) {
    const double x = 5.0*z/2.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&B3_lt1_cs, x, &c);
    return c.val;
  }
  else if(z < 1.2) {
    const double a = 1.0-z;
    const double c0 = -0.00047461779655995980754441833105;
    const double c1 = -0.00095572913429464297452176811898;
    const double c2 = -0.00080369634512082892655558133973;
    const double c3 = -0.00000727921669154784138080600339;
    const double c4 =  0.00093162500331581345235746518994;
    const double c5 =  0.00149848796913751497227188612403;
    const double c6 =  0.00148406039675949727870390426462;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*c6)))));
  }
  else {
    const double x   = 12.0/(5.0*z) - 1.0;
    const double zi2 = 1.0/(z*z);
    gsl_sf_result c;
    cheb_eval_e(&B3_gt1_cs, x, &c);
    return  c.val * zi2*zi2*zi2;
  }
}

/* assumes s >= 0 and s != 1.0 */
inline
static int
riemann_zeta_sgt0(double s, gsl_sf_result * result)
{
  if(s < 1.0) {
    gsl_sf_result c;
    cheb_eval_e(&zeta_xlt1_cs, 2.0*s - 1.0, &c);
    result->val = c.val / (s - 1.0);
    result->err = c.err / fabs(s-1.0) + GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(s <= 20.0) {
    double x = (2.0*s - 21.0)/19.0;
    gsl_sf_result c;
    cheb_eval_e(&zeta_xgt1_cs, x, &c);
    result->val = c.val / (s - 1.0);
    result->err = c.err / (s - 1.0) + GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    double f2 = 1.0 - pow(2.0,-s);
    double f3 = 1.0 - pow(3.0,-s);
    double f5 = 1.0 - pow(5.0,-s);
    double f7 = 1.0 - pow(7.0,-s);
    result->val = 1.0/(f2*f3*f5*f7);
    result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

static double olver_A4(double z, double abs_zeta)
{
  if(z < 0.8) {
    const double x = 5.0*z/2.0 - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&A4_lt1_cs, x, &c);
    return c.val;
  }
  else if(z < 1.2) {
    double a = 1.0-z;
    const double c0 =  0.00037819419920177291402661228437;
    const double c1 =  0.00040494390552363233477213857527;
    const double c2 = -0.00045764735528936113047289344569;
    const double c3 = -0.00165361044229650225813161341879;
    const double c4 = -0.00217527517983360049717137015539;
    const double c5 = -0.00152003287866490735107772795537;
    return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*c5))));
  }
  else {
    const double x   = 12.0/(5.0*z) - 1.0;
    const double zi2 = 1.0/(z*z);
    gsl_sf_result c;
    cheb_eval_e(&A4_gt1_cs, x, &c);
    return c.val * zi2*zi2*zi2*zi2;
  }
}

int gsl_sf_bessel_J0_e(const double x, gsl_sf_result * result)
{
  double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y < 2.0*GSL_SQRT_DBL_EPSILON) {
    result->val = 1.0;
    result->err = y*y;
    return GSL_SUCCESS;
  }
  else if(y <= 4.0) {
    return cheb_eval_e(&bj0_cs, 0.125*y*y - 1.0, result);
  }
  else {
    const double z = 32.0/(y*y) - 1.0;
    gsl_sf_result ca;
    gsl_sf_result ct;
    gsl_sf_result cp;
    const int stat_ca = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bm0_cs,  z, &ca);
    const int stat_ct = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bth0_cs, z, &ct);
    const int stat_cp = gsl_sf_bessel_cos_pi4_e(y, ct.val/y, &cp);
    const double sqrty = sqrt(y);
    const double ampl  = (0.75 + ca.val) / sqrty;
    result->val  = ampl * cp.val;
    result->err  = fabs(cp.val) * ca.err/sqrty + fabs(ampl) * cp.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_3(stat_ca, stat_ct, stat_cp);
  }
}

int gsl_sf_fermi_dirac_2_e(const double x, gsl_sf_result * result)
{
  if(x < GSL_LOG_DBL_MIN) {
    UNDERFLOW_ERROR(result);
  }
  else if(x < -1.0) {
    /* series [Goano (6)]
     */
    double ex   = exp(x);
    double term = ex;
    double sum  = term;
    int n;
    for(n=2; n<100 ; n++) {
      double rat = (n-1.0)/n;
      term *= -ex * rat * rat * rat;
      sum  += term;
      if(fabs(term/sum) < GSL_DBL_EPSILON) break;
    }
    result->val = sum;
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
    return GSL_SUCCESS;
  }
  else if(x < 1.0) {
    return cheb_eval_e(&fd_2_a_cs, x, result);
  }
  else if(x < 4.0) {
    double t = 2.0/3.0*(x-1.0) - 1.0;
    return cheb_eval_e(&fd_2_b_cs, t, result);
  }
  else if(x < 10.0) {
    double t = 1.0/3.0*(x-4.0) - 1.0;
    return cheb_eval_e(&fd_2_c_cs, t, result);
  }
  else if(x < 30.0) {
    double t = 0.1*x - 2.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_2_d_cs, t, &c);
    result->val  = c.val * x*x*x;
    result->err  = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < 1.0/GSL_ROOT3_DBL_EPSILON) {
    double t = 60.0/x - 1.0;
    gsl_sf_result c;
    cheb_eval_e(&fd_2_e_cs, t, &c);
    result->val = c.val * x*x*x;
    result->err = c.err * x*x*x + 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < GSL_ROOT3_DBL_MAX) {
    result->val = 1.0/6.0 * x*x*x;
    result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result)
{
  const double xmin    = 2.0 * GSL_DBL_MIN;
  const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON;
  const double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(y < xmin) {
    UNDERFLOW_ERROR(result);
  }
  else if(y < x_small) {
    result->val = 0.5*x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(y <= 3.0) {
    const double ey = exp(-y);
    gsl_sf_result c;
    cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c);
    result->val  = x * ey * (0.875 + c.val);
    result->err  = ey * c.err + y * GSL_DBL_EPSILON * fabs(result->val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(y <= 8.0) {
    const double sy = sqrt(y);
    gsl_sf_result c;
    double b;
    double s;
    cheb_eval_e(&ai1_cs, (48.0/y-11.0)/5.0, &c);
    b = (0.375 + c.val) / sy;
    s = (x > 0.0 ? 1.0 : -1.0);
    result->val  = s * b;
    result->err  = c.err / sy;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sy = sqrt(y);
    gsl_sf_result c;
    double b;
    double s;
    cheb_eval_e(&ai12_cs, 16.0/y-1.0, &c);
    b = (0.375 + c.val) / sy;
    s = (x > 0.0 ? 1.0 : -1.0);
    result->val  = s * b;
    result->err  = c.err / sy;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_Y1_e(const double x, gsl_sf_result * result)
{
  const double two_over_pi = 2.0/M_PI;
  const double xmin = 1.571*GSL_DBL_MIN; /*exp ( amax1(alog(r1mach(1)), -alog(r1mach(2)))+.01) */
  const double x_small = 2.0 * GSL_SQRT_DBL_EPSILON;
  const double xmax    = 1.0/GSL_DBL_EPSILON;

  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < xmin) {
    OVERFLOW_ERROR(result);
  }
  else if(x < x_small) {
    const double lnterm = log(0.5*x);
    gsl_sf_result J1;
    gsl_sf_result c;
    int status = gsl_sf_bessel_J1_e(x, &J1);
    cheb_eval_e(&by1_cs, -1.0, &c);
    result->val = two_over_pi * lnterm * J1.val + (0.5 + c.val)/x;
    result->err = fabs(lnterm) * (fabs(GSL_DBL_EPSILON * J1.val) + J1.err) + c.err/x;
    return status;
  }
  else if(x < 4.0) {
    const double lnterm = log(0.5*x);
    int status;
    gsl_sf_result J1;
    gsl_sf_result c;
    cheb_eval_e(&by1_cs, 0.125*x*x-1.0, &c);
    status = gsl_sf_bessel_J1_e(x, &J1);
    result->val = two_over_pi * lnterm * J1.val + (0.5 + c.val)/x;
    result->err = fabs(lnterm) * (fabs(GSL_DBL_EPSILON * J1.val) + J1.err) + c.err/x;
    return status;
  }
  else if(x < xmax) {
    const double z = 32.0/(x*x) - 1.0;
    gsl_sf_result ca;
    gsl_sf_result ct;
    gsl_sf_result cp;
    const int stat_ca = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bm1_cs,  z, &ca);
    const int stat_ct = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bth1_cs, z, &ct);
    const int stat_cp = gsl_sf_bessel_cos_pi4_e(x, ct.val/x, &cp);
    const double sqrtx = sqrt(x);
    const double ampl  = (0.75 + ca.val) / sqrtx;
    result->val  = -ampl * cp.val;
    result->err  = fabs(cp.val) * ca.err/sqrtx + fabs(ampl) * cp.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_3(stat_ca, stat_ct, stat_cp);
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

static
int
gsl_sf_temme_gamma(const double nu, double * g_1pnu, double * g_1mnu, double * g1, double * g2)
{
  const double anu = fabs(nu);    /* functions are even */
  const double x = 4.0*anu - 1.0;
  gsl_sf_result r_g1;
  gsl_sf_result r_g2;
  cheb_eval_e(&g1_cs, x, &r_g1);
  cheb_eval_e(&g2_cs, x, &r_g2);
  *g1 = r_g1.val;
  *g2 = r_g2.val;
  *g_1mnu = 1.0/(r_g2.val + nu * r_g1.val);
  *g_1pnu = 1.0/(r_g2.val - nu * r_g1.val);
  return GSL_SUCCESS;
}

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_sinc_e(double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  {
    const double ax = fabs(x);

    if(ax < 0.8) {
      /* Do not go to the limit of the fit since
       * there is a zero there and the Chebyshev
       * accuracy will go to zero.
       */
      return cheb_eval_e(&sinc_cs, 2.0*ax-1.0, result);
    }
    else if(ax < 100.0) {
      /* Small arguments are no problem.
       * We trust the library sin() to
       * roughly machine precision.
       */
      result->val = sin(M_PI * ax)/(M_PI * ax);
      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
    else {
      /* Large arguments must be handled separately.
       */
      const double r = M_PI*ax;
      gsl_sf_result s;
      int stat_s = gsl_sf_sin_e(r, &s);
      result->val = s.val/r;
      result->err = s.err/r + 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_s;
    }
  }
}

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

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

int gsl_sf_Si_e(const double x, gsl_sf_result * result)
{
  double ax = fabs(x);
  
  /* CHECK_POINTER(result) */

  if(ax < GSL_SQRT_DBL_EPSILON) {
    result->val = x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(ax <= 4.0) {
    gsl_sf_result result_c;
    cheb_eval_e(&si_cs, (x*x-8.0)*0.125, &result_c);
    result->val  =  x * (0.75 + result_c.val);
    result->err  = ax * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    /* Note there is no loss of precision
     * here bcause of the leading constant.
     */
    gsl_sf_result f;
    gsl_sf_result g;
    fg_asymp(ax, &f, &g);
    result->val  = 0.5 * M_PI - f.val*cos(ax) - g.val*sin(ax);
    result->err  = f.err + g.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    if(x < 0.0) result->val = -result->val;
    return GSL_SUCCESS;
  }
}