diff options
Diffstat (limited to 'libm/float/zetaf.c')
-rw-r--r-- | libm/float/zetaf.c | 175 |
1 files changed, 0 insertions, 175 deletions
diff --git a/libm/float/zetaf.c b/libm/float/zetaf.c deleted file mode 100644 index d01f1d2b2..000000000 --- a/libm/float/zetaf.c +++ /dev/null @@ -1,175 +0,0 @@ -/* zetaf.c - * - * Riemann zeta function of two arguments - * - * - * - * SYNOPSIS: - * - * float x, q, y, zetaf(); - * - * y = zetaf( 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: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE 0,25 10000 6.9e-7 1.0e-7 - * - * Large arguments may produce underflow in powf(), in which - * case the results are inaccurate. - * - * REFERENCE: - * - * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, - * Series, and Products, p. 1073; Academic Press, 1980. - * - */ - -/* -Cephes Math Library Release 2.2: July, 1992 -Copyright 1984, 1987, 1992 by Stephen L. Moshier -Direct inquiries to 30 Frost Street, Cambridge, MA 02140 -*/ - -#include <math.h> -extern float MAXNUMF, MACHEPF; - -/* Expansion coefficients - * for Euler-Maclaurin summation formula - * (2k)! / B2k - * where B2k are Bernoulli numbers - */ -static float 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 */ - - -#define fabsf(x) ( (x) < 0 ? -(x) : (x) ) - - -float powf( float, float ); -float zetaf(float xx, float qq) -{ -int i; -float x, q, a, b, k, s, w, t; - -x = xx; -q = qq; -if( x == 1.0 ) - return( MAXNUMF ); - -if( x < 1.0 ) - { - mtherr( "zetaf", DOMAIN ); - return(0.0); - } - - -/* Euler-Maclaurin summation formula */ -/* -if( x < 25.0 ) -{ -*/ -w = 9.0; -s = powf( q, -x ); -a = q; -for( i=0; i<9; i++ ) - { - a += 1.0; - b = powf( a, -x ); - s += b; - if( b/s < MACHEPF ) - 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 = fabsf(t/s); - if( t < MACHEPF ) - goto done; - k += 1.0; - a *= x + k; - b /= w; - k += 1.0; - } -done: -return(s); -/* -} -*/ - - -/* Basic sum of inverse powers */ -/* -pseres: - -s = powf( q, -x ); -a = q; -do - { - a += 2.0; - b = powf( a, -x ); - s += b; - } -while( b/s > MACHEPF ); - -b = powf( 2.0, -x ); -s = (s + b)/(1.0-b); -return(s); -*/ -} |