zeta.c

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY without ev

/* specfunc/zeta.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_elementary.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf_zeta.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

#define LogTwoPi_  1.8378770664093454835606594728111235279723


/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* chebyshev fit for (s(t)-1)Zeta[s(t)]
 * s(t)= (t+1)/2
 * -1 <= t <= 1
 */
static double zeta_xlt1_data[14] = {
  1.48018677156931561235192914649,
  0.25012062539889426471999938167,
  0.00991137502135360774243761467,
 -0.00012084759656676410329833091,
 -4.7585866367662556504652535281e-06,
  2.2229946694466391855561441361e-07,
 -2.2237496498030257121309056582e-09,
 -1.0173226513229028319420799028e-10,
  4.3756643450424558284466248449e-12,
 -6.2229632593100551465504090814e-14,
 -6.6116201003272207115277520305e-16,
  4.9477279533373912324518463830e-17,
 -1.0429819093456189719660003522e-18,
  6.9925216166580021051464412040e-21,
};
static cheb_series zeta_xlt1_cs = {
  zeta_xlt1_data,
  13,
  -1, 1,
  8
};

/* chebyshev fit for (s(t)-1)Zeta[s(t)]
 * s(t)= (19t+21)/2
 * -1 <= t <= 1
 */
static double zeta_xgt1_data[30] = {
  19.3918515726724119415911269006,
   9.1525329692510756181581271500,
   0.2427897658867379985365270155,
  -0.1339000688262027338316641329,
   0.0577827064065028595578410202,
  -0.0187625983754002298566409700,
   0.0039403014258320354840823803,
  -0.0000581508273158127963598882,
  -0.0003756148907214820704594549,
   0.0001892530548109214349092999,
  -0.0000549032199695513496115090,
   8.7086484008939038610413331863e-6,
   6.4609477924811889068410083425e-7,
  -9.6749773915059089205835337136e-7,
   3.6585400766767257736982342461e-7,
  -8.4592516427275164351876072573e-8,
   9.9956786144497936572288988883e-9,
   1.4260036420951118112457144842e-9,
  -1.1761968823382879195380320948e-9,
   3.7114575899785204664648987295e-10,
  -7.4756855194210961661210215325e-11,
   7.8536934209183700456512982968e-12,
   9.9827182259685539619810406271e-13,
  -7.5276687030192221587850302453e-13,
   2.1955026393964279988917878654e-13,
  -4.1934859852834647427576319246e-14,
   4.6341149635933550715779074274e-15,
   2.3742488509048340106830309402e-16,
  -2.7276516388124786119323824391e-16,
   7.8473570134636044722154797225e-17
};
static cheb_series zeta_xgt1_cs = {
  zeta_xgt1_data,
  29,
  -1, 1,
  17
};


/* chebyshev fit for Ln[Zeta[s(t)] - 1 - 2^(-s(t))]
 * s(t)= 10 + 5t
 * -1 <= t <= 1; 5 <= s <= 15
 */
static double zetam1_inter_data[24] = {
  -21.7509435653088483422022339374,
  -5.63036877698121782876372020472,
   0.0528041358684229425504861579635,
  -0.0156381809179670789342700883562,
   0.00408218474372355881195080781927,
  -0.0010264867349474874045036628282,
   0.000260469880409886900143834962387,
  -0.0000676175847209968878098566819447,
   0.0000179284472587833525426660171124,
  -4.83238651318556188834107605116e-6,
   1.31913788964999288471371329447e-6,
  -3.63760500656329972578222188542e-7,
   1.01146847513194744989748396574e-7,
  -2.83215225141806501619105289509e-8,
   7.97733710252021423361012829496e-9,
  -2.25850168553956886676250696891e-9,
   6.42269392950164306086395744145e-10,
  -1.83363861846127284505060843614e-10,
   5.25309763895283179960368072104e-11,
  -1.50958687042589821074710575446e-11,
   4.34997545516049244697776942981e-12,
  -1.25597782748190416118082322061e-12,
   3.61280740072222650030134104162e-13,
  -9.66437239205745207188920348801e-14
}; 
static cheb_series zetam1_inter_cs = {
  zetam1_inter_data,
  22,
  -1, 1,
  12
};



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


inline
static int
riemann_zeta1ms_slt0(double s, gsl_sf_result * result)
{
  if(s > -19.0) {
    double x = (-19 - 2.0*s)/19.0;
    gsl_sf_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;
  }
}


/* works for 5 < s < 15*/
static int
riemann_zeta_minus_1_intermediate_s(double s, gsl_sf_result * result)
{
  double t = (s - 10.0)/5.0;
  gsl_sf_result c;
  cheb_eval_e(&zetam1_inter_cs, t, &c);
  result->val = exp(c.val) + pow(2.0, -s);
  result->err = (c.err + 2.0*GSL_DBL_EPSILON)*result->val;
  return GSL_SUCCESS;
}


/* assumes s is large and positive
 * write: zeta(s) - 1 = zeta(s) * (1 - 1/zeta(s))
 * and expand a few terms of the product formula to evaluate 1 - 1/zeta(s)
 *
 * works well for s > 15
 */
