From dd7d23f137eda8262338adfc12d7de4dc2be86ed Mon Sep 17 00:00:00 2001 From: 0PrashantYadav0 Date: Fri, 23 Oct 2026 23:05:49 +0530 Subject: [PATCH] fix: avoid premature underflow in `math/base/special/gamma` --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown_pkg_readmes status: na - task: lint_markdown_docs status: na - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: passed - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: missing_dependencies - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../math/base/special/factorial/test/test.js | 4 +-- .../special/factorial/test/test.native.js | 4 +-- .../math/base/special/gamma/lib/main.js | 10 +++++- .../math/base/special/gamma/manifest.json | 3 ++ .../math/base/special/gamma/src/main.c | 9 ++++- .../math/base/special/gamma/test/test.js | 35 +++++++++++++++++-- .../base/special/gamma/test/test.native.js | 35 +++++++++++++++++-- 7 files changed, 90 insertions(+), 10 deletions(-) diff --git a/lib/node_modules/@stdlib/math/base/special/factorial/test/test.js b/lib/node_modules/@stdlib/math/base/special/factorial/test/test.js index 40949318d0f1..bcc4cf5d5b04 100644 --- a/lib/node_modules/@stdlib/math/base/special/factorial/test/test.js +++ b/lib/node_modules/@stdlib/math/base/special/factorial/test/test.js @@ -87,12 +87,12 @@ tape( 'if `x > 170.6144...`, the function returns positive infinity', function t t.end(); }); -tape( 'if `x < -171.56749...`, the function returns zero', function test( t ) { +tape( 'if `x < -185`, the function returns zero', function test( t ) { var values; var v; var i; - values = incrspace( -171.57, -1000.0, -10.1234 ); + values = incrspace( -185.01, -1000.0, -10.1234 ); for ( i = 0; i < values.length; i++ ) { v = factorial( values[ i ] ); t.strictEqual( v, 0.0, 'returns expected value when provided ' + values[ i ] ); diff --git a/lib/node_modules/@stdlib/math/base/special/factorial/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/factorial/test/test.native.js index c626c6693b6f..fbe8a6ece360 100644 --- a/lib/node_modules/@stdlib/math/base/special/factorial/test/test.native.js +++ b/lib/node_modules/@stdlib/math/base/special/factorial/test/test.native.js @@ -96,12 +96,12 @@ tape( 'if `x > 170.6144...`, the function returns positive infinity', opts, func t.end(); }); -tape( 'if `x < -171.56749...`, the function returns zero', opts, function test( t ) { +tape( 'if `x < -185`, the function returns zero', opts, function test( t ) { var values; var v; var i; - values = incrspace( -171.57, -1000.0, -10.1234 ); + values = incrspace( -185.01, -1000.0, -10.1234 ); for ( i = 0; i < values.length; i++ ) { v = factorial( values[ i ] ); t.strictEqual( v, 0.0, 'returns expected value when provided ' + values[ i ] ); diff --git a/lib/node_modules/@stdlib/math/base/special/gamma/lib/main.js b/lib/node_modules/@stdlib/math/base/special/gamma/lib/main.js index 66e5ba73f8cd..88d094a4601c 100644 --- a/lib/node_modules/@stdlib/math/base/special/gamma/lib/main.js +++ b/lib/node_modules/@stdlib/math/base/special/gamma/lib/main.js @@ -39,10 +39,12 @@ var isInteger = require( '@stdlib/math/base/assert/is-integer' ); var isNegativeZero = require( '@stdlib/math/base/assert/is-negative-zero' ); var abs = require( '@stdlib/math/base/special/abs' ); var floor = require( '@stdlib/math/base/special/floor' ); +var pow = require( '@stdlib/math/base/special/pow' ); var sin = require( '@stdlib/math/base/special/sin' ); var PINF = require( '@stdlib/constants/float64/pinf' ); var NINF = require( '@stdlib/constants/float64/ninf' ); var PI = require( '@stdlib/constants/float64/pi' ); +var SQRT_PI = require( '@stdlib/constants/float64/sqrt-pi' ); var stirlingApprox = require( './stirling_approximation.js' ); var smallApprox = require( './small_approximation.js' ); var rateval = require( './rational_pq.js' ); @@ -124,7 +126,7 @@ function gamma( x ) { if ( x > 171.61447887182298 ) { return PINF; } - if ( x < -170.5674972726612 ) { + if ( x < -184.0 ) { return 0.0; } q = abs( x ); @@ -146,6 +148,12 @@ function gamma( x ) { z = q - p; } z = q * sin( PI * z ); + if ( x < -170.5674972726612 ) { + // Avoid premature underflow by evaluating `gamma(q)` using the duplication formula `gamma(q) = gamma(q/2) * gamma(q/2 + 1/2) * 2^(q-1) / sqrt(pi)`: + p = 0.5 * q; + z = abs( z ) * stirlingApprox( p ) * stirlingApprox( p+0.5 ); + return sign * ( ( PI*SQRT_PI ) / z ) * pow( 2.0, 1.0-q ); + } return sign * PI / ( abs(z)*stirlingApprox(q) ); } // Reduce `x`... diff --git a/lib/node_modules/@stdlib/math/base/special/gamma/manifest.json b/lib/node_modules/@stdlib/math/base/special/gamma/manifest.json index 95be3c49993e..ffc3d0af9f6d 100644 --- a/lib/node_modules/@stdlib/math/base/special/gamma/manifest.json +++ b/lib/node_modules/@stdlib/math/base/special/gamma/manifest.json @@ -46,6 +46,7 @@ "@stdlib/constants/float64/pinf", "@stdlib/constants/float64/ninf", "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/sqrt-pi", "@stdlib/constants/float64/eulergamma", "@stdlib/constants/float64/sqrt-two-pi", "@stdlib/math/base/special/pow", @@ -73,6 +74,7 @@ "@stdlib/constants/float64/pinf", "@stdlib/constants/float64/ninf", "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/sqrt-pi", "@stdlib/constants/float64/eulergamma", "@stdlib/constants/float64/sqrt-two-pi", "@stdlib/math/base/special/pow", @@ -100,6 +102,7 @@ "@stdlib/constants/float64/pinf", "@stdlib/constants/float64/ninf", "@stdlib/constants/float64/pi", + "@stdlib/constants/float64/sqrt-pi", "@stdlib/constants/float64/eulergamma", "@stdlib/constants/float64/sqrt-two-pi", "@stdlib/math/base/special/pow", diff --git a/lib/node_modules/@stdlib/math/base/special/gamma/src/main.c b/lib/node_modules/@stdlib/math/base/special/gamma/src/main.c index 29b600d65a55..596e90a6967f 100644 --- a/lib/node_modules/@stdlib/math/base/special/gamma/src/main.c +++ b/lib/node_modules/@stdlib/math/base/special/gamma/src/main.c @@ -40,6 +40,7 @@ #include "stdlib/constants/float64/pinf.h" #include "stdlib/constants/float64/ninf.h" #include "stdlib/constants/float64/pi.h" +#include "stdlib/constants/float64/sqrt_pi.h" #include "stdlib/constants/float64/eulergamma.h" #include "stdlib/constants/float64/sqrt_two_pi.h" #include "stdlib/constants/float64/nan.h" @@ -199,7 +200,7 @@ double stdlib_base_gamma( const double x ) { if ( x > 171.61447887182298 ) { return STDLIB_CONSTANT_FLOAT64_PINF; } - if ( x < -170.5674972726612 ) { + if ( x < -184.0 ) { return 0.0; } q = stdlib_base_abs( x ); @@ -222,6 +223,12 @@ double stdlib_base_gamma( const double x ) { z = q - p; } z = q * stdlib_base_sin( STDLIB_CONSTANT_FLOAT64_PI * z ); + if ( x < -170.5674972726612 ) { + // Avoid premature underflow by evaluating `gamma(q)` using the duplication formula `gamma(q) = gamma(q/2) * gamma(q/2 + 1/2) * 2^(q-1) / sqrt(pi)`: + p = 0.5 * q; + z = stdlib_base_abs( z ) * stirlingApprox( p ) * stirlingApprox( p + 0.5 ); + return sign * ( ( STDLIB_CONSTANT_FLOAT64_PI*STDLIB_CONSTANT_FLOAT64_SQRT_PI ) / z ) * stdlib_base_pow( 2.0, 1.0-q ); + } return sign * STDLIB_CONSTANT_FLOAT64_PI / ( stdlib_base_abs( z ) * stirlingApprox( q ) ); } diff --git a/lib/node_modules/@stdlib/math/base/special/gamma/test/test.js b/lib/node_modules/@stdlib/math/base/special/gamma/test/test.js index 6509ae05d6f3..73e5001e7ad9 100644 --- a/lib/node_modules/@stdlib/math/base/special/gamma/test/test.js +++ b/lib/node_modules/@stdlib/math/base/special/gamma/test/test.js @@ -93,8 +93,8 @@ tape( 'if `x > 171.6144...`, the function returns positive infinity', function t t.end(); }); -tape( 'if `x < -170.56749...`, the function returns zero', function test( t ) { - var values = incrspace( -170.57, -1000.0, -10.1234 ); +tape( 'if `x < -184`, the function returns zero', function test( t ) { + var values = incrspace( -184.01, -1000.0, -10.1234 ); var v; var i; @@ -105,6 +105,37 @@ tape( 'if `x < -170.56749...`, the function returns zero', function test( t ) { t.end(); }); +tape( 'the function returns subnormal values for `-184 < x < -170.56749...`', function test( t ) { + var expected; + var x; + var v; + var i; + + x = [ + -170.6, + -171.5, + -172.5, + -175.25, + -177.125, + -180.99999999999997, + -183.00000000000003 + ]; + expected = [ + -2.0829008361502995e-308, + 1.9316265431712e-310, + -1.11978350329e-312, + 1.085373e-318, + 1.24e-322, + -9.6761e-319, + 3.0e-323 + ]; + for ( i = 0; i < x.length; i++ ) { + v = gamma( x[ i ] ); + t.strictEqual( isAlmostSameValue( v, expected[ i ], 1 ), true, 'returns expected value when provided ' + x[ i ] ); + } + t.end(); +}); + tape( 'the function evaluates the gamma function (positive integers)', function test( t ) { var v; var i; diff --git a/lib/node_modules/@stdlib/math/base/special/gamma/test/test.native.js b/lib/node_modules/@stdlib/math/base/special/gamma/test/test.native.js index 1c5cc61fe808..ed1a5335db79 100644 --- a/lib/node_modules/@stdlib/math/base/special/gamma/test/test.native.js +++ b/lib/node_modules/@stdlib/math/base/special/gamma/test/test.native.js @@ -102,8 +102,8 @@ tape( 'if `x > 171.6144...`, the function returns positive infinity', opts, func t.end(); }); -tape( 'if `x < -170.56749...`, the function returns zero', opts, function test( t ) { - var values = incrspace( -170.57, -1000.0, -10.1234 ); +tape( 'if `x < -184`, the function returns zero', opts, function test( t ) { + var values = incrspace( -184.01, -1000.0, -10.1234 ); var v; var i; @@ -114,6 +114,37 @@ tape( 'if `x < -170.56749...`, the function returns zero', opts, function test( t.end(); }); +tape( 'the function returns subnormal values for `-184 < x < -170.56749...`', opts, function test( t ) { + var expected; + var x; + var v; + var i; + + x = [ + -170.6, + -171.5, + -172.5, + -175.25, + -177.125, + -180.99999999999997, + -183.00000000000003 + ]; + expected = [ + -2.0829008361502995e-308, + 1.9316265431712e-310, + -1.11978350329e-312, + 1.085373e-318, + 1.24e-322, + -9.6761e-319, + 3.0e-323 + ]; + for ( i = 0; i < x.length; i++ ) { + v = gamma( x[ i ] ); + t.strictEqual( isAlmostSameValue( v, expected[ i ], 1 ), true, 'returns expected value when provided ' + x[ i ] ); + } + t.end(); +}); + tape( 'the function evaluates the gamma function (positive integers)', opts, function test( t ) { var v; var i;