/[gsl]/gsl/cdf/zeta.c
ViewVC logotype

Diff of /gsl/cdf/zeta.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1.2.1 by jstover, Tue Dec 31 16:16:47 2002 UTC revision 1.1.2.2 by jstover, Fri Jan 3 04:09:25 2003 UTC
# Line 1  Line 1 
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   *   *
# Line 14  Line 14 
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 */
# Line 109  static cheb_series zeta_xgt1_cs = { Line 109  static cheb_series zeta_xgt1_cs = {
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] = {
# Line 506  static double eta_neg_int_table[ETA_NEG_ Line 452  static double eta_neg_int_table[ETA_NEG_
452  -3.5978775704117283875784869570e+106    /* eta(-99)  */  -3.5978775704117283875784869570e+106    /* eta(-99)  */
453  };  };
454    
   
455  /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/  /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
456    
457    
# Line 571  int gsl_cdf_hzeta_e(const double s, cons Line 516  int gsl_cdf_hzeta_e(const double s, cons
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));  
 }  

Legend:
Removed from v.1.1.2.1  
changed lines
  Added in v.1.1.2.2

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26