123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189 |
- /* 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 <math.h>
- #ifdef ANSIPROT
- extern double fabs ( double );
- extern double pow ( double, double );
- extern double floor ( double );
- #else
- double fabs(), pow(), 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 == floor(q))
- {
- mtherr( "zeta", SING );
- retinf:
- return( MAXNUM );
- }
- if( x != 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 = pow( q, -x );
- a = q;
- i = 0;
- b = 0.0;
- while( (i < 9) || (a <= 9.0) )
- {
- i += 1;
- a += 1.0;
- b = pow( a, -x );
- s += b;
- if( 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 = 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 = pow( q, -x );
- a = q;
- do
- {
- a += 2.0;
- b = pow( a, -x );
- s += b;
- }
- while( b/s > MACHEP );
- b = pow( 2.0, -x );
- s = (s + b)/(1.0-b);
- return(s);
- */
- }
|