diff options
Diffstat (limited to 'libm/double/zeta.c')
-rw-r--r-- | libm/double/zeta.c | 189 |
1 files changed, 0 insertions, 189 deletions
diff --git a/libm/double/zeta.c b/libm/double/zeta.c deleted file mode 100644 index a49c619d5..000000000 --- a/libm/double/zeta.c +++ /dev/null @@ -1,189 +0,0 @@ -/* 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); -*/ -} |