libm/double/zeta.c

   1 /*                                                      zeta.c
   2  *
   3  *      Riemann zeta function of two arguments
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, q, y, zeta();
  10  *
  11  * y = zeta( x, q );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  *
  18  *
  19  *                 inf.
  20  *                  -        -x
  21  *   zeta(x,q)  =   >   (k+q)
  22  *                  -
  23  *                 k=0
  24  *
  25  * where x > 1 and q is not a negative integer or zero.
  26  * The Euler-Maclaurin summation formula is used to obtain
  27  * the expansion
  28  *
  29  *                n
  30  *                -       -x
  31  * zeta(x,q)  =   >  (k+q)
  32  *                -
  33  *               k=1
  34  *
  35  *           1-x                 inf.  B   x(x+1)...(x+2j)
  36  *      (n+q)           1         -     2j
  37  *  +  ---------  -  -------  +   >    --------------------
  38  *        x-1              x      -                   x+2j+1
  39  *                   2(n+q)      j=1       (2j)! (n+q)
  40  *
  41  * where the B2j are Bernoulli numbers.  Note that (see zetac.c)
  42  * zeta(x,1) = zetac(x) + 1.
  43  *
  44  *
  45  *
  46  * ACCURACY:
  47  *
  48  *
  49  *
  50  * REFERENCE:
  51  *
  52  * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
  53  * Series, and Products, p. 1073; Academic Press, 1980.
  54  *
  55  */
  56 \f
  57 /*
  58 Cephes Math Library Release 2.8:  June, 2000
  59 Copyright 1984, 1987, 2000 by Stephen L. Moshier
  60 */
  61
  62 #include <math.h>
  63 #ifdef ANSIPROT
  64 extern double fabs ( double );
  65 extern double pow ( double, double );
  66 extern double floor ( double );
  67 #else
  68 double fabs(), pow(), floor();
  69 #endif
  70 extern double MAXNUM, MACHEP;
  71
  72 /* Expansion coefficients
  73  * for Euler-Maclaurin summation formula
  74  * (2k)! / B2k
  75  * where B2k are Bernoulli numbers
  76  */
  77 static double A[] = {
  78 12.0,
  79 -720.0,
  80 30240.0,
  81 -1209600.0,
  82 47900160.0,
  83 -1.8924375803183791606e9, /*1.307674368e12/691*/
  84 7.47242496e10,
  85 -2.950130727918164224e12, /*1.067062284288e16/3617*/
  86 1.1646782814350067249e14, /*5.109094217170944e18/43867*/
  87 -4.5979787224074726105e15, /*8.028576626982912e20/174611*/
  88 1.8152105401943546773e17, /*1.5511210043330985984e23/854513*/
  89 -7.1661652561756670113e18 /*1.6938241367317436694528e27/236364091*/
  90 };
  91 /* 30 Nov 86 -- error in third coefficient fixed */
  92
  93
  94 double zeta(x,q)
  95 double x,q;
  96 {
  97 int i;
  98 double a, b, k, s, t, w;
  99
 100 if( x == 1.0 )
 101         goto retinf;
 102
 103 if( x < 1.0 )
 104         {
 105 domerr:
 106         mtherr( "zeta", DOMAIN );
 107         return(0.0);
 108         }
 109
 110 if( q <= 0.0 )
 111         {
 112         if(q == floor(q))
 113                 {
 114                 mtherr( "zeta", SING );
 115 retinf:
 116                 return( MAXNUM );
 117                 }
 118         if( x != floor(x) )
 119                 goto domerr; /* because q^-x not defined */
 120         }
 121
 122 /* Euler-Maclaurin summation formula */
 123 /*
 124 if( x < 25.0 )
 125 */
 126 {
 127 /* Permit negative q but continue sum until n+q > +9 .
 128  * This case should be handled by a reflection formula.
 129  * If q<0 and x is an integer, there is a relation to
 130  * the polygamma function.
 131  */
 132 s = pow( q, -x );
 133 a = q;
 134 i = 0;
 135 b = 0.0;
 136 while( (i < 9) || (a <= 9.0) )
 137         {
 138         i += 1;
 139         a += 1.0;
 140         b = pow( a, -x );
 141         s += b;
 142         if( fabs(b/s) < MACHEP )
 143                 goto done;
 144         }
 145
 146 w = a;
 147 s += b*w/(x-1.0);
 148 s -= 0.5 * b;
 149 a = 1.0;
 150 k = 0.0;
 151 for( i=0; i<12; i++ )
 152         {
 153         a *= x + k;
 154         b /= w;
 155         t = a*b/A[i];
 156         s = s + t;
 157         t = fabs(t/s);
 158         if( t < MACHEP )
 159                 goto done;
 160         k += 1.0;
 161         a *= x + k;
 162         b /= w;
 163         k += 1.0;
 164         }
 165 done:
 166 return(s);
 167 }
 168
 169
 170
 171 /* Basic sum of inverse powers */
 172 /*
 173 pseres:
 174
 175 s = pow( q, -x );
 176 a = q;
 177 do
 178         {
 179         a += 2.0;
 180         b = pow( a, -x );
 181         s += b;
 182         }
 183 while( b/s > MACHEP );
 184
 185 b = pow( 2.0, -x );
 186 s = (s + b)/(1.0-b);
 187 return(s);
 188 */
 189 }