diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/docs/types/index.d.ts b/lib/node_modules/@stdlib/math/base/special/gammainc/docs/types/index.d.ts
index 92a4b55936fe..698d225e7d5d 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/docs/types/index.d.ts
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/docs/types/index.d.ts
@@ -33,22 +33,27 @@
 * @returns function value
 *
 * @example
-* var y = gammainc( 6.0, 2.0 )
+* var y = gammainc( 6.0, 2.0 );
 * // returns ~0.9826
+*
 * @example
-* var y = gammainc( 1.0, 2.0, true, true )
+* var y = gammainc( 1.0, 2.0, true, true );
 * // returns ~0.7358
+*
 * @example
-* var y = gammainc( 7.0, 5.0 )
+* var y = gammainc( 7.0, 5.0 );
 * // returns ~0.8270
+*
 * @example
-* var y = gammainc( 7.0, 5.0, false )
+* var y = gammainc( 7.0, 5.0, false );
 * // returns ~19.8482
+*
 * @example
-* var y = gammainc( NaN, 2.0 )
+* var y = gammainc( NaN, 2.0 );
 * // returns NaN
+*
 * @example
-* var y = gammainc( 6.0, NaN )
+* var y = gammainc( 6.0, NaN );
 * // returns NaN
 */
 declare function gammainc( x: number, a: number, regularized?: boolean, upper?: boolean ): number;
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_gamma_q.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_gamma_q.js
index 6d450f58dd7c..b70347d4e3ad 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_gamma_q.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_gamma_q.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_half_gamma_q.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_half_gamma_q.js
index 1dcbd6e2c407..697594451b52 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_half_gamma_q.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/finite_half_gamma_q.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/full_igamma_prefix.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/full_igamma_prefix.js
index ab0ae1793a9b..c8281a747d1b 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/full_igamma_prefix.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/full_igamma_prefix.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
@@ -59,20 +59,16 @@ function fullIGammaPrefix( a, z ) {
 	if ( z >= 1.0 ) {
 		if ( ( alz < MAX_LN ) && ( -z > MIN_LN ) ) {
 			prefix = pow( z, a ) * exp( -z );
-		}
-		else if ( a >= 1.0 ) {
+		} else if ( a >= 1.0 ) {
 			prefix = pow( z / exp(z/a), a );
-		}
-		else {
+		} else {
 			prefix = exp( alz - z );
 		}
-	}
-	else {
+	} else {
 		/* eslint-disable no-lonely-if */
 		if ( alz > MIN_LN ) {
 			prefix = pow( z, a ) * exp( -z );
-		}
-		else if ( z/a < MAX_LN ) {
+		} else if ( z/a < MAX_LN ) {
 			prefix = pow( z / exp(z/a), a );
 		} else {
 			prefix = exp( alz - z );
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/igamma_final.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/igamma_final.js
new file mode 100644
index 000000000000..9407bb1dbb0c
--- /dev/null
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/igamma_final.js
@@ -0,0 +1,258 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2025 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 https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+*
+* ```text
+* (C) Copyright John Maddock 2006.
+* (C) Copyright Paul A. Bristow 2007.
+*
+* 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 floor = require( '@stdlib/math/base/special/floor' );
+var gamma = require( '@stdlib/math/base/special/gamma' );
+var abs = require( '@stdlib/math/base/special/abs' );
+var pow = require( '@stdlib/math/base/special/pow' );
+var ln = require( '@stdlib/math/base/special/ln' );
+var SQRT_EPSILON = require( '@stdlib/constants/float64/sqrt-eps' );
+var FLOAT64_MAX = require( '@stdlib/constants/float64/max' );
+var MAX_LN = require( '@stdlib/constants/float64/max-ln' );
+var tgammaILargeX = require( './tgamma_i_large_x.js' );
+var finiteGammaQ = require( './finite_gamma_q.js' );
+var finiteHalfGammaQ = require( './finite_half_gamma_q.js' );
+var fullIGammaPrefix = require( './full_igamma_prefix.js' );
+var igammaTemmeLarge = require( './igamma_temme_large.js' );
+var lowerGammaSeries = require( './lower_gamma_series.js' );
+var regularisedGammaPrefix = require( './regularised_gamma_prefix.js' );
+var tgammaSmallUpperPart = require( './tgamma_small_upper_part.js' );
+var upperGammaFraction = require( './upper_gamma_fraction.js' );
+
+
+// MAIN //
+
+/**
+* Main entry point for evaluating all the four incomplete gamma functions (lower, upper, regularized lower, and regularized upper).
+*
+* @private
+* @param {number} x - function parameter
+* @param {number} a - function parameter
+* @param {boolean} regularized - boolean indicating if the function should evaluate the regularized or non-regularized incomplete gamma functions
+* @param {boolean} upper - boolean indicating if the function should return the upper tail of the incomplete gamma function
+* @returns {number} function value
+*/
+function igammaFinal( x, a, regularized, upper ) {
+	var optimisedInvert;
+	var evalMethod;
+	var isHalfInt;
+	var initValue;
+	var useTemme;
+	var isSmallA;
+	var result;
+	var invert;
+	var isInt;
+	var sigma;
+	var res;
+	var gam;
+	var fa;
+	var g;
+
+	result = 0.0;
+	invert = upper;
+	isSmallA = ( a < 30 ) && ( a <= x + 1.0 ) && ( x < MAX_LN );
+	if ( isSmallA ) {
+		fa = floor( a );
+		isInt = ( fa === a );
+		isHalfInt = ( isInt ) ? false : ( abs( fa - a ) === 0.5 );
+	} else {
+		isInt = false;
+		isHalfInt = false;
+	}
+	if ( isInt && x > 0.6 ) {
+		// Calculate Q via finite sum:
+		invert = !invert;
+		evalMethod = 0;
+	} else if ( isHalfInt && x > 0.2 ) {
+		// Calculate Q via finite sum for half integer a:
+		invert = !invert;
+		evalMethod = 1;
+	} else if ( x < SQRT_EPSILON && a > 1.0 ) {
+		evalMethod = 6;
+	} else if ( ( x > 1000.0 ) && ( ( a < x ) || ( abs( a - 50.0 ) / x < 1.0 ) ) ) {
+		// Calculate Q via asymptotic approximation:
+		invert = !invert;
+		evalMethod = 7;
+	} else if ( x < 0.5 ) {
+		// Changeover criterion chosen to give a changeover at Q ~ 0.33:
+		if ( -0.4 / ln( x ) < a ) {
+			evalMethod = 2;
+		} else {
+			evalMethod = 3;
+		}
+	} else if ( x < 1.1 ) {
+		// Changeover here occurs when P ~ 0.75 or Q ~ 0.25:
+		if ( x * 0.75 < a ) {
+			evalMethod = 2;
+		} else {
+			evalMethod = 3;
+		}
+	} else {
+		// Begin by testing whether we're in the "bad" zone where the result will be near 0.5 and the usual series and continued fractions are slow to converge:
+		useTemme = false;
+		if ( regularized && a > 20 ) {
+			sigma = abs( (x-a)/a );
+			if ( a > 200 ) {
+				// Limit chosen so that we use Temme's expansion only if the result would be larger than about 10^-6. Below that the regular series and continued fractions converge OK, and if we use Temme's method we get increasing errors from the dominant erfc term as it's (inexact) argument increases in magnitude.
+				if ( 20 / a > sigma * sigma ) {
+					useTemme = true;
+				}
+			} else if ( sigma < 0.4 ) {
+				useTemme = true;
+			}
+		}
+		if ( useTemme ) {
+			evalMethod = 5;
+		}
+		// Regular case where the result will not be too close to 0.5: Changeover occurs at P ~ Q ~ 0.5. Note that series computation of P is about x2 faster than continued fraction calculation of Q, so try and use the CF only when really necessary, especially for small x.
+		else if ( x - ( 1.0 / (3.0 * x) ) < a ) {
+			evalMethod = 2;
+		} else {
+			evalMethod = 4;
+			invert = !invert;
+		}
+	}
+
+	/* eslint-disable default-case */
+	switch ( evalMethod ) {
+	case 0:
+		result = finiteGammaQ( a, x );
+		if ( regularized === false ) {
+			result *= gamma( a );
+		}
+		break;
+	case 1:
+		result = finiteHalfGammaQ( a, x );
+		if ( regularized === false ) {
+			result *= gamma( a );
+		}
+		break;
+	case 2:
+		// Compute P:
+		result = ( regularized ) ?
+			regularisedGammaPrefix( a, x ) :
+			fullIGammaPrefix( a, x );
+		if ( result !== 0.0 ) {
+			// If the result will be inverted, we can potentially save computation by initializing the series sum closer to the final result. This reduces the number of iterations needed in the series evaluation. However, care must be taken to avoid spurious numeric overflows. In practice, the more expensive overflow checks below are typically bypassed for well-behaved input values.
+			initValue = 0.0;
+			optimisedInvert = false;
+			if ( invert ) {
+				initValue = ( regularized ) ? 1.0 : gamma( a );
+				if (
+					regularized ||
+					result >= 1.0 ||
+					FLOAT64_MAX * result > initValue
+				) {
+					initValue /= result;
+					if (
+						regularized ||
+						a < 1.0 ||
+						( FLOAT64_MAX / a > initValue )
+					) {
+						initValue *= -a;
+						optimisedInvert = true;
+					} else {
+						initValue = 0.0;
+					}
+				} else {
+					initValue = 0.0;
+				}
+			}
+			result *= lowerGammaSeries( a, x, initValue ) / a;
+			if ( optimisedInvert ) {
+				invert = false;
+				result = -result;
+			}
+		}
+		break;
+	case 3:
+		// Compute Q:
+		invert = !invert;
+		res = tgammaSmallUpperPart( a, x, invert );
+		result = res[ 0 ];
+		g = res[ 1 ];
+		invert = false;
+		if ( regularized ) {
+			result /= g;
+		}
+		break;
+	case 4:
+		// Compute Q:
+		result = ( regularized ) ?
+			regularisedGammaPrefix( a, x ) :
+			fullIGammaPrefix( a, x );
+		if ( result !== 0 ) {
+			result *= upperGammaFraction( a, x );
+		}
+		break;
+	case 5:
+		result = igammaTemmeLarge( a, x );
+		if ( x >= a ) {
+			invert = !invert;
+		}
+		break;
+	case 6:
+		// Since `x` is so small that P is necessarily very small too, use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
+		result = ( regularized ) ?
+			pow( x, a ) / gamma( a + 1.0 ) :
+			pow( x, a ) / a;
+		result *= 1.0 - ( a * x / ( a + 1.0 ) );
+		break;
+	case 7:
+		// `x` is large, so compute Q:
+		result = ( regularized ) ?
+			regularisedGammaPrefix( a, x ) :
+			fullIGammaPrefix( a, x );
+		result /= x;
+		if ( result !== 0.0 ) {
+			result *= tgammaILargeX( a, x );
+		}
+		break;
+	}
+	if ( regularized && result > 1.0 ) {
+		result = 1.0;
+	}
+	if ( invert ) {
+		gam = ( regularized ) ? 1.0 : gamma( a );
+		result = gam - result;
+	}
+	return result;
+}
+
+
+// EXPORTS //
+
+module.exports = igammaFinal;
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/igamma_temme_large.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/igamma_temme_large.js
index 96a6db6e7d1f..e8640c5dd124 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/igamma_temme_large.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/igamma_temme_large.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_gamma_series.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_gamma_series.js
index c73d37112e9a..2f2680fca9d3 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_gamma_series.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_gamma_series.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
@@ -46,7 +46,7 @@ var lowerIncompleteGammaSeries = require( './lower_incomplete_gamma_series.js' )
 * ## Method
 *
 * -   Multiply result by `((z^a) * (e^-z) / a)` to get the full lower incomplete integral.
-* -   Divide by `tgamma(a)` to get the normalized value.
+* -   Divide by `gamma(a)` to get the normalized value.
 *
 * @private
 * @param {number} a - function parameter
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_incomplete_gamma_series.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_incomplete_gamma_series.js
index f75f3d77d11e..45813d48c616 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_incomplete_gamma_series.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/lower_incomplete_gamma_series.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
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 4f185c76f383..81ec1d9e96a9 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
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_62_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006-7, 2013-14.
@@ -37,32 +37,21 @@
 // MODULES //
 
 var gammaln = require( '@stdlib/math/base/special/gammaln' );
-var floor = require( '@stdlib/math/base/special/floor' );
-var gamma = require( '@stdlib/math/base/special/gamma' );
-var abs = require( '@stdlib/math/base/special/abs' );
 var exp = require( '@stdlib/math/base/special/exp' );
-var pow = require( '@stdlib/math/base/special/pow' );
 var ln = require( '@stdlib/math/base/special/ln' );
-var SQRT_EPSILON = require( '@stdlib/constants/float64/sqrt-eps' );
-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_NTH_FACTORIAL = require( '@stdlib/constants/float64/max-nth-factorial' );
-var finiteGammaQ = require( './finite_gamma_q.js' );
-var finiteHalfGammaQ = require( './finite_half_gamma_q.js' );
-var fullIGammaPrefix = require( './full_igamma_prefix.js' );
-var igammaTemmeLarge = require( './igamma_temme_large.js' );
+var igammaFinal = require( './igamma_final.js' );
 var lowerGammaSeries = require( './lower_gamma_series.js' );
-var regularisedGammaPrefix = require( './regularised_gamma_prefix.js' );
-var tgammaSmallUpperPart = require( './tgamma_small_upper_part.js' );
 var upperGammaFraction = require( './upper_gamma_fraction.js' );
 
 
 // MAIN //
 
 /**
-* Computes the regularized incomplete gamma function.
+* Computes the incomplete gamma function.
 *
 * ## Notes
 *
@@ -74,58 +63,67 @@ var upperGammaFraction = require( './upper_gamma_fraction.js' );
 * @param {boolean} [regularized=true] - boolean indicating if the function should evaluate the regularized or non-regularized incomplete gamma functions
 * @param {boolean} [upper=false] - boolean indicating if the function should return the upper tail of the incomplete gamma function
 * @returns {number} function value
+*
+* @example
+* var y = gammainc( 6.0, 2.0 );
+* // returns ~0.9826
+*
+* @example
+* var y = gammainc( 1.0, 2.0, true, true );
+* // returns ~0.7358
+*
+* @example
+* var y = gammainc( 7.0, 5.0 );
+* // returns ~0.8270
+*
+* @example
+* var y = gammainc( 7.0, 5.0, false );
+* // returns ~19.8482
+*
+* @example
+* var y = gammainc( NaN, 2.0 );
+* // returns NaN
+*
+* @example
+* var y = gammainc( 6.0, NaN );
+* // returns NaN
 */
 function gammainc( x, a, regularized, upper ) {
-	var optimisedInvert;
 	var normalized;
-	var evalMethod;
 	var initValue;
-	var isHalfInt;
-	var useTemme;
-	var isSmallA;
 	var invert;
 	var result;
-	var isInt;
-	var sigma;
-	var gam;
-	var res;
-	var fa;
-	var g;
 
 	if ( x < 0.0 || a <= 0.0 ) {
 		return NaN;
 	}
 	normalized = ( regularized === void 0 ) ? true : regularized;
 	invert = upper;
-	result = 0.0;
 	if ( a >= FLOAT64_MAX_NTH_FACTORIAL && !normalized ) {
 		if ( invert && ( a * 4.0 < x ) ) {
-			// This is method 4 below, done in logs:
+			// This is method 4 in `igammaFinal`, done in logs:
 			result = ( a * ln(x) ) - x;
 			result += ln( upperGammaFraction( a, x ) );
-		}
-		else if ( !invert && ( a > 4.0 * x ) ) {
-			// This is method 2 below, done in logs:
+		} else if ( !invert && ( a > 4.0 * x ) ) {
+			// This is method 2 in `igammaFinal`, done in logs:
 			result = ( a * ln(x) ) - x;
-			initValue = 0;
+			initValue = 0.0;
 			result += ln( lowerGammaSeries( a, x, initValue ) / a );
-		}
-		else {
-			result = gammainc( a, x, true, invert );
+		} else {
+			result = igammaFinal( x, a, true, invert );
 			if ( result === 0.0 ) {
 				if ( invert ) {
-					// Try http://functions.wolfram.com/06.06.06.0039.01
+					// Try http://functions.wolfram.com/06.06.06.0039.01:
 					result = 1.0 + ( 1.0 / (12.0*a) ) + ( 1.0 / (288.0*a*a) );
 					result = ln( result ) - a + ( ( a-0.5 ) * ln(a) );
 					result += ln( SQRT_TWO_PI );
 				} else {
-					// This is method 2 below, done in logs, we're really outside the range of this method, but since the result is almost certainly infinite, we should probably be OK:
+					// This is method 2 in `igammaFinal`, done in logs, we're really outside the range of this method, but since the result is almost certainly infinite, we should probably be OK:
 					result = ( a * ln( x ) ) - x;
 					initValue = 0.0;
 					result += ln( lowerGammaSeries( a, x, initValue ) / a );
 				}
-			}
-			else {
+			} else {
 				result = ln( result ) + gammaln( a );
 			}
 		}
@@ -134,164 +132,8 @@ function gammainc( x, a, regularized, upper ) {
 		}
 		return exp( result );
 	}
-	isSmallA = ( a < 30 ) && ( a <= x + 1.0 ) && ( x < MAX_LN );
-	if ( isSmallA ) {
-		fa = floor( a );
-		isInt = ( fa === a );
-		isHalfInt = ( isInt ) ? false : ( abs( fa - a ) === 0.5 );
-	} else {
-		isInt = isHalfInt = false;
-	}
-	if ( isInt && x > 0.6 ) {
-		// Calculate Q via finite sum:
-		invert = !invert;
-		evalMethod = 0;
-	}
-	else if ( isHalfInt && x > 0.2 ) {
-		// Calculate Q via finite sum for half integer a:
-		invert = !invert;
-		evalMethod = 1;
-	}
-	else if ( x < SQRT_EPSILON && a > 1.0 ) {
-		evalMethod = 6;
-	}
-	else if ( x < 0.5 ) {
-		// Changeover criterion chosen to give a changeover at Q ~ 0.33:
-		if ( -0.4 / ln( x ) < a ) {
-			evalMethod = 2;
-		} else {
-			evalMethod = 3;
-		}
-	}
-	else if ( x < 1.1 ) {
-		// Changeover here occurs when P ~ 0.75 or Q ~ 0.25:
-		if ( x * 0.75 < a ) {
-			evalMethod = 2;
-		} else {
-			evalMethod = 3;
-		}
-	}
-	else {
-		// Begin by testing whether we're in the "bad" zone where the result will be near 0.5 and the usual series and continued fractions are slow to converge:
-		useTemme = false;
-		if ( normalized && a > 20 ) {
-			sigma = abs( (x-a)/a );
-			if ( a > 200 ) {
-				// Limit chosen so that we use Temme's expansion only if the result would be larger than about 10^-6. Below that the regular series and continued fractions converge OK, and if we use Temme's method we get increasing errors from the dominant erfc term as it's (inexact) argument increases in magnitude.
-				if ( 20 / a > sigma * sigma ) {
-					useTemme = true;
-				}
-			} else if ( sigma < 0.4 ) {
-				useTemme = true;
-			}
-		}
-		if ( useTemme ) {
-			evalMethod = 5;
-		}
-		// Regular case where the result will not be too close to 0.5: Changeover occurs at P ~ Q ~ 0.5. Note that series computation of P is about x2 faster than continued fraction calculation of Q, so try and use the CF only when really necessary, especially for small x.
-		else if ( x - ( 1.0 / (3.0 * x) ) < a ) {
-			evalMethod = 2;
-		} else {
-			evalMethod = 4;
-			invert = !invert;
-		}
-	}
-
-	/* eslint-disable default-case */
-	switch ( evalMethod ) {
-	case 0:
-		result = finiteGammaQ( a, x );
-		if (normalized === false ) {
-			result *= gamma( a );
-		}
-		break;
-	case 1:
-		result = finiteHalfGammaQ( a, x );
-		if ( normalized === false ) {
-			result *= gamma( a );
-		}
-		break;
-	case 2:
-		// Compute P:
-		result = ( normalized ) ?
-			regularisedGammaPrefix( a, x ) :
-			fullIGammaPrefix( a, x );
-		if ( result !== 0.0 ) {
-			initValue = 0.0;
-			optimisedInvert = false;
-			if ( invert ) {
-				initValue = ( normalized ) ? 1.0 : gamma( a );
-				if (
-					normalized ||
-					result >= 1.0 ||
-					FLOAT64_MAX * result > initValue
-				) {
-					initValue /= result;
-					if (
-						normalized ||
-						a < 1.0 ||
-						( FLOAT64_MAX / a > initValue )
-					) {
-						initValue *= -a;
-						optimisedInvert = true;
-					}
-					else {
-						initValue = 0.0;
-					}
-				}
-				else {
-					initValue = 0.0;
-				}
-			}
-		}
-		result *= lowerGammaSeries( a, x, initValue ) / a;
-		if ( optimisedInvert ) {
-			invert = false;
-			result = -result;
-		}
-		break;
-	case 3:
-		// Compute Q:
-		invert = !invert;
-		res = tgammaSmallUpperPart( a, x, invert );
-		result = res[ 0 ];
-		g = res[ 1 ];
-		invert = false;
-		if ( normalized ) {
-			result /= g;
-		}
-		break;
-	case 4:
-		// Compute Q:
-		result = ( normalized ) ?
-			regularisedGammaPrefix( a, x ) :
-			fullIGammaPrefix( a, x );
-		if ( result !== 0 ) {
-			result *= upperGammaFraction( a, x );
-		}
-		break;
-	case 5:
-		result = igammaTemmeLarge( a, x );
-		if ( x >= a ) {
-			invert = !invert;
-		}
-		break;
-	case 6:
-		// Since x is so small that P is necessarily very small too, use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
-		result = ( normalized ) ?
-			pow(x, a) / gamma( a + 1.0 ) :
-			pow( x, a ) / a;
-		result *= 1.0 - ( a * x / ( a + 1.0 ) );
-		break;
-	}
-	if ( normalized && result > 1.0 ) {
-		result = 1.0;
-	}
-	if ( invert ) {
-		gam = ( normalized ) ? 1.0 : gamma( a );
-		result = gam - result;
-	}
-	return result;
+	// If no special handling is required then we proceed as normal:
+	return igammaFinal( x, a, normalized, invert );
 }
 
 
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/regularised_gamma_prefix.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/regularised_gamma_prefix.js
index 3df7edaa2554..fc790c6415b5 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/regularised_gamma_prefix.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/regularised_gamma_prefix.js
@@ -39,7 +39,7 @@
 var lanczosSumExpGScaled = require( '@stdlib/math/base/special/gamma-lanczos-sum-expg-scaled' );
 var gammaln = require( '@stdlib/math/base/special/gammaln' );
 var gamma = require( '@stdlib/math/base/special/gamma' );
-var log1p = require( '@stdlib/math/base/special/log1p' );
+var log1pmx = require( '@stdlib/math/base/special/log1pmx' );
 var sqrt = require( '@stdlib/math/base/special/sqrt' );
 var abs = require( '@stdlib/math/base/special/abs' );
 var exp = require( '@stdlib/math/base/special/exp' );
@@ -47,6 +47,7 @@ var pow = require( '@stdlib/math/base/special/pow' );
 var max = require( '@stdlib/math/base/special/max' );
 var min = require( '@stdlib/math/base/special/min' );
 var ln = require( '@stdlib/math/base/special/ln' );
+var FLOAT64_MAX = require( '@stdlib/constants/float64/max' );
 var MAX_LN = require( '@stdlib/constants/float64/max-ln' );
 var MIN_LN = require( '@stdlib/constants/float64/min-ln' );
 var G = require( '@stdlib/constants/float64/gamma-lanczos-g' );
@@ -75,22 +76,21 @@ function regularisedGammaPrefix( a, z ) {
 	agh = a + G - 0.5;
 	d = ( (z - a) - G + 0.5 ) / agh;
 	if ( a < 1.0 ) {
-		// Treat a < 1 as a special case because our Lanczos approximations are optimized against the factorials with a > 1, and for high precision types very small values of `a` can give rather erroneous results for gamma:
-		if ( z <= MIN_LN ) {
+		// Treat a < 1.0 as a special case because our Lanczos approximations are optimized against the factorials with a > 1.0, and for high precision types very small values of `a` can give rather erroneous results for gamma:
+		if ( z <= MIN_LN || a < 1.0 / FLOAT64_MAX ) {
 			// Use logs, so should be free of cancellation errors:
 			return exp( ( a * ln(z) ) - z - gammaln( a ) );
 		}
-		// No danger of overflow as gamma(a) < 1/a for small a, so direct calculation:
+		// No danger of overflow as `gamma(a) < 1/a` for small `a`, so direct calculation:
 		return pow( z, a ) * exp( -z ) / gamma( a );
 	}
 	if ( abs(d*d*a) <= 100.0 && a > 150.0 ) {
 		// Special case for large a and a ~ z:
-		prefix = ( a * ( log1p( d ) - d ) ) + ( z * ( 0.5-G ) / agh );
+		prefix = ( a * log1pmx( d ) ) + ( z * ( 0.5-G ) / agh );
 		prefix = exp( prefix );
-	}
-	else {
+	} else {
 		// General case. Direct computation is most accurate, but use various fallbacks for different parts of the problem domain:
-		alz = a * ln(z / agh);
+		alz = a * ln( z / agh );
 		amz = a - z;
 		if (
 			min(alz, amz) <= MIN_LN ||
@@ -104,8 +104,7 @@ function regularisedGammaPrefix( a, z ) {
 				// Compute square root of the result and then square it:
 				sq = pow( z / agh, a / 2.0 ) * exp( amz / 2.0 );
 				prefix = sq * sq;
-			}
-			else if (
+			} else if (
 				min(alz, amz)/4.0 > MIN_LN &&
 				max(alz, amz)/4.0 < MAX_LN &&
 				z > a
@@ -114,19 +113,15 @@ function regularisedGammaPrefix( a, z ) {
 				sq = pow( z / agh, a / 4.0 ) * exp( amz / 4.0 );
 				prefix = sq * sq;
 				prefix *= prefix;
-			}
-			else if (
+			} else if (
 				amza > MIN_LN &&
 				amza < MAX_LN
 			) {
 				prefix = pow( (z * exp(amza)) / agh, a );
-			}
-			else {
+			} else {
 				prefix = exp( alz + amz );
 			}
-		}
-		else
-		{
+		} else {
 			prefix = pow( z / agh, a ) * exp( amz );
 		}
 	}
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/small_gamma2_series.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/small_gamma2_series.js
index 54470236925d..f7bd4b67e1d2 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/small_gamma2_series.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/small_gamma2_series.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
@@ -33,7 +33,7 @@
 'use strict';
 
 /**
-* Series representation for upper fraction when `z` is small.
+* Series representation for the full upper fraction (Q) when `a` is very small.
 *
 * @private
 * @param {number} a - function parameter
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_i_large_x.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_i_large_x.js
new file mode 100644
index 000000000000..df665b99e13b
--- /dev/null
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_i_large_x.js
@@ -0,0 +1,63 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2025 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 https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+*
+* ```text
+* (C) Copyright John Maddock 2006.
+* (C) Copyright Paul A. Bristow 2007.
+*
+* 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 sumSeries = require( '@stdlib/math/base/tools/sum-series' );
+var tgammaILargeXSeries = require( './tgamma_i_large_x_series.js' );
+
+
+// MAIN //
+
+/**
+* Sums the elements of the series given by the full upper fraction (Q) when `x` is large.
+*
+* @private
+* @param {number} a - function parameter
+* @param {number} x - function parameter
+* @returns {number} sum of series
+*/
+function tgammaILargeX( a, x ) {
+	var result;
+	var s;
+
+	s = tgammaILargeXSeries( a, x );
+	result = sumSeries( s );
+	return result;
+}
+
+
+// EXPORTS //
+
+module.exports = tgammaILargeX;
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_i_large_x_series.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_i_large_x_series.js
new file mode 100644
index 000000000000..8d20922b8fd8
--- /dev/null
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_i_large_x_series.js
@@ -0,0 +1,68 @@
+/**
+* @license Apache-2.0
+*
+* Copyright (c) 2025 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 https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+*
+* ```text
+* (C) Copyright John Maddock 2006.
+* (C) Copyright Paul A. Bristow 2007.
+*
+* 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';
+
+// MAIN //
+
+/**
+* Creates a function to evaluate a series expansion of the full upper fraction (Q) when `x` is large.
+*
+* @private
+* @param {number} a1 - function parameter
+* @param {number} x1 - function parameter
+* @returns {Function} series function
+*/
+function tgammaILargeXSeries( a1, x1 ) {
+	var result = 1.0;
+	var a = a1;
+	var x = x1;
+	return next;
+
+	/**
+	* Calculate the next term of the series.
+	*
+	* @private
+	* @returns {number} series expansion term
+	*/
+	function next() {
+		var r = result;
+		result *= a / x;
+		a -= 1.0;
+		return r;
+	}
+}
+
+
+// EXPORTS //
+
+module.exports = tgammaILargeXSeries;
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_small_upper_part.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_small_upper_part.js
index bc675fa57878..be340ecb85b4 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_small_upper_part.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/tgamma_small_upper_part.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_gamma_fraction.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_gamma_fraction.js
index 4b88050cd610..af142add10cd 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_gamma_fraction.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_gamma_fraction.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_incomplete_gamma_fract.js b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_incomplete_gamma_fract.js
index 46d0e0f58565..51f15ccd9fb9 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_incomplete_gamma_fract.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/lib/upper_incomplete_gamma_fract.js
@@ -18,7 +18,7 @@
 *
 * ## Notice
 *
-* The original C++ code and copyright notice are from the [Boost library]{@link http://www.boost.org/doc/libs/1_37_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
+* The original C++ code and copyright notice are from the [Boost library]{@link https://www.boost.org/doc/libs/1_88_0/boost/math/special_functions/gamma.hpp}. The implementation has been modified for JavaScript.
 *
 * ```text
 * (C) Copyright John Maddock 2006.
@@ -55,7 +55,7 @@ function upperIncompleteGammaFract( a1, z1 ) {
 	* @returns {Array} series expansion terms
 	*/
 	function next() {
-		k += 1;
+		k += 1.0;
 		z += 2.0;
 		return [
 			k * (a - k),
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c b/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c
index b52f8b524d62..0a1db14da11a 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/src/main.c
@@ -348,7 +348,7 @@ static double upperGammaFraction( const double a, const double z ) {
 * Computes the next term in the series expansion of the lower incomplete gamma function.
 *
 * @param a          pointer to function parameter
-* @param x          pointer to function parameter
+* @param z          function parameter
 * @param result     pointer to result so far in series
 * @return           next term in series
 */
@@ -558,8 +558,7 @@ static double fullIGammaPrefix( const double a, const double z ) {
 		} else {
 			prefix = stdlib_base_exp( alz - z );
 		}
-	}
-	else {
+	} else {
 		if ( alz > MIN_LN ) {
 			prefix = stdlib_base_pow( z, a ) * stdlib_base_exp( -z );
 		} else if ( z/a < MAX_LN ) {
@@ -791,7 +790,7 @@ static double igammaTemmeLarge( const double a, const double x ) {
 * @param upper              boolean indicating if the function should return the upper tail of the incomplete gamma function
 * @return                   function value
 */
-static double igammaFinal( const double a, const double x, const bool regularized, const bool upper ) {
+static double igammaFinal( const double x, const double a, const bool regularized, const bool upper ) {
 	bool optimisedInvert;
 	int32_t evalMethod;
 	double initValue;
@@ -829,7 +828,7 @@ static double igammaFinal( const double a, const double x, const bool regularize
 	} else if ( x < SQRT_EPS && a > 1.0 ) {
 		evalMethod = 6;
 	} else if ( ( x > 1000.0 ) && ( ( a < x ) || ( stdlib_base_abs( a - 50.0 ) / x < 1.0 ) ) ) {
-		// calculate Q via asymptotic approximation:
+		// Calculate Q via asymptotic approximation:
 		invert = !invert;
 		evalMethod = 7;
 	} else if ( x < 0.5 ) {
@@ -1016,7 +1015,7 @@ double stdlib_base_gammainc( const double x, const double a, const bool regulari
 			initValue = 0.0;
 			result += stdlib_base_ln( lowerGammaSeries( a, x, initValue ) / a );
 		} else {
-			result = igammaFinal( a, x, true, invert );
+			result = igammaFinal( x, a, true, invert );
 			if ( result == 0.0 ) {
 				if ( invert ) {
 					// Try http://functions.wolfram.com/06.06.06.0039.01:
@@ -1039,5 +1038,5 @@ double stdlib_base_gammainc( const double x, const double a, const bool regulari
 		return stdlib_base_exp( result );
 	}
 	// If no special handling is required then we proceed as normal:
-	return igammaFinal( a, x, regularized, invert );
+	return igammaFinal( x, a, regularized, invert );
 }
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.js b/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.js
index b7fce2ab6a6c..1af976abb3b3 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.js
@@ -23,6 +23,7 @@
 var tape = require( 'tape' );
 var isSameValue = require( '@stdlib/assert/is-same-value' );
 var ulpdiff = require( '@stdlib/number/float64/base/ulp-difference' );
+var PINF = require( '@stdlib/constants/float64/pinf' );
 var gammainc = require( './../lib' );
 
 
@@ -32,6 +33,7 @@ var small = require( './fixtures/cpp/small.json' );
 var medium = require( './fixtures/cpp/medium.json' );
 var largeXSmallS = require( './fixtures/cpp/large_x_small_s.json' );
 var largeXMediumS = require( './fixtures/cpp/large_x_medium_s.json' );
+var largeXLargeS = require( './fixtures/cpp/large_x_large_s.json' );
 
 
 // TESTS //
@@ -369,3 +371,96 @@ tape( 'the function correctly evaluates the unregularized upper incomplete gamma
 	}
 	t.end();
 });
+
+tape( 'the function correctly evaluates the lower incomplete gamma function for large `x` and large `s`', function test( t ) {
+	var expected;
+	var actual;
+	var x;
+	var s;
+	var i;
+
+	expected = largeXLargeS.lower_regularized;
+	x = largeXLargeS.x;
+	s = largeXLargeS.s;
+
+	for ( i = 0; i < x.length; i++ ) {
+		actual = gammainc( x[ i ], s[ i ], true, false );
+
+		if ( expected[ i ] === 'PINF' ) {
+			t.strictEqual( isSameValue( actual, PINF ), true, 'returns expected value' );
+			continue;
+		}
+		// NOTE: The difference is high because some of the expected values are very small.
+		t.strictEqual( ulpdiff( actual, expected[ i ] ) <= 1370, true, 'returns expected value within 1370 ulp' );
+	}
+	t.end();
+});
+
+tape( 'the function correctly evaluates the upper incomplete gamma function for large `x` and large `s`', function test( t ) {
+	var expected;
+	var actual;
+	var x;
+	var s;
+	var i;
+
+	expected = largeXLargeS.upper_regularized;
+	x = largeXLargeS.x;
+	s = largeXLargeS.s;
+
+	for ( i = 0; i < x.length; i++ ) {
+		actual = gammainc( x[ i ], s[ i ], true, true );
+		if ( expected[ i ] === 'PINF' ) {
+			t.strictEqual( isSameValue( actual, PINF ), true, 'returns expected value' );
+			continue;
+		}
+		// NOTE: The difference is high because some of the expected values are very large.
+		t.strictEqual( ulpdiff( actual, expected[ i ] ) <= 658, true, 'returns expected value within 658 ulp' );
+	}
+	t.end();
+});
+
+tape( 'the function correctly evaluates the unregularized lower incomplete gamma function for large `x` and large `s`', function test( t ) {
+	var expected;
+	var actual;
+	var x;
+	var s;
+	var i;
+
+	expected = largeXLargeS.lower_unregularized;
+	x = largeXLargeS.x;
+	s = largeXLargeS.s;
+
+	for ( i = 0; i < x.length; i++ ) {
+		actual = gammainc( x[ i ], s[ i ], false, false );
+		if ( expected[ i ] === 'PINF' ) {
+			t.strictEqual( isSameValue( actual, PINF ), true, 'returns expected value' );
+			continue;
+		}
+		// NOTE: The difference is high because some of the expected values are very large and the compiler optimizations have been disabled.
+		t.strictEqual( ulpdiff( actual, expected[ i ] ) <= 947, true, 'returns expected value within 947 ulp' );
+	}
+	t.end();
+});
+
+tape( 'the function correctly evaluates the unregularized upper incomplete gamma function for large `x` and large `s`', function test( t ) {
+	var expected;
+	var actual;
+	var x;
+	var s;
+	var i;
+
+	expected = largeXLargeS.upper_unregularized;
+	x = largeXLargeS.x;
+	s = largeXLargeS.s;
+
+	for ( i = 0; i < x.length; i++ ) {
+		actual = gammainc( x[ i ], s[ i ], false, true );
+		if ( expected[ i ] === 'PINF' ) {
+			t.strictEqual( isSameValue( actual, PINF ), true, 'returns expected value' );
+			continue;
+		}
+		// NOTE: The difference is high because some of the expected values are very large and the compiler optimizations have been disabled.
+		t.strictEqual( ulpdiff( actual, expected[ i ] ) <= 1544, true, 'returns expected value within 1544 ulp' );
+	}
+	t.end();
+});
diff --git a/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.native.js
index f9af56aede7e..9b0c32af0192 100644
--- a/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.native.js
+++ b/lib/node_modules/@stdlib/math/base/special/gammainc/test/test.native.js
@@ -42,8 +42,6 @@ var small = require( './fixtures/cpp/small.json' );
 var medium = require( './fixtures/cpp/medium.json' );
 var largeXSmallS = require( './fixtures/cpp/large_x_small_s.json' );
 var largeXMediumS = require( './fixtures/cpp/large_x_medium_s.json' );
-
-// TODO: Add this to `test.js` once the JS implementation matches the latest Boost C++ implementation.
 var largeXLargeS = require( './fixtures/cpp/large_x_large_s.json' );