From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001
From: Eric Andersen <andersen@codepoet.org>
Date: Thu, 22 Nov 2001 14:04:29 +0000
Subject: 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

---
 libm/double/zeta.c | 189 -----------------------------------------------------
 1 file changed, 189 deletions(-)
 delete mode 100644 libm/double/zeta.c

(limited to 'libm/double/zeta.c')

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);
-*/
-}
-- 
cgit v1.2.3