summaryrefslogtreecommitdiff
path: root/libm/double/zeta.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/zeta.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/zeta.c')
-rw-r--r--libm/double/zeta.c189
1 files changed, 189 insertions, 0 deletions
diff --git a/libm/double/zeta.c b/libm/double/zeta.c
new file mode 100644
index 000000000..a49c619d5
--- /dev/null
+++ b/libm/double/zeta.c
@@ -0,0 +1,189 @@
+/* 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);
+*/
+}