/**
* @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;