← Back to branch summary

~ubuntu-branches/ubuntu/karmic/grace/karmic

« back to all changes in this revision

Viewing changes to cephes/zeta.c

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-19 14:19:58 UTC
  • Revision ID: james.westby@ubuntu.com-20020319141958-5gxna6vo1ek3zjml
Tags: upstream-5.1.7
Import upstream version 5.1.7

Show diffs side-by-side

added added

removed removed

Lines of Context:
cephes/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
 
 
57
/*
 
58
Cephes Math Library Release 2.0:  April, 1987
 
59
Copyright 1984, 1987 by Stephen L. Moshier
 
60
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
 
61
*/
 
62
 
 
63
#include "mconf.h"
 
64
#include "cephes.h"
 
65
 
 
66
extern double MAXNUM, MACHEP;
 
67
 
 
68
/* Expansion coefficients
 
69
 * for Euler-Maclaurin summation formula
 
70
 * (2k)! / B2k
 
71
 * where B2k are Bernoulli numbers
 
72
 */
 
73
static double A[] = {
 
74
12.0,
 
75
-720.0,
 
76
30240.0,
 
77
-1209600.0,
 
78
47900160.0,
 
79
-1.8924375803183791606e9, /*1.307674368e12/691*/
 
80
7.47242496e10,
 
81
-2.950130727918164224e12, /*1.067062284288e16/3617*/
 
82
1.1646782814350067249e14, /*5.109094217170944e18/43867*/
 
83
-4.5979787224074726105e15, /*8.028576626982912e20/174611*/
 
84
1.8152105401943546773e17, /*1.5511210043330985984e23/854513*/
 
85
-7.1661652561756670113e18 /*1.6938241367317436694528e27/236364091*/
 
86
};
 
87
/* 30 Nov 86 -- error in third coefficient fixed */
 
88
 
 
89
 
 
90
double zeta(x,q)
 
91
double x,q;
 
92
{
 
93
int i;
 
94
double a, b, k, s, t, w;
 
95
 
 
96
if( x == 1.0 )
 
97
        goto retinf;
 
98
 
 
99
if( x < 1.0 )
 
100
        {
 
101
domerr:
 
102
        mtherr( "zeta", DOMAIN );
 
103
        return(0.0);
 
104
        }
 
105
 
 
106
if( q <= 0.0 )
 
107
        {
 
108
        if(q == floor(q))
 
109
                {
 
110
                mtherr( "zeta", SING );
 
111
retinf:
 
112
                return( MAXNUM );
 
113
                }
 
114
        if( x != floor(x) )
 
115
                goto domerr; /* because q^-x not defined */
 
116
        }
 
117
 
 
118
/* Euler-Maclaurin summation formula */
 
119
/*
 
120
if( x < 25.0 )
 
121
*/
 
122
{
 
123
/* Permit negative q but continue sum until n+q > +9 .
 
124
 * This case should be handled by a reflection formula.
 
125
 * If q<0 and x is an integer, there is a relation to
 
126
 * the polygamma function.
 
127
 */
 
128
s = pow( q, -x );
 
129
a = q;
 
130
i = 0;
 
131
b = 0.0;
 
132
while( (i < 9) || (a <= 9.0) )
 
133
        {
 
134
        i += 1;
 
135
        a += 1.0;
 
136
        b = pow( a, -x );
 
137
        s += b;
 
138
        if( fabs(b/s) < MACHEP )
 
139
                goto done;
 
140
        }
 
141
 
 
142
w = a;
 
143
s += b*w/(x-1.0);
 
144
s -= 0.5 * b;
 
145
a = 1.0;
 
146
k = 0.0;
 
147
for( i=0; i<12; i++ )
 
148
        {
 
149
        a *= x + k;
 
150
        b /= w;
 
151
        t = a*b/A[i];
 
152
        s = s + t;
 
153
        t = fabs(t/s);
 
154
        if( t < MACHEP )
 
155
                goto done;
 
156
        k += 1.0;
 
157
        a *= x + k;
 
158
        b /= w;
 
159
        k += 1.0;
 
160
        }
 
161
done:
 
162
return(s);
 
163
}
 
164
 
 
165
 
 
166
 
 
167
/* Basic sum of inverse powers */
 
168
/*
 
169
pseres:
 
170
 
 
171
s = pow( q, -x );
 
172
a = q;
 
173
do
 
174
        {
 
175
        a += 2.0;
 
176
        b = pow( a, -x );
 
177
        s += b;
 
178
        }
 
179
while( b/s > MACHEP );
 
180
 
 
181
b = pow( 2.0, -x );
 
182
s = (s + b)/(1.0-b);
 
183
return(s);
 
184
*/
 
185
}