/**
* @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_64_0/boost/math/special_functions/factorials.hpp}. The implementation has been modified for JavaScript.
*
* ```text
* (C) Copyright John Maddock 2006, 2010.
*
* 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 isNonNegativeInteger = require( './../../../../base/assert/is-nonnegative-integer' );
var isInteger = require( './../../../../base/assert/is-integer' );
var isnan = require( './../../../../base/assert/is-nan' );
var gammaDeltaRatio = require( './../../../../base/special/gamma-delta-ratio' );
var floor = require( './../../../../base/special/floor' );
var abs = require( './../../../../base/special/abs' );
var FLOAT64_MAX = require( '@stdlib/constants/float64/max' );
var PINF = require( '@stdlib/constants/float64/pinf' );


// VARIABLES //

var MAX_FACTORIAL = 170; // TODO: consider packaging as constant


// FUNCTIONS //

/**
* Computes the rising factorial of `x` and `n`.
*
* ## Notes
*
* -   The rising factorial is defined as
*
*     ```tex
*     \operatorname{risingFactorial}(x, n) = x (x-1) (x-2) (x-3) \ldots (x-n+1)
*     ```
*
*     or equivalently
*
*     ```tex
*     \operatorname{risingFactorial}(x, n) = \frac{ \Gamma(x + n) }{ \Gamma(x) };
*     ```
*
* @param {number} x - first function parameter
* @param {integer} n - second function parameter
* @returns {number} function value
*
* @example
* var v = risingFactorial( 0.9, 5 );
* // returns ~94.766
*
* @example
* var v = risingFactorial( -9.0, 3 );
* // returns -504.0
*
* @example
* var v = risingFactorial( 0.0, 2 );
* // returns 0.0
*
* @example
* var v = risingFactorial( 3.0, -2 );
* // returns 0.5
*/
function risingFactorial( x, n ) {
	var result;
	var inv;

	if ( isnan( x ) || !isInteger( n ) ) {
		return NaN;
	}
	if ( x < 0.0 ) {
		// For `x < 0`, we really have a falling factorial, modulo a possible change of sign. Note that the falling factorial isn't defined for negative `n`, so we'll get rid of that case first:
		if ( n < 0.0 ) {
			x += n;
			n = -n;
			inv = true;
		}
		result = ( (n&1) ? -1.0 : 1.0 ) * fallingFactorial( -x, n );
		if ( inv ) {
			result = 1.0 / result;
		}
		return result;
	}
	if ( n === 0 ) {
		return 1.0;
	}
	if ( x === 0.0 ) {
		if ( n < 0 ) {
			return -gammaDeltaRatio( x+1.0, -n );
		}
		return 0.0;
	}
	if ( x < 1.0 && x+n < 0.0 ) {
		result = gammaDeltaRatio( 1.0-x, -n );
		return ( n&1 ) ? -result : result;
	}
	// We don't optimize this for small `n`, because `gammaDeltaRatio` is already optimized for that use case:
	return 1.0 / gammaDeltaRatio( x, n );
}


// MAIN //

/**
* Computes the falling factorial of `x` and `n`.
*
* ## Notes
*
* -   The falling factorial is defined as
*
*     ```tex
*     \operatorname{fallingFactorial}(x, n) = x (x-1) (x-2) (x-3) \ldots (x-n+1)
*     ```
*
* @param {number} x - first function parameter
* @param {NonNegativeInteger} n - second function parameter
* @returns {number} function value
*
* @example
* var v = fallingFactorial( 0.9, 5 );
* // returns ~0.644
*
* @example
* var v = fallingFactorial( -9.0, 3 );
* // returns -990.0
*
* @example
* var v = fallingFactorial( 0.0, 2 );
* // returns 0.0
*
* @example
* var v = fallingFactorial( 3.0, -2 );
* // returns NaN
*/
function fallingFactorial( x, n ) {
	var result;
	var xp1;
	var n2;
	var t1;
	var t2;

	if ( isnan( x ) || !isNonNegativeInteger( n ) ) {
		return NaN;
	}
	if ( x === 0.0 ) {
		return 0.0;
	}
	if ( x < 0.0 ) {
		// For `x < 0`, we really have a rising factorial modulo a possible change of sign:
		return ( ( n&1 ) ? -1.0 : 1.0 ) * risingFactorial( -x, n );
	}
	if ( n === 0 ) {
		return 1.0;
	}
	if ( x < 0.5 ) {
		// Computing `1 + x` will throw away digits, so split up calculation...
		if ( n > MAX_FACTORIAL-2 ) {
			// Given a ratio of two very large numbers, we need to split the calculation up into two blocks:
			t1 = x * fallingFactorial( x-1.0, MAX_FACTORIAL-2 );
			t2 = fallingFactorial( x-MAX_FACTORIAL+1.0, n-MAX_FACTORIAL+1 );
			if ( FLOAT64_MAX/abs(t1) < abs(t2) ) {
				return PINF;
			}
			return t1 * t2;
		}
		return x * fallingFactorial( x-1.0, n-1.0 );
	}
	if ( x <= n-1.0 ) {
		// `x+1-n` will be negative and computing the ratio of two gammas will not work, so split the product up into three parts:
		xp1 = x + 1.0;
		n2 = abs( floor( xp1 ) );
		if ( n2 === xp1 ) {
			return 0.0;
		}
		result = gammaDeltaRatio( xp1, -n2 );
		x -= n2;
		result *= x;
		n2 += 1.0;
		if ( n2 < n ) {
			result *= fallingFactorial( x-1.0, n-n2 );
		}
		return result;
	}
	// Simple case: just the ratio of two (positive argument) gamma functions. Note that we don't optimize this for small `n`, because `gammaDeltaRatio` is already optimized for that use case:
	return gammaDeltaRatio( x+1.0, -n );
}


// EXPORTS //

module.exports = fallingFactorial;