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

intgsl_sf_ellint_RC_e(double x, double y, gsl_mode_t mode, gsl_sf_result * result){  const double lolim = 5.0 * GSL_DBL_MIN;  const double uplim = 0.2 * GSL_DBL_MAX;  const gsl_prec_t goal = GSL_MODE_PREC(mode);  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );  const double prec   = gsl_prec_eps[goal];  if(x < 0.0 || y < 0.0 || x + y < lolim) {    DOMAIN_ERROR(result);  }  else if(GSL_MAX(x, y) < uplim) {     const double c1 = 1.0 / 7.0;    const double c2 = 9.0 / 22.0;    double xn = x;    double yn = y;    double mu, sn, lamda, s;    while(1) {      mu = (xn + yn + yn) / 3.0;      sn = (yn + mu) / mu - 2.0;      if (fabs(sn) < errtol) break;      lamda = 2.0 * sqrt(xn) * sqrt(yn) + yn;      xn = (xn + lamda) * 0.25;      yn = (yn + lamda) * 0.25;    }    s = sn * sn * (0.3 + sn * (c1 + sn * (0.375 + sn * c2)));    result->val = (1.0 + s) / sqrt(mu);    result->err = prec * fabs(result->val);    return GSL_SUCCESS;  }  else {    DOMAIN_ERROR(result);  }}

intgsl_sf_ellint_RD_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result){  const gsl_prec_t goal = GSL_MODE_PREC(mode);  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );  const double prec   = gsl_prec_eps[goal];  const double lolim = 2.0/pow(GSL_DBL_MAX, 2.0/3.0);  const double uplim = pow(0.1*errtol/GSL_DBL_MIN, 2.0/3.0);  if(GSL_MIN(x,y) < 0.0 || GSL_MIN(x+y,z) < lolim) {    DOMAIN_ERROR(result);  }  else if(locMAX3(x,y,z) < uplim) {    const double c1 = 3.0 / 14.0;    const double c2 = 1.0 /  6.0;    const double c3 = 9.0 / 22.0;    const double c4 = 3.0 / 26.0;    double xn = x;    double yn = y;    double zn = z;    double sigma  = 0.0;    double power4 = 1.0;    double ea, eb, ec, ed, ef, s1, s2;    double mu, xndev, yndev, zndev;    while(1) {      double xnroot, ynroot, znroot, lamda;      double epslon;      mu = (xn + yn + 3.0 * zn) * 0.2;      xndev = (mu - xn) / mu;      yndev = (mu - yn) / mu;      zndev = (mu - zn) / mu;      epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));      if (epslon < errtol) break;      xnroot = sqrt(xn);      ynroot = sqrt(yn);      znroot = sqrt(zn);      lamda = xnroot * (ynroot + znroot) + ynroot * znroot;      sigma  += power4 / (znroot * (zn + lamda));      power4 *= 0.25;      xn = (xn + lamda) * 0.25;      yn = (yn + lamda) * 0.25;      zn = (zn + lamda) * 0.25;    }    ea = xndev * yndev;    eb = zndev * zndev;    ec = ea - eb;    ed = ea - 6.0 * eb;    ef = ed + ec + ec;    s1 = ed * (- c1 + 0.25 * c3 * ed - 1.5 * c4 * zndev * ef);    s2 = zndev * (c2 * ef + zndev * (- c3 * ec + zndev * c4 * ea));    result->val = 3.0 * sigma + power4 * (1.0 + s1 + s2) / (mu * sqrt(mu));    result->err = prec * fabs(result->val);    return GSL_SUCCESS;  }  else {    DOMAIN_ERROR(result);  }}

