1 |
/* specfunc/zeta.c |
/* cdf/zeta.c |
2 |
* |
* |
3 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman |
4 |
* |
* |
14 |
* |
* |
15 |
* You should have received a copy of the GNU General Public License |
* You should have received a copy of the GNU General Public License |
16 |
* along with this program; if not, write to the Free Software |
* along with this program; if not, write to the Free Software |
17 |
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. |
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA. |
18 |
*/ |
*/ |
19 |
|
|
20 |
/* Author: G. Jungman */ |
/* Author: G. Jungman */ |
109 |
}; |
}; |
110 |
|
|
111 |
|
|
|
/* assumes s >= 0 and s != 1.0 */ |
|
|
inline |
|
|
static int |
|
|
riemann_zeta_sgt0(double s, gsl_cdf_result * result) |
|
|
{ |
|
|
if(s < 1.0) { |
|
|
gsl_cdf_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_cdf_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; |
|
|
} |
|
|
} |
|
|
|
|
|
inline |
|
|
static int |
|
|
riemann_zeta1m_slt0(double s, gsl_cdf_result * result) |
|
|
{ |
|
|
if(s > -19.0) { |
|
|
double x = (-19 - 2.0*s)/19.0; |
|
|
gsl_cdf_result c; |
|
|
cheb_eval_e(&zeta_xgt1_cs, x, &c); |
|
|
result->val = c.val / (-s); |
|
|
result->err = c.err / (-s) + GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else { |
|
|
double f2 = 1.0 - pow(2.0,-(1.0-s)); |
|
|
double f3 = 1.0 - pow(3.0,-(1.0-s)); |
|
|
double f5 = 1.0 - pow(5.0,-(1.0-s)); |
|
|
double f7 = 1.0 - pow(7.0,-(1.0-s)); |
|
|
result->val = 1.0/(f2*f3*f5*f7); |
|
|
result->err = 3.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
} |
|
|
|
|
112 |
/* zeta(n) */ |
/* zeta(n) */ |
113 |
#define ZETA_POS_TABLE_NMAX 100 |
#define ZETA_POS_TABLE_NMAX 100 |
114 |
static double zeta_pos_int_table[ZETA_POS_TABLE_NMAX+1] = { |
static double zeta_pos_int_table[ZETA_POS_TABLE_NMAX+1] = { |
452 |
-3.5978775704117283875784869570e+106 /* eta(-99) */ |
-3.5978775704117283875784869570e+106 /* eta(-99) */ |
453 |
}; |
}; |
454 |
|
|
|
|
|
455 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
456 |
|
|
457 |
|
|
516 |
} |
} |
517 |
} |
} |
518 |
|
|
|
|
|
|
int gsl_cdf_zeta_e(const double s, gsl_cdf_result * result) |
|
|
{ |
|
|
/* CHECK_POINTER(result) */ |
|
|
|
|
|
if(s == 1.0) { |
|
|
DOMAIN_ERROR(result); |
|
|
} |
|
|
else if(s >= 0.0) { |
|
|
return riemann_zeta_sgt0(s, result); |
|
|
} |
|
|
else { |
|
|
/* reflection formula, [Abramowitz+Stegun, 23.2.5] */ |
|
|
|
|
|
gsl_cdf_result zeta_one_minus_s; |
|
|
const int stat_zoms = riemann_zeta1m_slt0(s, &zeta_one_minus_s); |
|
|
const double sin_term = sin(0.5*M_PI*s)/M_PI; |
|
|
|
|
|
if(sin_term == 0.0) { |
|
|
result->val = 0.0; |
|
|
result->err = 0.0; |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else if(s > -170) { |
|
|
/* We have to be careful about losing digits |
|
|
* in calculating pow(2 Pi, s). The gamma |
|
|
* function is fine because we were careful |
|
|
* with that implementation. |
|
|
* We keep an array of (2 Pi)^(10 n). |
|
|
*/ |
|
|
const double twopi_pow[18] = { 1.0, |
|
|
9.589560061550901348e+007, |
|
|
9.195966217409212684e+015, |
|
|
8.818527036583869903e+023, |
|
|
8.456579467173150313e+031, |
|
|
8.109487671573504384e+039, |
|
|
7.776641909496069036e+047, |
|
|
7.457457466828644277e+055, |
|
|
7.151373628461452286e+063, |
|
|
6.857852693272229709e+071, |
|
|
6.576379029540265771e+079, |
|
|
6.306458169130020789e+087, |
|
|
6.047615938853066678e+095, |
|
|
5.799397627482402614e+103, |
|
|
5.561367186955830005e+111, |
|
|
5.333106466365131227e+119, |
|
|
5.114214477385391780e+127, |
|
|
4.904306689854036836e+135 |
|
|
}; |
|
|
const int n = floor((-s)/10.0); |
|
|
const double fs = s + 10.0*n; |
|
|
const double p = pow(2.0*M_PI, fs) / twopi_pow[n]; |
|
|
|
|
|
gsl_cdf_result g; |
|
|
const int stat_g = gsl_cdf_gamma_e(1.0-s, &g); |
|
|
result->val = p * g.val * sin_term * zeta_one_minus_s.val; |
|
|
result->err = fabs(p * g.val * sin_term) * zeta_one_minus_s.err; |
|
|
result->err += fabs(p * sin_term * zeta_one_minus_s.val) * g.err; |
|
|
result->err += GSL_DBL_EPSILON * (fabs(s)+2.0) * fabs(result->val); |
|
|
return GSL_ERROR_SELECT_2(stat_g, stat_zoms); |
|
|
} |
|
|
else { |
|
|
/* The actual zeta function may or may not |
|
|
* overflow here. But we have no easy way |
|
|
* to calculate it when the prefactor(s) |
|
|
* overflow. Trying to use log's and exp |
|
|
* is no good because we loose a couple |
|
|
* digits to the exp error amplification. |
|
|
* When we gather a little more patience, |
|
|
* we can implement something here. Until |
|
|
* then just give up. |
|
|
*/ |
|
|
OVERFLOW_ERROR(result); |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
int gsl_cdf_zeta_int_e(const int n, gsl_cdf_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_cdf_zeta_e((double)n, result); |
|
|
} |
|
|
} |
|
|
else if(n == 1){ |
|
|
DOMAIN_ERROR(result); |
|
|
} |
|
|
else if(n <= ZETA_POS_TABLE_NMAX){ |
|
|
result->val = zeta_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_cdf_eta_int_e(int n, gsl_cdf_result * result) |
|
|
{ |
|
|
if(n > ETA_POS_TABLE_NMAX) { |
|
|
result->val = 1.0; |
|
|
result->err = GSL_DBL_EPSILON; |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else if(n >= 0) { |
|
|
result->val = eta_pos_int_table[n]; |
|
|
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else { |
|
|
/* n < 0 */ |
|
|
|
|
|
if(!GSL_IS_ODD(n)) { |
|
|
/* exactly zero at even negative integers */ |
|
|
result->val = 0.0; |
|
|
result->err = 0.0; |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else if(n > -ETA_NEG_TABLE_NMAX) { |
|
|
result->val = eta_neg_int_table[-(n+1)/2]; |
|
|
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else { |
|
|
gsl_cdf_result z; |
|
|
gsl_cdf_result p; |
|
|
int stat_z = gsl_cdf_zeta_int_e(n, &z); |
|
|
int stat_p = gsl_cdf_exp_e((1.0-n)*M_LN2, &p); |
|
|
int stat_m = gsl_cdf_multiply_e(-p.val, z.val, result); |
|
|
result->err = fabs(p.err * (M_LN2*(1.0-n)) * z.val) + z.err * fabs(p.val); |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_ERROR_SELECT_3(stat_m, stat_p, stat_z); |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
int gsl_cdf_eta_e(const double s, gsl_cdf_result * result) |
|
|
{ |
|
|
/* CHECK_POINTER(result) */ |
|
|
|
|
|
if(s > 100.0) { |
|
|
result->val = 1.0; |
|
|
result->err = GSL_DBL_EPSILON; |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else if(fabs(s-1.0) < 10.0*GSL_ROOT5_DBL_EPSILON) { |
|
|
double del = s-1.0; |
|
|
double c0 = M_LN2; |
|
|
double c1 = M_LN2 * (M_EULER - 0.5*M_LN2); |
|
|
double c2 = -0.0326862962794492996; |
|
|
double c3 = 0.0015689917054155150; |
|
|
double c4 = 0.00074987242112047532; |
|
|
result->val = c0 + del * (c1 + del * (c2 + del * (c3 + del * c4))); |
|
|
result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_SUCCESS; |
|
|
} |
|
|
else { |
|
|
gsl_cdf_result z; |
|
|
gsl_cdf_result p; |
|
|
int stat_z = gsl_cdf_zeta_e(s, &z); |
|
|
int stat_p = gsl_cdf_exp_e((1.0-s)*M_LN2, &p); |
|
|
int stat_m = gsl_cdf_multiply_e(1.0-p.val, z.val, result); |
|
|
result->err = fabs(p.err * (M_LN2*(1.0-s)) * z.val) + z.err * fabs(p.val); |
|
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
|
|
return GSL_ERROR_SELECT_3(stat_m, stat_p, stat_z); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
|
|
|
|
519 |
#include "eval.h" |
#include "eval.h" |
520 |
|
|
|
double gsl_cdf_zeta(const double s) |
|
|
{ |
|
|
EVAL_RESULT(gsl_cdf_zeta_e(s, &result)); |
|
|
} |
|
|
|
|
|
double gsl_cdf_hzeta(const double s, const double a) |
|
|
{ |
|
|
EVAL_RESULT(gsl_cdf_hzeta_e(s, a, &result)); |
|
|
} |
|
|
|
|
|
double gsl_cdf_zeta_int(const int s) |
|
|
{ |
|
|
EVAL_RESULT(gsl_cdf_zeta_int_e(s, &result)); |
|
|
} |
|
|
|
|
|
double gsl_cdf_eta_int(const int s) |
|
|
{ |
|
|
EVAL_RESULT(gsl_cdf_eta_int_e(s, &result)); |
|
|
} |
|
|
|
|
|
double gsl_cdf_eta(const double s) |
|
|
{ |
|
|
EVAL_RESULT(gsl_cdf_eta_e(s, &result)); |
|
|
} |
|