time-to-botec

zeta.js (8087B)

---

 /**
 * @license Apache-2.0
 *
 * Copyright (c) 2018 The Stdlib Authors.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 *
 * ## Notice
 *
 * The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_60_0/boost/math/special_functions/zeta.hpp}. The implementation follows the original, but has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
 *
 * Use, modification and distribution are subject to the
 * Boost Software License, Version 1.0. (See accompanying file
 * LICENSE or copy at http://www.boost.org/LICENSE_1_0.txt)
 * ```
 */
 
 'use strict';
 
 // MODULES //
 
 var isnan = require( './../../../../base/assert/is-nan' );
 var isInteger = require( './../../../../base/assert/is-integer' );
 var abs = require( './../../../../base/special/abs' );
 var exp = require( './../../../../base/special/exp' );
 var floor = require( './../../../../base/special/floor' );
 var gamma = require( './../../../../base/special/gamma' );
 var gammaln = require( './../../../../base/special/gammaln' );
 var sinpi = require( './../../../../base/special/sinpi' );
 var pow = require( './../../../../base/special/pow' );
 var ln = require( './../../../../base/special/ln' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
 var NINF = require( '@stdlib/constants/float64/ninf' );
 var TWO_PI = require( '@stdlib/constants/float64/two-pi' );
 var SQRT_EPSILON = require( '@stdlib/constants/float64/sqrt-eps' );
 var LN_SQRT_TWO_PI = require( '@stdlib/constants/float64/ln-sqrt-two-pi' );
 var ODD_POSITIVE_INTEGERS = require( './odd_positive_integers.json' );
 var EVEN_NONNEGATIVE_INTEGERS = require( './even_nonnegative_integers.json' );
 var BERNOULLI = require( './bernoulli.json' );
 var rateval1 = require( './rational_p1q1.js' );
 var rateval2 = require( './rational_p2q2.js' );
 var rateval3 = require( './rational_p3q3.js' );
 var rateval4 = require( './rational_p4q4.js' );
 var rateval5 = require( './rational_p5q5.js' );
 var rateval6 = require( './rational_p6q6.js' );
 
 
 // VARIABLES //
 
 var MAX_BERNOULLI_2N = 129;
 var MAX_FACTORIAL = 170; // TODO: consider making external constant
 var MAX_LN = 709; // TODO: consider making external constant
 var Y1 = 1.2433929443359375;
 var Y3 = 0.6986598968505859375;
 
 
 // MAIN //
 
 /**
 * Evaluates the Riemann zeta function.
 *
 * ## Method
 *
 * 1.  First, we use the reflection formula
 *
 *     ```tex
 *     \zeta(1-s) = 2 \sin\biggl(\frac{\pi(1-s)}{2}\biggr)(2\pi^{-s})\Gamma(s)\zeta(s)
 *     ```
 *
 *     to make \\(s\\) positive.
 *
 * 2.  For \\(s \in (0,1)\\), we use the approximation
 *
 *     ```tex
 *     \zeta(s) = \frac{C + \operatorname{R}(1-s) - s}{1-s}
 *     ```
 *
 *     with rational approximation \\(\operatorname{R}(1-z)\\) and constant \\(C\\).
 *
 * 3.  For \\(s \in (1,4)\\), we use the approximation
 *
 *     ```tex
 *     \zeta(s) = C + \operatorname{R}(s-n) + \frac{1}{s-1}
 *     ```
 *
 *     with rational approximation \\(\operatorname{R}(z-n)\\), constant \\(C\\), and integer \\(n\\).
 *
 * 4.  For \\(s > 4\\), we use the approximation
 *
 *     ```tex
 *     \zeta(s) = 1 + e^{\operatorname{R}(z-n)}
 *     ```
 *
 *     with rational approximation \\(\operatorname{R}(z-n)\\) and integer \\(n\\).
 *
 * 5.  For negative odd integers, we use the closed form
 *
 *     ```tex
 *     \zeta(-n) = \frac{(-1)^n}{n+1} B_{n+1}
 *     ```
 *
 *     where \\(B_{n+1}\\) is a Bernoulli number.
 *
 * 6.  For negative even integers, we use the closed form
 *
 *     ```tex
 *     \zeta(-2n) = 0
 *     ```
 *
 * 7.  For nonnegative even integers, we could use the closed form
 *
 *     ```tex
 *     \zeta(2n) = \frac{(-1)^{n-1}2^{2n-1}\pi^{2n}}{(2n)!} B_{2n}
 *     ```
 *
 *     where \\(B_{2n}\\) is a Bernoulli number. However, to speed computation, we use precomputed values (Wolfram Alpha).
 *
 * 8.  For positive negative integers, we use precomputed values (Wolfram Alpha), as these values are useful for certain infinite series calculations.
 *
 *
 * ## Notes
 *
 * -   \\(\[\approx 1.5\mbox{e-}8, 1)\\)
 *
 *     -   max deviation: \\(2.020\mbox{e-}18\\)
 *     -   expected error: \\(-2.020\mbox{e-}18\\)
 *     -   max error found (double): \\(3.994987\mbox{e-}17\\)
 *
 * -   \\(\[1,2\]\\)
 *
 *     -   max deviation: \\(9.007\mbox{e-}20\\)
 *     -   expected error: \\(9.007\mbox{e-}20\\)
 *
 * -   \\((2,4\]\\)
 *
 *     -   max deviation: \\(5.946\mbox{e-}22\\)
 *     -   expected error: \\(-5.946\mbox{e-}22\\)
 *
 * -   \\((4,7\]\\)
 *
 *     -   max deviation: \\(2.955\mbox{e-}17\\)
 *     -   expected error: \\(2.955\mbox{e-}17\\)
 *     -   max error found (double): \\(2.009135\mbox{e-}16\\)
 *
 * -   \\((7,15)\\)
 *
 *     -   max deviation: \\(7.117\mbox{e-}16\\)
 *     -   expected error: \\(7.117\mbox{e-}16\\)
 *     -   max error found (double): \\(9.387771\mbox{e-}16\\)
 *
 * -   \\(\[15,36)\\)
 *
 *     -   max error (in interpolated form): \\(1.668\mbox{e-}17\\)
 *     -   max error found (long double): \\(1.669714\mbox{e-}17\\)
 *
 *
 * @param {number} s - input value
 * @returns {number} function value
 *
 * @example
 * var v = zeta( 1.1 );
 * // returns ~10.584
 *
 * @example
 * var v = zeta( -4.0 );
 * // returns 0.0
 *
 * @example
 * var v = zeta( 70.0 );
 * // returns 1.0
 *
 * @example
 * var v = zeta( 0.5 );
 * // returns ~-1.46
 *
 * @example
 * var v = zeta( 1.0 ); // pole
 * // returns NaN
 *
 * @example
 * var v = zeta( NaN );
 * // returns NaN
 */
 function zeta( s ) {
 	var tmp;
 	var sc;
 	var as;
 	var is;
 	var r;
 	var n;
 
 	// Check for `NaN`:
 	if ( isnan( s ) ) {
 		return NaN;
 	}
 	// Check for a pole:
 	if ( s === 1.0 ) {
 		return NaN;
 	}
 	// Check for large value:
 	if ( s >= 56.0 ) {
 		return 1.0;
 	}
 	// Check for a closed form (integers):
 	if ( isInteger( s ) ) {
 		// Cast `s` to a 32-bit signed integer:
 		is = s|0; // asm type annotation
 
 		// Check that `s` does not exceed MAX_INT32:
 		if ( is === s ) {
 			if ( is < 0 ) {
 				as = (-is)|0; // asm type annotation
 
 				// Check if even negative integer:
 				if ( (as&1) === 0 ) {
 					return 0.0;
 				}
 				n = ( (as+1) / 2 )|0; // asm type annotation
 
 				// Check if less than max Bernoulli number:
 				if ( n <= MAX_BERNOULLI_2N ) {
 					return -BERNOULLI[ n ] / (as+1.0);
 				}
 				// fall through...
 			}
 			// Check if even nonnegative integer:
 			else if ( (is&1) === 0 ) {
 				return EVEN_NONNEGATIVE_INTEGERS[ is/2 ];
 			}
 			// Must be a odd positive integer:
 			else {
 				return ODD_POSITIVE_INTEGERS[ (is-3)/2 ];
 			}
 		}
 		// fall through...
 	}
 	if ( abs(s) < SQRT_EPSILON ) {
 		return -0.5 - (LN_SQRT_TWO_PI * s);
 	}
 	sc = 1.0 - s;
 	if ( s < 0.0 ) {
 		// Check if even negative integer:
 		if ( floor(s/2.0) === s/2.0 ) {
 			return 0.0;
 		}
 		// Swap `s` and `sc`:
 		tmp = s;
 		s = sc;
 		sc = tmp;
 
 		// Determine if computation will overflow:
 		if ( s > MAX_FACTORIAL ) {
 			tmp = sinpi( 0.5*sc ) * 2.0 * zeta( s );
 			r = gammaln( s );
 			r -= s * ln( TWO_PI );
 			if ( r > MAX_LN ) {
 				return ( tmp < 0.0 ) ? NINF : PINF;
 			}
 			return tmp * exp( r );
 		}
 		return sinpi( 0.5*sc ) * 2.0 * pow( TWO_PI, -s ) * gamma( s ) * zeta( s ); // eslint-disable-line max-len
 	}
 	if ( s < 1.0 ) {
 		tmp = rateval1( sc );
 		tmp -= Y1;
 		tmp += sc;
 		tmp /= sc;
 		return tmp;
 	}
 	if ( s <= 2.0 ) {
 		sc = -sc;
 		tmp = 1.0 / sc;
 		return tmp + rateval2( sc );
 	}
 	if ( s <= 4.0 ) {
 		tmp = Y3 + ( 1.0 / (-sc) );
 		return tmp + rateval3( s-2.0 );
 	}
 	if ( s <= 7.0 ) {
 		tmp = rateval4( s-4.0 );
 		return 1.0 + exp( tmp );
 	}
 	if ( s < 15.0 ) {
 		tmp = rateval5( s-7.0 );
 		return 1.0 + exp( tmp );
 	}
 	if ( s < 36.0 ) {
 		tmp = rateval6( s-15.0 );
 		return 1.0 + exp( tmp );
 	}
 	// s < 56
 	return 1.0 + pow( 2.0, -s );
 }
 
 
 // EXPORTS //
 
 module.exports = zeta;