static int get_trait(SEXP variable, struct design *s, struct design *r, struct variable *v){	SEXP dim, cn, dn, nm;	int j, n, p;	double *xt, *x;	size_t dims[VAR_RANK_MAX];	size_t rank;	dim = GET_DIM(variable);	n = INTEGER(dim)[0];	p = INTEGER(dim)[1];	xt = NUMERIC_POINTER(variable);	x = (void *)R_alloc(n * p, sizeof(*x));	if (p == 0)		DOMAIN_ERROR("variable dimensions must be nonzero");	if (n > 0)		matrix_dtrans(n, p, xt, n, x, p);	if (p == 1) {		rank = 0;	} else {		rank = 1;		dims[0] = (size_t)p;	}	if (inherits(variable, "send")) {		if (n != design_count(s))			DOMAIN_ERROR("variable dimension must match sender count");		v->design = VARIABLE_DESIGN_SEND;		v->var.send = design_add_trait(s, NULL, x, dims, rank);	} else {		if (n != design_count(r))			DOMAIN_ERROR("variable dimension must match receiver count");		v->design = VARIABLE_DESIGN_RECV;		v->var.recv = design_add_trait(r, NULL, x, dims, rank);	}	v->type = VARIABLE_TYPE_TRAIT;	dn = GET_DIMNAMES(variable);	cn = GET_COLNAMES(dn);	v->names = (void *)R_alloc(p, sizeof(*v->names));	for (j = 0; j < p; j++) {		nm = STRING_ELT(cn, j);		v->names[j] = CHAR(nm);	}	return 0;}

intgsl_sf_ellint_RF_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result){  const double lolim = 5.0 * GSL_DBL_MIN;  const double uplim = 0.2 * GSL_DBL_MAX;  const gsl_prec_t goal = GSL_MODE_PREC(mode);  const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );  const double prec   = gsl_prec_eps[goal];  if(x < 0.0 || y < 0.0 || z < 0.0) {    DOMAIN_ERROR(result);  }  else if(x+y < lolim || x+z < lolim || y+z < lolim) {    DOMAIN_ERROR(result);  }  else if(locMAX3(x,y,z) < uplim) {     const double c1 = 1.0 / 24.0;    const double c2 = 3.0 / 44.0;    const double c3 = 1.0 / 14.0;    double xn = x;    double yn = y;    double zn = z;    double mu, xndev, yndev, zndev, e2, e3, s;    while(1) {      double epslon, lamda;      double xnroot, ynroot, znroot;      mu = (xn + yn + zn) / 3.0;      xndev = 2.0 - (mu + xn) / mu;      yndev = 2.0 - (mu + yn) / mu;      zndev = 2.0 - (mu + zn) / mu;      epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));      if (epslon < errtol) break;      xnroot = sqrt(xn);      ynroot = sqrt(yn);      znroot = sqrt(zn);      lamda = xnroot * (ynroot + znroot) + ynroot * znroot;      xn = (xn + lamda) * 0.25;      yn = (yn + lamda) * 0.25;      zn = (zn + lamda) * 0.25;    }    e2 = xndev * yndev - zndev * zndev;    e3 = xndev * yndev * zndev;    s = 1.0 + (c1 * e2 - 0.1 - c2 * e3) * e2 + c3 * e3;    result->val = s / sqrt(mu);    result->err = prec * fabs(result->val);    return GSL_SUCCESS;  }  else {    DOMAIN_ERROR(result);  }}

