/**
* @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_53_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf_gamma/digamma.html}. 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 floor = require( './../../../../base/special/floor' );
var tan = require( './../../../../base/special/tan' );
var PI = require( '@stdlib/constants/float64/pi' );
var asymptoticApprox = require( './asymptotic_expansion.js' );
var rationalApprox = require( './rational_approximation.js' );


// VARIABLES //

var MIN_SAFE_ASYMPTOTIC = 10.0; // BIG!


// MAIN //

/**
* Evaluates the digamma function.
*
* ## Method
*
* 1.  For \\(x < 0\\), we use the reflection formula
*
*     ```tex
*     \psi(1-x) = \psi(x) + \frac{\pi}{\tan(\pi x)}
*     ```
*
*     to make \\(x\\) positive.
*
* 2.  For \\(x \in \[0,1]\\), we use the recurrence relation
*
*     ```tex
*     \psi(x) = \psi(x+1) - \frac{1}{x}
*     ```
*
*     to shift the evaluation range to \\(\[1,2]\\).
*
* 3.  For \\(x \in \[1,2]\\), we use a rational approximation of the form
*
*     ```tex
*     \psi(x) = (x - \mathrm{root})(Y + \operatorname{R}(x-1))
*     ```
*
*     where \\(\mathrm{root}\\) is the location of the positive root of \\(\psi\\), \\(Y\\) is a constant, and \\(R\\) is optimized for low absolute error compared to \\(Y\\).
*
*     <!-- <note>-->
*
*     Note that, since \\(\mathrm{root}\\) is irrational, we need twice as many digits in \\(\mathrm{root}\\) as in \\(x\\) in order to avoid cancellation error during subtraction, assuming \\(x\\) has an exact value. This means that, even if \\(x\\) is rounded to the next representable value, the result of \\(\psi(x)\\) will not be zero.
*
*     <!-- </note> -->
*
*     <!-- <note> -->
*
*     This approach gives 17-digit precision.
*
*     <!-- </note> -->
*
* 4.  For \\(x \in \[2,\mathrm{BIG}]\\), we use the recurrence relation
*
*     ```tex
*     \psi(x+1) = \psi(x) + \frac{1}{x}
*     ```
*
*     to shift the evaluation range to \\(\[1,2]\\).
*
* 5.  For \\(x > \mathrm{BIG}\\), we use the asymptotic expression
*
*     ```tex
*     \psi(x) = \ln(x) + \frac{1}{2x} - \biggl( \frac{B_{21}}{2x^2} + \frac{B_{22}}{4x^4} + \frac{B_{23}}{6x^6} + \ldots \biggr)
*     ```
*
*     This expansion, however, is divergent after a few terms. The number of terms depends on \\(x\\). Accordingly, we must choose a value of \\(\mathrm{BIG}\\) which allows us to truncate the series at a term that is too small to have an effect on the result. Setting \\(\mathrm{BIG} = 10\\), allows us to truncate the series early and evaluate as \\(1/x^2\\).
*
*     <!-- <note> -->
*
*     This approach gives 17-digit precision for \\(x \geq 10\\).
*
*     <!-- </note> -->
*
* ## Notes
*
* -   Maximum deviation found: \\(1.466\\mbox{e-}18\\)
* -   Max error found: \\(2.452\mbox{e-}17\\) (double precision)
*
*
* @param {number} x - input value
* @returns {number} function value
*
* @example
* var v = digamma( -2.5 );
* // returns ~1.103
*
* @example
* var v = digamma( 1.0 );
* // returns ~-0.577
*
* @example
* var v = digamma( 10.0 );
* // returns ~2.252
*
* @example
* var v = digamma( NaN );
* // returns NaN
*
* @example
* var v = digamma( -1.0 );
* // returns NaN
*/
function digamma( x ) {
	var rem;
	var tmp;
	if ( isnan( x ) || x === 0.0 ) {
		return NaN;
	}
	// If `x` is negative, use reflection...
	if ( x <= -1.0 ) {
		// Reflect:
		x = 1.0 - x;

		// Argument reduction for tan:
		rem = x - floor(x);

		// Shift to negative if > 0.5:
		if ( rem > 0.5 ) {
			rem -= 1.0;
		}
		// Check for evaluation at a negative pole:
		if ( rem === 0.0 ) {
			return NaN;
		}
		tmp = PI / tan( PI * rem );
	} else {
		tmp = 0.0;
	}
	// If we're above the lower-limit for the asymptotic expansion, then use it...
	if ( x >= MIN_SAFE_ASYMPTOTIC ) {
		tmp += asymptoticApprox( x );
		return tmp;
	}
	// If x > 2, reduce to the interval [1,2]...
	while ( x > 2.0 ) {
		x -= 1.0;
		tmp += 1.0/x;
	}
	// If x < 1, use recurrence to shift to > 1..
	while ( x < 1.0 ) {
		tmp -= 1.0/x;
		x += 1.0;
	}
	tmp += rationalApprox( x );
	return tmp;
}


// EXPORTS //

module.exports = digamma;