Skip to content

fix: avoid premature underflow in math/base/special/gamma #16298

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
0PrashantYadav0 wants to merge 1 commit into
base: develop
Choose a base branch
from
+90 −10
Open

fix: avoid premature underflow in math/base/special/gamma #16298

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions lib/node_modules/@stdlib/math/base/special/factorial/test/test.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 ] );
Expand Down
4 changes: 2 additions & 2 deletions lib/node_modules/@stdlib/math/base/special/factorial/test/test.native.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 ] );
Expand Down
10 changes: 9 additions & 1 deletion lib/node_modules/@stdlib/math/base/special/gamma/lib/main.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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' );
Expand Down Expand Up @@ -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 );
Expand All @@ -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`...
Expand Down
3 changes: 3 additions & 0 deletions lib/node_modules/@stdlib/math/base/special/gamma/manifest.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down
9 changes: 8 additions & 1 deletion lib/node_modules/@stdlib/math/base/special/gamma/src/main.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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 );
Expand All @@ -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 ) );
}

Expand Down
35 changes: 33 additions & 2 deletions lib/node_modules/@stdlib/math/base/special/gamma/test/test.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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;
Expand Down
35 changes: 33 additions & 2 deletions lib/node_modules/@stdlib/math/base/special/gamma/test/test.native.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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;
Expand Down
Loading