int gsl_sf_zeta_int_e(const int n, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(n < 0) {    if(!GSL_IS_ODD(n)) {      result->val = 0.0; /* exactly zero at even negative integers */      result->err = 0.0;      return GSL_SUCCESS;    }    else if(n > -ZETA_NEG_TABLE_NMAX) {      result->val = zeta_neg_int_table[-(n+1)/2];      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);      return GSL_SUCCESS;    }    else {      return gsl_sf_zeta_e((double)n, result);    }  }  else if(n == 1){    DOMAIN_ERROR(result);  }  else if(n <= ZETA_POS_TABLE_NMAX){    result->val = 1.0 + zetam1_pos_int_table[n];    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else {    result->val = 1.0;    result->err = GSL_DBL_EPSILON;    return GSL_SUCCESS;  }}

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_psi_int_e(const int n, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(n <= 0) {    DOMAIN_ERROR(result);  }  else if(n <= PSI_TABLE_NMAX) {    result->val = psi_table[n];    result->err = GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else {    /* Abramowitz+Stegun 6.3.18 */    const double c2 = -1.0/12.0;    const double c3 =  1.0/120.0;    const double c4 = -1.0/252.0;    const double c5 =  1.0/240.0;    const double ni2 = (1.0/n)*(1.0/n);    const double ser = ni2 * (c2 + ni2 * (c3 + ni2 * (c4 + ni2*c5)));    result->val  = log(n) - 0.5/n + ser;    result->err  = GSL_DBL_EPSILON * (fabs(log(n)) + fabs(0.5/n) + fabs(ser));    result->err += GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }}

int gsl_sf_psi_1_int_e(const int n, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(n <= 0) {    DOMAIN_ERROR(result);  }  else if(n <= PSI_1_TABLE_NMAX) {    result->val = psi_1_table[n];    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else {    /* Abramowitz+Stegun 6.4.12     * double-precision for n > 100     */    const double c0 = -1.0/30.0;    const double c1 =  1.0/42.0;    const double c2 = -1.0/30.0;    const double ni2 = (1.0/n)*(1.0/n);    const double ser =  ni2*ni2 * (c0 + ni2*(c1 + c2*ni2));    result->val = (1.0 + 0.5/n + 1.0/(6.0*n*n) + ser) / n;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }}

intgsl_sf_hyperg_2F0_e(const double a, const double b, const double x, gsl_sf_result * result){  if(x < 0.0) {    /* Use "definition" 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x).     */    gsl_sf_result U;    double pre = pow(-1.0/x, a);    int stat_U = gsl_sf_hyperg_U_e(a, 1.0+a-b, -1.0/x, &U);    result->val = pre * U.val;    result->err = GSL_DBL_EPSILON * fabs(result->val) + pre * U.err;    return stat_U;  }  else if(x == 0.0) {    result->val = 1.0;    result->err = 0.0;    return GSL_SUCCESS;  }  else {    /* Use asymptotic series. ??     */    /* return hyperg_2F0_series(a, b, x, -1, result, &prec); */    DOMAIN_ERROR(result);  }}

intgsl_sf_legendre_H3d_0_e(const double lambda, const double eta, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(eta < 0.0) {    DOMAIN_ERROR(result);  }  else if(eta == 0.0 || lambda == 0.0) {    result->val = 1.0;    result->err = 0.0;    return GSL_SUCCESS;  }  else {    const double lam_eta = lambda * eta;    gsl_sf_result s;    gsl_sf_sin_err_e(lam_eta, 2.0*GSL_DBL_EPSILON * fabs(lam_eta), &s);    if(eta > -0.5*GSL_LOG_DBL_EPSILON) {      double f = 2.0 / lambda * exp(-eta);      result->val  = f * s.val;      result->err  = fabs(f * s.val) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;      result->err += fabs(f) * s.err;      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    }    else {      double f = 1.0/(lambda*sinh(eta));      result->val  = f * s.val;      result->err  = fabs(f * s.val) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;      result->err += fabs(f) * s.err;      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);    }    return GSL_SUCCESS;  }}

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

/* [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_psi_n_e(const int n, const double x, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(n == 0)  {    return gsl_sf_psi_e(x, result);  }  else if(n == 1)  {    return gsl_sf_psi_1_e(x, result);  }  else if(n < 0 || x <= 0.0) {    DOMAIN_ERROR(result);  }  else {    gsl_sf_result ln_nf;    gsl_sf_result hzeta;    int stat_hz = gsl_sf_hzeta_e(n+1.0, x, &hzeta);    int stat_nf = gsl_sf_lnfact_e((unsigned int) n, &ln_nf);    int stat_e  = gsl_sf_exp_mult_err_e(ln_nf.val, ln_nf.err,                                           hzeta.val, hzeta.err,                                           result);    if(GSL_IS_EVEN(n)) result->val = -result->val;    return GSL_ERROR_SELECT_3(stat_e, stat_nf, stat_hz);  }}

intgsl_sf_lnsinh_e(const double x, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(x <= 0.0) {    DOMAIN_ERROR(result);  }  else if(fabs(x) < 1.0) {    double eps;    sinh_series(x, &eps);    result->val = log(eps);    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else if(x < -0.5*GSL_LOG_DBL_EPSILON) {    result->val = x + log(0.5*(1.0 - exp(-2.0*x)));    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else {    result->val = -M_LN2 + x;    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }}

int gsl_sf_zetam1_int_e(const int n, gsl_sf_result * result){  if(n < 0) {    if(!GSL_IS_ODD(n)) {      result->val = -1.0; /* at even negative integers zetam1 == -1 since zeta is exactly zero */      result->err = 0.0;      return GSL_SUCCESS;    }    else if(n > -ZETA_NEG_TABLE_NMAX) {      result->val = zeta_neg_int_table[-(n+1)/2] - 1.0;      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);      return GSL_SUCCESS;    }    else {      /* could use gsl_sf_zetam1_e here but subtracting 1 makes no difference         for such large values, so go straight to the result */      return gsl_sf_zeta_e((double)n, result);      }  }  else if(n == 1){    DOMAIN_ERROR(result);  }  else if(n <= ZETA_POS_TABLE_NMAX){    result->val = zetam1_pos_int_table[n];    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else {    return gsl_sf_zetam1_e(n, result);  }}

