/* zeta.c
*
* Riemann zeta function of two arguments
*
*
*
* SYNOPSIS:
*
* double x, q, y, zeta();
*
* y = zeta( x, q );
*
*
*
* DESCRIPTION:
*
*
*
* inf.
* - -x
* zeta(x,q) = > (k+q)
* -
* k=0
*
* where x > 1 and q is not a negative integer or zero.
* The Euler-Maclaurin summation formula is used to obtain
* the expansion
*
* n
* - -x
* zeta(x,q) = > (k+q)
* -
* k=1
*
* 1-x inf. B x(x+1)...(x+2j)
* (n+q) 1 - 2j
* + --------- - ------- + > --------------------
* x-1 x - x+2j+1
* 2(n+q) j=1 (2j)! (n+q)
*
* where the B2j are Bernoulli numbers. Note that (see zetac.c)
* zeta(x,1) = zetac(x) + 1.
*
*
*
* ACCURACY:
*
*
*
* REFERENCE:
*
* Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
* Series, and Products, p. 1073; Academic Press, 1980.
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1987, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double md_fabs ( double );
extern double md_pow ( double, double );
extern double md_floor ( double );
#else
double md_fabs(), md_pow(), md_floor();
#endif
extern double MAXNUM, MACHEP;
/* Expansion coefficients
* for Euler-Maclaurin summation formula
* (2k)! / B2k
* where B2k are Bernoulli numbers
*/
static double A[] = {
12.0,
-720.0,
30240.0,
-1209600.0,
47900160.0,
-1.8924375803183791606e9, /*1.307674368e12/691*/
7.47242496e10,
-2.950130727918164224e12, /*1.067062284288e16/3617*/
1.1646782814350067249e14, /*5.109094217170944e18/43867*/
-4.5979787224074726105e15, /*8.028576626982912e20/174611*/
1.8152105401943546773e17, /*1.5511210043330985984e23/854513*/
-7.1661652561756670113e18 /*1.6938241367317436694528e27/236364091*/
};
/* 30 Nov 86 -- error in third coefficient fixed */
double zeta(x,q)
double x,q;
{
int i;
double a, b, k, s, t, w;
if( x == 1.0 )
goto retinf;
if( x < 1.0 )
{
domerr:
mtherr( "zeta", DOMAIN );
return(0.0);
}
if( q <= 0.0 )
{
if(q == md_floor(q))
{
mtherr( "zeta", SING );
retinf:
return( MAXNUM );
}
if( x != md_floor(x) )
goto domerr; /* because q^-x not defined */
}
/* Euler-Maclaurin summation formula */
/*
if( x < 25.0 )
*/
{
/* Permit negative q but continue sum until n+q > +9 .
* This case should be handled by a reflection formula.
* If q<0 and x is an integer, there is a relation to
* the polygamma function.
*/
s = md_pow( q, -x );
a = q;
i = 0;
b = 0.0;
while( (i < 9) || (a <= 9.0) )
{
i += 1;
a += 1.0;
b = md_pow( a, -x );
s += b;
if( md_fabs(b/s) < MACHEP )
goto done;
}
w = a;
s += b*w/(x-1.0);
s -= 0.5 * b;
a = 1.0;
k = 0.0;
for( i=0; i<12; i++ )
{
a *= x + k;
b /= w;
t = a*b/A[i];
s = s + t;
t = md_fabs(t/s);
if( t < MACHEP )
goto done;
k += 1.0;
a *= x + k;
b /= w;
k += 1.0;
}
done:
return(s);
}
/* Basic sum of inverse powers */
/*
pseres:
s = md_pow( q, -x );
a = q;
do
{
a += 2.0;
b = md_pow( a, -x );
s += b;
}
while( b/s > MACHEP );
b = md_pow( 2.0, -x );
s = (s + b)/(1.0-b);
return(s);
*/
}