summaryrefslogtreecommitdiff
path: root/libm/float/zetaf.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/float/zetaf.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/float/zetaf.c')
-rw-r--r--libm/float/zetaf.c175
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);
-*/
-}