static int get_variables(SEXP variables,			 struct design *s, struct design *r, struct design2 *d,			 struct variables *v){	SEXP var;	int i, n;	int err = 0;	if (!IS_VECTOR(variables))		DOMAIN_ERROR("'variables' should be a list");	n = LENGTH(variables);	v->count = (size_t)n;	v->item = (void *)R_alloc(n, sizeof(*v->item));	for (i = 0; i < n; i++) {		var = VECTOR_ELT(variables, i);		err = get_variable(var, s, r, d, &v->item[i]);		if (err < 0)			goto out;	}out:	return err;}

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

int gsl_sf_hzeta_e(const double s, const double q, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(s <= 1.0 || q <= 0.0) {    DOMAIN_ERROR(result);  }  else {    const double max_bits = 54.0;    const double ln_term0 = -s * log(q);      if(ln_term0 < GSL_LOG_DBL_MIN + 1.0) {      UNDERFLOW_ERROR(result);    }    else if(ln_term0 > GSL_LOG_DBL_MAX - 1.0) {      OVERFLOW_ERROR (result);    }    else if((s > max_bits && q < 1.0) || (s > 0.5*max_bits && q < 0.25)) {      result->val = pow(q, -s);      result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);      return GSL_SUCCESS;    }    else if(s > 0.5*max_bits && q < 1.0) {      const double p1 = pow(q, -s);      const double p2 = pow(q/(1.0+q), s);      const double p3 = pow(q/(2.0+q), s);      result->val = p1 * (1.0 + p2 + p3);      result->err = GSL_DBL_EPSILON * (0.5*s + 2.0) * fabs(result->val);      return GSL_SUCCESS;    }    else {      /* Euler-Maclaurin summation formula        * [Moshier, p. 400, with several typo corrections]       */      const int jmax = 12;      const int kmax = 10;      int j, k;      const double pmax  = pow(kmax + q, -s);      double scp = s;      double pcp = pmax / (kmax + q);      double ans = pmax*((kmax+q)/(s-1.0) + 0.5);      for(k=0; k<kmax; k++) {        ans += pow(k + q, -s);      }      for(j=0; j<=jmax; j++) {        double delta = hzeta_c[j+1] * scp * pcp;        ans += delta;        if(fabs(delta/ans) < 0.5*GSL_DBL_EPSILON) break;        scp *= (s+2*j+1)*(s+2*j+2);        pcp /= (kmax + q)*(kmax + q);      }      result->val = ans;      result->err = 2.0 * (jmax + 1.0) * GSL_DBL_EPSILON * fabs(ans);      return GSL_SUCCESS;    }  }}

