zeta.c

开放gsl矩阵运算

0.99999999999997157829090808339,
0.99999999999998578914539762720,
0.99999999999999289457268000875,
0.99999999999999644728633373609,
0.99999999999999822364316477861,
0.99999999999999911182158169283,
0.99999999999999955591079061426,
0.99999999999999977795539522974,
0.99999999999999988897769758908,
0.99999999999999994448884878594,
0.99999999999999997224442439010,
0.99999999999999998612221219410,
0.99999999999999999306110609673,
0.99999999999999999653055304826,
0.99999999999999999826527652409,
0.99999999999999999913263826204,
0.99999999999999999956631913101,
0.99999999999999999978315956551,
0.99999999999999999989157978275,
0.99999999999999999994578989138,
0.99999999999999999997289494569,
0.99999999999999999998644747284,
0.99999999999999999999322373642,
0.99999999999999999999661186821,
0.99999999999999999999830593411,
0.99999999999999999999915296705,
0.99999999999999999999957648353,
0.99999999999999999999978824176,
0.99999999999999999999989412088,
0.99999999999999999999994706044,
0.99999999999999999999997353022,
0.99999999999999999999998676511,
0.99999999999999999999999338256,
0.99999999999999999999999669128,
0.99999999999999999999999834564,
0.99999999999999999999999917282,
0.99999999999999999999999958641,
0.99999999999999999999999979320,
0.99999999999999999999999989660,
0.99999999999999999999999994830,
0.99999999999999999999999997415,
0.99999999999999999999999998708,
0.99999999999999999999999999354,
0.99999999999999999999999999677,
0.99999999999999999999999999838,
0.99999999999999999999999999919,
0.99999999999999999999999999960,
0.99999999999999999999999999980,
0.99999999999999999999999999990,
0.99999999999999999999999999995,
0.99999999999999999999999999997,
0.99999999999999999999999999999,
0.99999999999999999999999999999,
1.00000000000000000000000000000,
1.00000000000000000000000000000,
1.00000000000000000000000000000,
};


#define ETA_NEG_TABLE_NMAX  99
#define ETA_NEG_TABLE_SIZE  50
static double eta_neg_int_table[ETA_NEG_TABLE_SIZE] = {
 0.25000000000000000000000000000,   /* eta(-1) */
-0.12500000000000000000000000000,   /* eta(-3) */
 0.25000000000000000000000000000,   /* ...	*/
-1.06250000000000000000000000000,
 7.75000000000000000000000000000,
-86.3750000000000000000000000000,
 1365.25000000000000000000000000,
-29049.0312500000000000000000000,
 800572.750000000000000000000000,
-2.7741322625000000000000000000e+7,
 1.1805291302500000000000000000e+9,
-6.0523980051687500000000000000e+10,
 3.6794167785377500000000000000e+12,
-2.6170760990658387500000000000e+14,
 2.1531418140800295250000000000e+16,
-2.0288775575173015930156250000e+18,
 2.1708009902623770590275000000e+20,
-2.6173826968455814932120125000e+22,
 3.5324148876863877826668602500e+24,
-5.3042033406864906641493838981e+26,
 8.8138218364311576767253114668e+28,
-1.6128065107490778547354654864e+31,
 3.2355470001722734208527794569e+33,
-7.0876727476537493198506645215e+35,
 1.6890450341293965779175629389e+38,
-4.3639690731216831157655651358e+40,
 1.2185998827061261322605065672e+43,
-3.6670584803153006180101262324e+45,
 1.1859898526302099104271449748e+48,
-4.1120769493584015047981746438e+50,
 1.5249042436787620309090168687e+53,
-6.0349693196941307074572991901e+55,
 2.5437161764210695823197691519e+58,
-1.1396923802632287851130360170e+61,
 5.4180861064753979196802726455e+63,
-2.7283654799994373847287197104e+66,
 1.4529750514918543238511171663e+69,
-8.1705519371067450079777183386e+71,
 4.8445781606678367790247757259e+74,
-3.0246694206649519336179448018e+77,
 1.9858807961690493054169047970e+80,
-1.3694474620720086994386818232e+83,
 9.9070382984295807826303785989e+85,
-7.5103780796592645925968460677e+88,
 5.9598418264260880840077992227e+91,
-4.9455988887500020399263196307e+94,
 4.2873596927020241277675775935e+97,
-3.8791952037716162900707994047e+100,
 3.6600317773156342245401829308e+103,
-3.5978775704117283875784869570e+106    /* eta(-99)  */
};


/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/


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


int gsl_sf_zeta_e(const double s, gsl_sf_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_sf_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_sf_result g;
      const int stat_g = gsl_sf_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_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 = 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_sf_eta_int_e(int n, gsl_sf_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_sf_result z;
      gsl_sf_result p;
      int stat_z = gsl_sf_zeta_int_e(n, &z);
      int stat_p = gsl_sf_exp_e((1.0-n)*M_LN2, &p);
      int stat_m = gsl_sf_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_sf_eta_e(const double s, gsl_sf_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_sf_result z;
    gsl_sf_result p;
    int stat_z = gsl_sf_zeta_e(s, &z);
    int stat_p = gsl_sf_exp_e((1.0-s)*M_LN2, &p);
    int stat_m = gsl_sf_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 *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_zeta(const double s)
{
  EVAL_RESULT(gsl_sf_zeta_e(s, &result));
}

double gsl_sf_hzeta(const double s, const double a)
{
  EVAL_RESULT(gsl_sf_hzeta_e(s, a, &result));
}

double gsl_sf_zeta_int(const int s)
{
  EVAL_RESULT(gsl_sf_zeta_int_e(s, &result));
}

double gsl_sf_eta_int(const int s)
{
  EVAL_RESULT(gsl_sf_eta_int_e(s, &result));
}

double gsl_sf_eta(const double s)
{
  EVAL_RESULT(gsl_sf_eta_e(s, &result));
}