static int
riemann_zeta_minus1_large_s(double s, gsl_sf_result * result)
{
  double a = pow( 2.0,-s);
  double b = pow( 3.0,-s);
  double c = pow( 5.0,-s);
  double d = pow( 7.0,-s);
  double e = pow(11.0,-s);
  double f = pow(13.0,-s);
  double t1 = a + b + c + d + e + f;
  double t2 = a*(b+c+d+e+f) + b*(c+d+e+f) + c*(d+e+f) + d*(e+f) + e*f;
  /*
  double t3 = a*(b*(c+d+e+f) + c*(d+e+f) + d*(e+f) + e*f) +
              b*(c*(d+e+f) + d*(e+f) + e*f) +
              c*(d*(e+f) + e*f) +
              d*e*f;
  double t4 = a*(b*(c*(d + e + f) + d*(e + f) + e*f) + c*(d*(e+f) + e*f) + d*e*f) +
              b*(c*(d*(e+f) + e*f) + d*e*f) +
              c*d*e*f;
  double t5 = b*c*d*e*f + a*c*d*e*f+ a*b*d*e*f+ a*b*c*e*f+ a*b*c*d*f+ a*b*c*d*e;
  double t6 = a*b*c*d*e*f;
  */
  double numt = t1 - t2 /* + t3 - t4 + t5 - t6 */;
  double zeta = 1.0/((1.0-a)*(1.0-b)*(1.0-c)*(1.0-d)*(1.0-e)*(1.0-f));
  result->val = numt*zeta;
  result->err = (15.0/s + 1.0) * 6.0*GSL_DBL_EPSILON*result->val;
  return GSL_SUCCESS;
}


#if 0
/* zeta(n) */
#define ZETA_POS_TABLE_NMAX   100
static double zeta_pos_int_table_OLD[ZETA_POS_TABLE_NMAX+1] = {
 -0.50000000000000000000000000000,       /* zeta(0) */
  0.0 /* FIXME: DirectedInfinity() */,   /* zeta(1) */
  1.64493406684822643647241516665,       /* ...     */
  1.20205690315959428539973816151,
  1.08232323371113819151600369654,
  1.03692775514336992633136548646,
  1.01734306198444913971451792979,
  1.00834927738192282683979754985,
  1.00407735619794433937868523851,
  1.00200839282608221441785276923,
  1.00099457512781808533714595890,
  1.00049418860411946455870228253,
  1.00024608655330804829863799805,
  1.00012271334757848914675183653,
  1.00006124813505870482925854511,
  1.00003058823630702049355172851,
  1.00001528225940865187173257149,
  1.00000763719763789976227360029,
  1.00000381729326499983985646164,
  1.00000190821271655393892565696,
  1.00000095396203387279611315204,
  1.00000047693298678780646311672,
  1.00000023845050272773299000365,
  1.00000011921992596531107306779,
  1.00000005960818905125947961244,
  1.00000002980350351465228018606,
  1.00000001490155482836504123466,
  1.00000000745071178983542949198,
  1.00000000372533402478845705482,
  1.00000000186265972351304900640,
  1.00000000093132743241966818287,
  1.00000000046566290650337840730,
  1.00000000023283118336765054920,
  1.00000000011641550172700519776,
  1.00000000005820772087902700889,
  1.00000000002910385044497099687,
  1.00000000001455192189104198424,
  1.00000000000727595983505748101,
  1.00000000000363797954737865119,
  1.00000000000181898965030706595,
  1.00000000000090949478402638893,
  1.00000000000045474737830421540,
  1.00000000000022737368458246525,
  1.00000000000011368684076802278,
  1.00000000000005684341987627586,
  1.00000000000002842170976889302,
  1.00000000000001421085482803161,
  1.00000000000000710542739521085,
  1.00000000000000355271369133711,
  1.00000000000000177635684357912,
  1.00000000000000088817842109308,
  1.00000000000000044408921031438,
  1.00000000000000022204460507980,
  1.00000000000000011102230251411,
  1.00000000000000005551115124845,
  1.00000000000000002775557562136,
  1.00000000000000001387778780973,
  1.00000000000000000693889390454,
  1.00000000000000000346944695217,
  1.00000000000000000173472347605,
  1.00000000000000000086736173801,
  1.00000000000000000043368086900,
  1.00000000000000000021684043450,
  1.00000000000000000010842021725,
  1.00000000000000000005421010862,
  1.00000000000000000002710505431,
  1.00000000000000000001355252716,
  1.00000000000000000000677626358,
  1.00000000000000000000338813179,
  1.00000000000000000000169406589,
  1.00000000000000000000084703295,
  1.00000000000000000000042351647,
  1.00000000000000000000021175824,
  1.00000000000000000000010587912,
  1.00000000000000000000005293956,
  1.00000000000000000000002646978,
  1.00000000000000000000001323489,
  1.00000000000000000000000661744,
  1.00000000000000000000000330872,
  1.00000000000000000000000165436,
  1.00000000000000000000000082718,
  1.00000000000000000000000041359,
  1.00000000000000000000000020680,
  1.00000000000000000000000010340,
  1.00000000000000000000000005170,
  1.00000000000000000000000002585,
  1.00000000000000000000000001292,
  1.00000000000000000000000000646,
  1.00000000000000000000000000323,
  1.00000000000000000000000000162,
  1.00000000000000000000000000081,