intgsl_sf_bessel_lnKnu_e(const double nu, const double x, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(x <= 0.0 || nu < 0.0) {    DOMAIN_ERROR(result);  }  else if(nu == 0.0) {    gsl_sf_result K_scaled;    /* This cannot underflow, and     * it will not throw GSL_EDOM     * since that is already checked.     */    gsl_sf_bessel_K0_scaled_e(x, &K_scaled);    result->val  = -x + log(fabs(K_scaled.val));    result->err  = GSL_DBL_EPSILON * fabs(x) + fabs(K_scaled.err/K_scaled.val);    result->err += GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }  else if(x < 2.0 && nu > 1.0) {    /* Make use of the inequality     * Knu(x) <= 1/2 (2/x)^nu Gamma(nu),     * which follows from the integral representation     * [Abramowitz+Stegun, 9.6.23 (2)]. With this     * we decide whether or not there is an overflow     * problem because x is small.     */    double ln_bound;    gsl_sf_result lg_nu;    gsl_sf_lngamma_e(nu, &lg_nu);    ln_bound = -M_LN2 - nu*log(0.5*x) + lg_nu.val;    if(ln_bound > GSL_LOG_DBL_MAX - 20.0) {      /* x must be very small or nu very large (or both).       */      double xi  = 0.25*x*x;      double sum = 1.0 - xi/(nu-1.0);      if(nu > 2.0) sum +=  (xi/(nu-1.0)) * (xi/(nu-2.0));      result->val  = ln_bound + log(sum);      result->err  = lg_nu.err;      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);      return GSL_SUCCESS;    }    /* can drop-through here */  }  {    /* We passed the above tests, so no problem.     * Evaluate as usual. Note the possible drop-through     * in the above code!     */    gsl_sf_result K_scaled;    gsl_sf_bessel_Knu_scaled_e(nu, x, &K_scaled);    result->val  = -x + log(fabs(K_scaled.val));    result->err  = GSL_DBL_EPSILON * fabs(x) + fabs(K_scaled.err/K_scaled.val);    result->err += GSL_DBL_EPSILON * fabs(result->val);    return GSL_SUCCESS;  }}

/* [Carlson, Numer. Math. 33 (1979) 1, (4.5)] */intgsl_sf_ellint_Kcomp_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.33] */    const double y = 1.0 - k*k;    const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };    const double b[] = { 0.5, 0.12498593597, 0.06880248576 };    const double ta = a[0] + y*(a[1] + y*a[2]);    const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));    result->val = ta + tb;    result->err = 2.0 * GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else {    /* This was previously computed as,         return gsl_sf_ellint_RF_e(0.0, 1.0 - k*k, 1.0, mode, result);       but this underestimated the total error for small k, since the        argument y=1-k^2 is not exact (there is an absolute error of       GSL_DBL_EPSILON near y=0 due to cancellation in the subtraction).       Taking the singular behavior of -log(y) above gives an error       of 0.5*epsilon/y near y=0. (BJG) */    double y = 1.0 - k*k;    int status = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, result);    result->err += 0.5 * GSL_DBL_EPSILON / y;    return status ;  }}

intgsl_sf_hyperg_1F1_series_e(const double a, const double b, const double x,                           gsl_sf_result * result                          ){    double an  = a;    double bn  = b;    double n   = 1.0;    double del = 1.0;    double abs_del = 1.0;    double max_abs_del = 1.0;    double sum_val = 1.0;    double sum_err = 0.0;    while(abs_del/fabs(sum_val) > GSL_DBL_EPSILON) {        double u, abs_u;        if(bn == 0.0) {            DOMAIN_ERROR(result);        }        if(an == 0.0 || n > 1000.0) {            result->val  = sum_val;            result->err  = sum_err;            result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);            return GSL_SUCCESS;        }        u = x * (an/(bn*n));        abs_u = fabs(u);        if(abs_u > 1.0 && max_abs_del > GSL_DBL_MAX/abs_u) {            result->val = sum_val;            result->err = fabs(sum_val);            GSL_ERROR ("overflow", GSL_EOVRFLW);        }        del *= u;        sum_val += del;        if(fabs(sum_val) > SUM_LARGE) {            result->val = sum_val;            result->err = fabs(sum_val);            GSL_ERROR ("overflow", GSL_EOVRFLW);        }        abs_del = fabs(del);        max_abs_del = GSL_MAX_DBL(abs_del, max_abs_del);        sum_err += 2.0*GSL_DBL_EPSILON*abs_del;        an += 1.0;        bn += 1.0;        n  += 1.0;    }    result->val  = sum_val;    result->err  = sum_err;    result->err += abs_del;    result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);    return GSL_SUCCESS;}

intgsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result * result){  if(x < 0.0) {    DOMAIN_ERROR(result);  }  else if(x == 0.0) {    return gsl_sf_gamma_e(a, result);  }  else if(a == 0.0)  {    return GAMMA_INC_A_0(x, result);  }  else if(a > 0.0)  {    return gamma_inc_a_gt_0(a, x, result);  }  else if(x > 0.25)  {    /* continued fraction seems to fail for x too small; otherwise       it is ok, independent of the value of |x/a|, because of the       non-oscillation in the expansion, i.e. the CF is       un-conditionally convergent for a < 0 and x > 0     */    return gamma_inc_CF(a, x, result);  }  else if(fabs(a) < 0.5)  {    return gamma_inc_series(a, x, result);  }  else  {    /* a = fa + da; da >= 0 */    const double fa = floor(a);    const double da = a - fa;    gsl_sf_result g_da;    const int stat_g_da = ( da > 0.0 ? gamma_inc_a_gt_0(da, x, &g_da)                                     : GAMMA_INC_A_0(x, &g_da));    double alpha = da;    double gax = g_da.val;    /* Gamma(alpha-1,x) = 1/(alpha-1) (Gamma(a,x) - x^(alpha-1) e^-x) */    do    {      const double shift = exp(-x + (alpha-1.0)*log(x));      gax = (gax - shift) / (alpha - 1.0);      alpha -= 1.0;    } while(alpha > a);    result->val = gax;    result->err = 2.0*(1.0 + fabs(a))*GSL_DBL_EPSILON*fabs(gax);    return stat_g_da;  }}

