gonum.org/v1/gonum@v0.14.0/mathext/internal/cephes/zeta.go (about)

     1  // Derived from SciPy's special/cephes/zeta.c
     2  // https://github.com/scipy/scipy/blob/master/scipy/special/cephes/zeta.c
     3  // Made freely available by Stephen L. Moshier without support or guarantee.
     4  
     5  // Use of this source code is governed by a BSD-style
     6  // license that can be found in the LICENSE file.
     7  // Copyright ©1984, ©1987 by Stephen L. Moshier
     8  // Portions Copyright ©2016 The Gonum Authors. All rights reserved.
     9  
    10  package cephes
    11  
    12  import "math"
    13  
    14  // zetaCoegs are the expansion coefficients for Euler-Maclaurin summation
    15  // formula:
    16  //
    17  //	\frac{(2k)!}{B_{2k}}
    18  //
    19  // where
    20  //
    21  //	B_{2k}
    22  //
    23  // are Bernoulli numbers.
    24  var zetaCoefs = [...]float64{
    25  	12.0,
    26  	-720.0,
    27  	30240.0,
    28  	-1209600.0,
    29  	47900160.0,
    30  	-1.307674368e12 / 691,
    31  	7.47242496e10,
    32  	-1.067062284288e16 / 3617,
    33  	5.109094217170944e18 / 43867,
    34  	-8.028576626982912e20 / 174611,
    35  	1.5511210043330985984e23 / 854513,
    36  	-1.6938241367317436694528e27 / 236364091,
    37  }
    38  
    39  // Zeta computes the Riemann zeta function of two arguments.
    40  //
    41  //	Zeta(x,q) = \sum_{k=0}^{\infty} (k+q)^{-x}
    42  //
    43  // Note that Zeta returns +Inf if x is 1 and will panic if x is less than 1,
    44  // q is either zero or a negative integer, or q is negative and x is not an
    45  // integer.
    46  //
    47  // Note that:
    48  //
    49  //	zeta(x,1) = zetac(x) + 1
    50  func Zeta(x, q float64) float64 {
    51  	// REFERENCE: Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, Series,
    52  	// and Products, p. 1073; Academic Press, 1980.
    53  	if x == 1 {
    54  		return math.Inf(1)
    55  	}
    56  
    57  	if x < 1 {
    58  		panic(paramOutOfBounds)
    59  	}
    60  
    61  	if q <= 0 {
    62  		if q == math.Floor(q) {
    63  			panic(errParamFunctionSingularity)
    64  		}
    65  		if x != math.Floor(x) {
    66  			panic(paramOutOfBounds) // Because q^-x not defined
    67  		}
    68  	}
    69  
    70  	// Asymptotic expansion: http://dlmf.nist.gov/25.11#E43
    71  	if q > 1e8 {
    72  		return (1/(x-1) + 1/(2*q)) * math.Pow(q, 1-x)
    73  	}
    74  
    75  	// The Euler-Maclaurin summation formula is used to obtain the expansion:
    76  	//  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}}
    77  	// where
    78  	//  B_{2j}
    79  	// are Bernoulli numbers.
    80  	// Permit negative q but continue sum until n+q > 9. This case should be
    81  	// handled by a reflection formula. If q<0 and x is an integer, there is a
    82  	// relation to the polyGamma function.
    83  	s := math.Pow(q, -x)
    84  	a := q
    85  	i := 0
    86  	b := 0.0
    87  	for i < 9 || a <= 9 {
    88  		i++
    89  		a += 1.0
    90  		b = math.Pow(a, -x)
    91  		s += b
    92  		if math.Abs(b/s) < machEp {
    93  			return s
    94  		}
    95  	}
    96  
    97  	w := a
    98  	s += b * w / (x - 1)
    99  	s -= 0.5 * b
   100  	a = 1.0
   101  	k := 0.0
   102  	for _, coef := range zetaCoefs {
   103  		a *= x + k
   104  		b /= w
   105  		t := a * b / coef
   106  		s = s + t
   107  		t = math.Abs(t / s)
   108  		if t < machEp {
   109  			return s
   110  		}
   111  		k += 1.0
   112  		a *= x + k
   113  		b /= w
   114  		k += 1.0
   115  	}
   116  	return s
   117  }