fix: prevent non-terminating searches in stats/base/dists/negative-binomial/quantile

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/README.md

@@ -124,6 +124,13 @@ y = quantile( 0.3, 20.0, 1.5 );
 // returns NaN
 ```
 
+If provided a success probability `p` equal to `0`, the function returns `NaN`, as the number of failures before observing `r` successes is unbounded.
+
+```javascript
+var y = quantile( 0.5, 20.0, 0.0 );
+// returns NaN
+```
+
 #### quantile.factory( r, p )
 
 Returns a function for evaluating the [quantile function][quantile-function] for a [negative binomial][negative-binomial-distribution] distribution with number of successes until experiment is stopped `r` and success probability `p`.

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/docs/repl.txt

@@ -14,6 +14,9 @@
     If provided a success probability `p` outside of `[0,1]`, the function
     returns `NaN`.
 
+    If provided a success probability `p` equal to `0`, the function returns
+    `NaN`.
+
     Parameters
     ----------
     k: number
@@ -46,10 +49,10 @@
     > y = {{alias}}( -0.1, 20.0, 0.5 )
     NaN
 
-    > y = {{alias}}( 21.0, 15.5, 0.5 )
-    12
-    > y = {{alias}}( 5.0, 7.4, 0.4 )
-    10
+    > y = {{alias}}( 0.8, 15.5, 0.5 )
+    20
+    > y = {{alias}}( 0.4, 7.4, 0.4 )
+    9
 
     > y = {{alias}}( 0.5, 0.0, 0.5 )
     NaN

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/lib/factory.js

@@ -57,7 +57,7 @@ function factory( r, p ) {
 		isnan( r ) ||
 		isnan( p ) ||
 		r <= 0.0 ||
-		p < 0.0 ||
+		p <= 0.0 ||
 		p > 1.0
 	) {
 		return constantFunction( NaN );
@@ -94,7 +94,10 @@ function factory( r, p ) {
 		if ( k === 1.0 ) {
 			return PINF;
 		}
-
+		// If the mean overflows double-precision floating-point format, the sought quantile saturates at positive infinity:
+		if ( mu === PINF ) {
+			return PINF;
+		}
 		// Cornish-Fisher expansion:
 		if ( k < 0.5 ) {
 			x = -erfcinv( 2.0 * k ) * SQRT2;

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/lib/main.js

@@ -73,6 +73,10 @@ var search = require( './search.js' );
 * // returns NaN
 *
 * @example
+* var y = quantile( 0.5, 20.0, 0.0 );
+* // returns NaN
+*
+* @example
 * var y = quantile( 0.3, 20.0, -1.0 );
 * // returns NaN
 *
@@ -107,7 +111,7 @@ function quantile( k, r, p ) {
 		isnan( p ) ||
 		isnan( k ) ||
 		r <= 0.0 ||
-		p < 0.0 ||
+		p <= 0.0 ||
 		p > 1.0 ||
 		k < 0.0 ||
 		k > 1.0
@@ -122,6 +126,11 @@ function quantile( k, r, p ) {
 	}
 	q = 1.0 - p;
 	mu = ( r * q ) / p;
+
+	// If the mean overflows double-precision floating-point format, the sought quantile saturates at positive infinity:
+	if ( mu === PINF ) {
+		return PINF;
+	}
 	sigma = sqrt( r * q ) / p;
 	sigmaInv = 1.0 / sigma;
 

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/lib/native.js

@@ -67,6 +67,10 @@ var addon = require( './../src/addon.node' );
 * // returns NaN
 *
 * @example
+* var y = quantile( 0.5, 20.0, 0.0 );
+* // returns NaN
+*
+* @example
 * var y = quantile( 0.3, 20.0, -1.0 );
 * // returns NaN
 *

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/lib/search.js

@@ -40,11 +40,18 @@ var methods;
 * @returns {NonNegativeInteger} `k` quantile of the specified distribution
 */
 function searchLeft( x, k, r, p ) {
+	var xp;
 	while ( true ) {
 		if ( x === 0 || cdf( x - 1.0, r, p ) < k ) {
 			return x;
 		}
-		x -= 1;
+		xp = x - 1.0;
+
+		// Guard against non-advancing updates (e.g., when `x` is `NaN` or is so large that decrementing does not change its value) in order to avoid an infinite loop:
+		if ( !( xp < x ) ) {
+			return x;
+		}
+		x = xp;
 	}
 }
 
@@ -58,8 +65,15 @@ function searchLeft( x, k, r, p ) {
 * @returns {NonNegativeInteger} `k` quantile of the specified distribution
 */
 function searchRight( x, k, r, p ) {
+	var xn;
 	while ( true ) {
-		x += 1;
+		xn = x + 1.0;
+
+		// Guard against non-advancing updates (e.g., when `x` is `NaN` or is so large that incrementing does not change its value) in order to avoid an infinite loop:
+		if ( !( xn > x ) ) {
+			return x;
+		}
+		x = xn;
 		if ( cdf( x, r, p ) >= k ) {
 			return x;
 		}

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/src/main.c

@@ -35,8 +35,15 @@
 * @return     `k` quantile of the specified distribution
 */
 static double search_left( double x, const double k, const double r, const double p ) {
+	double xp;
 	while ( x > 0.0 && stdlib_base_dists_negative_binomial_cdf( x - 1.0, r, p ) >= k ) {
-		x -= 1.0;
+		xp = x - 1.0;
+
+		// Guard against non-advancing updates (e.g., when `x` is `NaN` or is so large that decrementing does not change its value) in order to avoid an infinite loop:
+		if ( !( xp < x ) ) {
+			return x;
+		}
+		x = xp;
 	}
 	return x;
 }
@@ -51,9 +58,16 @@ static double search_left( double x, const double k, const double r, const doubl
 * @return     `k` quantile of the specified distribution
 */
 static double search_right( double x, const double k, const double r, const double p ) {
-	x += 1.0;
-	while ( stdlib_base_dists_negative_binomial_cdf( x, r, p ) < k ) {
-		x += 1.0;
+	double xn;
+
+	// Guard against non-advancing updates (e.g., when `x` is `NaN` or is so large that incrementing does not change its value) in order to avoid an infinite loop:
+	xn = x + 1.0;
+	while ( xn > x ) {
+		x = xn;
+		if ( stdlib_base_dists_negative_binomial_cdf( x, r, p ) >= k ) {
+			return x;
+		}
+		xn = x + 1.0;
 	}
 	return x;
 }
@@ -85,7 +99,7 @@ double stdlib_base_dists_negative_binomial_quantile( const double k, const doubl
 		stdlib_base_is_nan( p ) ||
 		stdlib_base_is_nan( k ) ||
 		r <= 0.0 ||
-		p < 0.0 ||
+		p <= 0.0 ||
 		p > 1.0 ||
 		k < 0.0 ||
 		k > 1.0
@@ -100,6 +114,11 @@ double stdlib_base_dists_negative_binomial_quantile( const double k, const doubl
 	}
 	q = 1.0 - p;
 	mu = ( r * q ) / p;
+
+	// If the mean overflows double-precision floating-point format, the sought quantile saturates at positive infinity:
+	if ( mu == STDLIB_CONSTANT_FLOAT64_PINF ) {
+		return STDLIB_CONSTANT_FLOAT64_PINF;
+	}
 	sigma = stdlib_base_sqrt( r * q ) / p;
 	sigmaInv = 1.0 / sigma;
 

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/test/test.factory.js

@@ -141,6 +141,51 @@ tape( 'if provided a `r` which is not a positive integer, the created function a
 	t.end();
 });
 
+tape( 'if provided a success probability `p` equal to `0`, the created function always returns `NaN`', function test( t ) {
+	var quantile;
+	var y;
+
+	quantile = factory( 20.0, 0.0 );
+
+	y = quantile( 0.5 );
+	t.strictEqual( isnan( y ), true, 'returns expected value' );
+
+	y = quantile( 0.3 );
+	t.strictEqual( isnan( y ), true, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'if the distribution mean overflows double-precision floating-point format, the created function returns `+infinity`', function test( t ) {
+	var quantile;
+	var y;
+
+	quantile = factory( 1.0e308, 0.1 );
+
+	y = quantile( 0.3 );
+	t.strictEqual( y, PINF, 'returns expected value' );
+
+	y = quantile( 0.9 );
+	t.strictEqual( y, PINF, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'the created function terminates for very large `r` for which unit increments can no longer advance the search', function test( t ) {
+	var quantile;
+	var y;
+
+	quantile = factory( 1.0e300, 0.5 );
+
+	y = quantile( 0.9 );
+	t.strictEqual( y, 1.0e300, 'returns expected value' );
+
+	y = quantile( 0.1 );
+	t.strictEqual( y, 1.0e300, 'returns expected value' );
+
+	t.end();
+});
+
 tape( 'if provided a valid `r` and `p`, the function returns a function which accurately computes the 0% and 100% quantiles', function test( t ) {
 	var quantile;
 	var y;

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/test/test.main.js

@@ -94,6 +94,42 @@ tape( 'if provided a success probability `p` outside `[0,1]`, the function retur
 	t.end();
 });
 
+tape( 'if provided a success probability `p` equal to `0`, the function returns `NaN`', function test( t ) {
+	var y;
+
+	y = quantile( 0.5, 20, 0.0 );
+	t.strictEqual( isnan( y ), true, 'returns expected value' );
+
+	y = quantile( 0.3, 10, 0.0 );
+	t.strictEqual( isnan( y ), true, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'if the distribution mean overflows double-precision floating-point format, the function returns `+infinity`', function test( t ) {
+	var y;
+
+	y = quantile( 0.3, 1.0e308, 0.1 );
+	t.strictEqual( y, PINF, 'returns expected value' );
+
+	y = quantile( 0.9, 1.0e308, 0.1 );
+	t.strictEqual( y, PINF, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'the function terminates for very large `r` for which unit increments can no longer advance the search', function test( t ) {
+	var y;
+
+	y = quantile( 0.9, 1.0e300, 0.5 );
+	t.strictEqual( y, 1.0e300, 'returns expected value' );
+
+	y = quantile( 0.1, 1.0e300, 0.5 );
+	t.strictEqual( y, 1.0e300, 'returns expected value' );
+
+	t.end();
+});
+
 tape( 'if provided a valid `r` and `p`, the function accurately computes the 0% and 100% quantiles', function test( t ) {
 	var y;
 

lib/node_modules/@stdlib/stats/base/dists/negative-binomial/quantile/test/test.native.js

@@ -103,6 +103,42 @@ tape( 'if provided a success probability `p` outside `[0,1]`, the function retur
 	t.end();
 });
 
+tape( 'if provided a success probability `p` equal to `0`, the function returns `NaN`', opts, function test( t ) {
+	var y;
+
+	y = quantile( 0.5, 20, 0.0 );
+	t.strictEqual( isnan( y ), true, 'returns expected value' );
+
+	y = quantile( 0.3, 10, 0.0 );
+	t.strictEqual( isnan( y ), true, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'if the distribution mean overflows double-precision floating-point format, the function returns `+infinity`', opts, function test( t ) {
+	var y;
+
+	y = quantile( 0.3, 1.0e308, 0.1 );
+	t.strictEqual( y, PINF, 'returns expected value' );
+
+	y = quantile( 0.9, 1.0e308, 0.1 );
+	t.strictEqual( y, PINF, 'returns expected value' );
+
+	t.end();
+});
+
+tape( 'the function terminates for very large `r` for which unit increments can no longer advance the search', opts, function test( t ) {
+	var y;
+
+	y = quantile( 0.9, 1.0e300, 0.5 );
+	t.strictEqual( y, 1.0e300, 'returns expected value' );
+
+	y = quantile( 0.1, 1.0e300, 0.5 );
+	t.strictEqual( y, 1.0e300, 'returns expected value' );
+
+	t.end();
+});
+
 tape( 'if provided a valid `r` and `p`, the function accurately computes the 0% and 100% quantiles', opts, function test( t ) {
 	var y;