int gsl_sf_debye_5_e(const double x, gsl_sf_result * result){  const double val_infinity = 610.405837190669483828710757875 ;  const double xcut = -GSL_LOG_DBL_MIN;  /* CHECK_POINTER(result) */  if(x < 0.0) {    DOMAIN_ERROR(result);  }  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {    result->val = 1.0 - 5.0*x/12.0 + 5.0*x*x/84.0;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else if(x <= 4.0) {    const double t = x*x/8.0 - 1.0;    gsl_sf_result c;    cheb_eval_e(&adeb5_cs, t, &c);    result->val = c.val - 5.0*x/12.0;    result->err = c.err + GSL_DBL_EPSILON * 5.0*x/12.0;    return GSL_SUCCESS;  }  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {    const int nexp = floor(xcut/x);    const double ex  = exp(-x);    double xk  = nexp * x;    double rk  = nexp;    double sum = 0.0;    int i;    for(i=nexp; i>=1; i--) {      double xk_inv = 1.0/xk;      sum *= ex;      sum += (((((120.0*xk_inv + 120.0)*xk_inv + 60.0)*xk_inv + 20.0)*xk_inv + 5.0)*xk_inv+ 1.0) / rk;      rk -= 1.0;      xk -= x;    }    result->val = val_infinity/(x*x*x*x*x) - 5.0 * sum * ex;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else if(x < xcut) {    const double x2 = x*x;    const double x4 = x2*x2;    const double x5 = x4*x;    const double sum = 120.0 + 120.0*x + 60.0*x2 + 20.0*x2*x + 5.0*x4 + x5;    result->val = (val_infinity - 5.0 * sum * exp(-x)) / x5;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else {    result->val = ((((val_infinity/x)/x)/x)/x)/x;    result->err = GSL_DBL_EPSILON * result->val;    CHECK_UNDERFLOW(result);    return GSL_SUCCESS;  }}

int gsl_sf_debye_6_e(const double x, gsl_sf_result * result){  const double val_infinity = 4356.06887828990661194792541535 ;  const double xcut = -GSL_LOG_DBL_MIN;  /* CHECK_POINTER(result) */  if(x < 0.0) {    DOMAIN_ERROR(result);  }  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {    result->val = 1.0 - 3.0*x/7.0 + x*x/16.0;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else if(x <= 4.0) {    const double t = x*x/8.0 - 1.0;    gsl_sf_result c;    cheb_eval_e(&adeb6_cs, t, &c);    result->val = c.val - 3.0*x/7.0;    result->err = c.err + GSL_DBL_EPSILON * 3.0*x/7.0;    return GSL_SUCCESS;  }  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {    const int nexp = floor(xcut/x);    const double ex  = exp(-x);    double xk  = nexp * x;    double rk  = nexp;    double sum = 0.0;    int i;    for(i=nexp; i>=1; i--) {      double xk_inv = 1.0/xk;      sum *= ex;      sum += ((((((720.0*xk_inv + 720.0)*xk_inv + 360.0)*xk_inv + 120.0)*xk_inv + 30.0)*xk_inv+ 6.0)*xk_inv+ 1.0) / rk;      rk -= 1.0;      xk -= x;    }    result->val = val_infinity/(x*x*x*x*x*x) - 6.0 * sum * ex;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else if(x < xcut) {    const double x2 = x*x;    const double x4 = x2*x2;    const double x6 = x4*x2;    const double sum = 720.0 + 720.0*x + 360.0*x2 + 120.0*x2*x + 30.0*x4 + 6.0*x4*x +x6 ;    result->val = (val_infinity - 6.0 * sum * exp(-x)) / x6;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else {    result->val = (((((val_infinity/x)/x)/x)/x)/x)/x ;    result->err = GSL_DBL_EPSILON * result->val;    CHECK_UNDERFLOW(result);    return GSL_SUCCESS;  }}

int gsl_sf_psi_1_e(const double x, gsl_sf_result * result){  /* CHECK_POINTER(result) */  if(x == 0.0 || x == -1.0 || x == -2.0) {    DOMAIN_ERROR(result);  }  else if(x > 0.0)  {    return psi_n_xg0(1, x, result);  }  else if(x > -5.0)  {    /* Abramowitz + Stegun 6.4.6 */    int M = -floor(x);    double fx = x + M;    double sum = 0.0;    int m;    if(fx == 0.0)      DOMAIN_ERROR(result);    for(m = 0; m < M; ++m)      sum += 1.0/((x+m)*(x+m));    {      int stat_psi = psi_n_xg0(1, fx, result);      result->val += sum;      result->err += M * GSL_DBL_EPSILON * sum;      return stat_psi;    }  }  else  {    /* Abramowitz + Stegun 6.4.7 */    const double sin_px = sin(M_PI * x);    const double d = M_PI*M_PI/(sin_px*sin_px);    gsl_sf_result r;    int stat_psi = psi_n_xg0(1, 1.0-x, &r);    result->val = d - r.val;    result->err = r.err + 2.0*GSL_DBL_EPSILON*d;    return stat_psi;  }}

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_debye_4_e(const double x, gsl_sf_result * result){  const double val_infinity = 99.5450644937635129;  const double xcut = -GSL_LOG_DBL_MIN;  /* CHECK_POINTER(result) */  if(x < 0.0) {    DOMAIN_ERROR(result);  }  else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {    result->val = 1.0 - 2.0*x/5.0 + x*x/18.0;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else if(x <= 4.0) {    const double t = x*x/8.0 - 1.0;    gsl_sf_result c;    cheb_eval_e(&adeb4_cs, t, &c);    result->val = c.val - 2.0*x/5.0;    result->err = c.err + GSL_DBL_EPSILON * 2.0*x/5.0;    return GSL_SUCCESS;  }  else if(x < -(M_LN2 + GSL_LOG_DBL_EPSILON)) {    const int nexp = floor(xcut/x);    const double ex  = exp(-x);    double xk  = nexp * x;    double rk  = nexp;    double sum = 0.0;    int i;    for(i=nexp; i>=1; i--) {      double xk_inv = 1.0/xk;      sum *= ex;      sum += ((((24.0*xk_inv + 24.0)*xk_inv + 12.0)*xk_inv + 4.0)*xk_inv + 1.0) / rk;      rk -= 1.0;      xk -= x;    }    result->val = val_infinity/(x*x*x*x) - 4.0 * sum * ex;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else if(x < xcut) {    const double x2 = x*x;    const double x4 = x2*x2;    const double sum = 24.0 + 24.0*x + 12.0*x2 + 4.0*x2*x + x4;    result->val = (val_infinity - 4.0 * sum * exp(-x)) / x4;    result->err = GSL_DBL_EPSILON * result->val;    return GSL_SUCCESS;  }  else {    result->val = (((val_infinity/x)/x)/x)/x;    result->err = GSL_DBL_EPSILON * result->val;    CHECK_UNDERFLOW(result);    return GSL_SUCCESS;  }}

intgsl_sf_lnbeta_e(const double x, const double y, gsl_sf_result * result){  double sgn;  int status = gsl_sf_lnbeta_sgn_e(x,y,result,&sgn);  if (sgn == -1) {    DOMAIN_ERROR(result);  }  return status;}