Revising Zeta comments and copyright

Author: zeroviscosity

internal/cephes/zeta.go

@@ -1,81 +1,53 @@
-// Copyright ©2016 The gonum Authors. All rights reserved.
+// Derived from SciPy's special/cephes/zeta.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/zeta.c
+
 // Use of this source code is governed by a BSD-style
 // license that can be found in the LICENSE file.
-
-/*
- * Cephes Math Library Release 2.1:  January, 1989
- * Copyright 1984, 1987, 1989 by Stephen L. Moshier
- * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
- */
+// Copyright ©1984, ©1987 by Stephen L. Moshier
+// Portions Copyright ©2016 The gonum Authors. All rights reserved.
 
 package cephes
 
 import "math"
 
-/*
- * Adapted from scipy's cephes zeta.c
- */
-
-/*
- *     Riemann zeta function of two arguments
- *
- *
- * 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.
- *
- *
- * REFERENCE:
- *
- * Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
- * Series, and Products, p. 1073; Academic Press, 1980.
- *
- */
-
-/* Expansion coefficients
- * for Euler-Maclaurin summation formula
- * (2k)! / B2k
- * where B2k are Bernoulli numbers
- */
+// zetaCoegs are the expansion coefficients for Euler-Maclaurin summation
+// formula:
+//  \frac{(2k)!}{B_{2k}}
+// where
+//  B_{2k}
+// are Bernoulli numbers.
 var zetaCoefs = []float64{
 	12.0,
 	-720.0,
 	30240.0,
 	-1209600.0,
 	47900160.0,
-	-1.8924375803183791606e9, /*1.307674368e12/691 */
+	-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 */
+	-1.067062284288e16 / 3617,
+	5.109094217170944e18 / 43867,
+	-8.028576626982912e20 / 174611,
+	1.5511210043330985984e23 / 854513,
+	-1.6938241367317436694528e27 / 236364091,
 }
 
-// Zeta calculates the Riemann zeta function of two arguments
+// Zeta computes the Riemann zeta function of two arguments.
+//  Zeta(x,q) = \sum_{k=0}^{\infty} (k+q)^{-x}
+// Note that Zeta returns +Inf if x is 1 and will panic if x is less than 1,
+// q is either zero or a negative integer, or q is negative and x is not an
+// integer.
+//
+// The Euler-Maclaurin summation formula is used to obtain the expansion:
+//  Zeta(x,q) = \sum_{k=1}^n (k+q)^{-x} + \frac{(n+q)^{1-x}}{x-1} - \frac{1}{2(n+q)^x} + \sum_{j=1}^{\infty} \frac{B_{2j}x(x+1)...(x+2j)}{(2j)! (n+q)^{x+2j+1}}
+// where
+//  B_{2j}
+// are Bernoulli numbers.
+//
+// Note that:
+//  zeta(x,1) = zetac(x) + 1
+//
+// REFERENCE: Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, Series,
+// and Products, p. 1073; Academic Press, 1980.
 func Zeta(x, q float64) float64 {
 	if x == 1 {
 		return math.MaxFloat64
@@ -94,20 +66,16 @@ func Zeta(x, q float64) float64 {
 		}
 	}
 
-	/* Asymptotic expansion
-	 * http://dlmf.nist.gov/25.11#E43
-	 */
+	// Asymptotic expansion: http://dlmf.nist.gov/25.11#E43
 	if q > 1e8 {
 		return (1/(x-1) + 1/(2*q)) * math.Pow(q, 1-x)
 	}
 
 	// Euler-Maclaurin summation formula
 
-	/* 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.
-	 */
+	// 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 := math.Pow(q, -x)
 	a := q
 	i := 0

zeta.go

@@ -8,7 +8,9 @@ import "github.com/gonum/mathext/internal/cephes"
 
 // Zeta computes the Riemann zeta function of two arguments.
 //  Zeta(x,q) = \sum_{k=0}^{\infty} (k+q)^{-x}
-// where x > 1 and q is not a negative integer or zero.
+// Note that Zeta returns +Inf if x is 1 and will panic if x is less than 1,
+// q is either zero or a negative integer, or q is negative and x is not an
+// integer.
 //
 // See http://mathworld.wolfram.com/HurwitzZetaFunction.html
 // or https://en.wikipedia.org/wiki/Multiple_zeta_function#Two_parameters_case