diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/main.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/main.js
index 1916cd504e3f..59073e62d91e 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/main.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/main.js
@@ -48,6 +48,7 @@ var FLOAT64_MAX = require( '@stdlib/constants/float64/max' );
 var SQRT_TWO_PI = require( '@stdlib/constants/float64/sqrt-two-pi' );
 var MAX_LN = require( '@stdlib/constants/float64/max-ln' );
 var PINF = require( '@stdlib/constants/float64/pinf' );
+var FLOAT64_MAX_SAFE_NTH_FACTORIAL = require( '@stdlib/constants/float64/max-safe-nth-factorial' ); // eslint-disable-line id-length
 var finiteGammaQ = require( './finite_gamma_q.js' );
 var finiteHalfGammaQ = require( './finite_half_gamma_q.js' );
 var fullIGammaPrefix = require( './full_igamma_prefix.js' );
@@ -58,11 +59,6 @@ var tgammaSmallUpperPart = require( './tgamma_small_upper_part.js' );
 var upperGammaFraction = require( './upper_gamma_fraction.js' );
 
 
-// VARIABLES //
-
-var MAX_FACTORIAL = 170; // TODO: consider extracting as a constant
-
-
 // MAIN //
 
 /**
@@ -70,7 +66,7 @@ var MAX_FACTORIAL = 170; // TODO: consider extracting as a constant
 *
 * ## Notes
 *
-* -   When `a >= MAX_FACTORIAL` and computing the non-normalized incomplete gamma, result is rather hard to compute unless we use logs. There are really two options a) if `x` is a long way from `a` in value then we can reliably use methods 2 and 4 below in logarithmic form and go straight to the result. Otherwise we let the regularized gamma take the strain (the result is unlikely to underflow in the central region anyway) and combine with `lgamma` in the hopes that we get a finite result.
+* -   When `a >= FLOAT64_MAX_SAFE_NTH_FACTORIAL` and computing the non-normalized incomplete gamma, result is rather hard to compute unless we use logs. There are really two options a) if `x` is a long way from `a` in value then we can reliably use methods 2 and 4 below in logarithmic form and go straight to the result. Otherwise we let the regularized gamma take the strain (the result is unlikely to underflow in the central region anyway) and combine with `lgamma` in the hopes that we get a finite result.
 *
 * @param {NonNegativeNumber} x - function parameter
 * @param {PositiveNumber} a - function parameter
@@ -101,7 +97,7 @@ function gammainc( x, a, regularized, upper ) {
 	normalized = ( regularized === void 0 ) ? true : regularized;
 	invert = upper;
 	result = 0.0;
-	if ( a >= MAX_FACTORIAL && !normalized ) {
+	if ( a >= FLOAT64_MAX_SAFE_NTH_FACTORIAL && !normalized ) {
 		if ( invert && ( a * 4.0 < x ) ) {
 			// This is method 4 below, done in logs:
 			result = ( a * ln(x) ) - x;