Skip to content

refactor: improve the implementation of math/base/special/hyp2f1 #11325

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

Merged
kgryte merged 5 commits into from
Apr 10, 2026
Merged

refactor: improve the implementation of math/base/special/hyp2f1 #11325

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b359e62
chore: improve the implementation
officiallyanee Apr 7, 2026
c5fec0a
Merge remote-tracking branch 'upstream/develop' into improve/hyp2f1
stdlib-bot Apr 9, 2026
cc3d8a1
Merge remote-tracking branch 'upstream/develop' into improve/hyp2f1
stdlib-bot Apr 10, 2026
c9ec417
chore: update test to use ULP-based testing & move hyp2f1NegCEqualBC …
officiallyanee Apr 10, 2026
f39d34d
Apply suggestions from code review
kgryte Apr 10, 2026
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
8 changes: 6 additions & 2 deletions lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/hyt2f1.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@
* @param {number} loss - starting loss of significance
* @returns {Object} the function value and error
*/
function hyt2f1( a, b, c, x, loss ) {

Check warning on line 73 in lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/hyt2f1.js

View workflow job for this annotation

GitHub Actions / Lint Changed Files 

Function 'hyt2f1' has too many statements (120). Maximum allowed is 100
var negIntA;
var negIntB;
var qVal;
Expand Down Expand Up @@ -122,7 +122,7 @@
d = c - a - b;
id = round( d );

if ( x > 0.9 && !negIntA && !negIntB ) {
if ( x > 0.85 && !negIntA && !negIntB ) {
Copy link
Copy Markdown
Contributor Author

@officiallyanee officiallyanee Apr 7, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

0.85 was chosen arbitrarily through trial and error, this helped the case I mentioned in the description.

Copy link
Copy Markdown
Member

@kgryte kgryte Apr 9, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@officiallyanee Do you happen to know where/why the 0.9 was chosen in the first place?

Copy link
Copy Markdown
Contributor Author

@officiallyanee officiallyanee Apr 9, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Its not explicitly mentioned in the cephes/scipy code but I went through some resources and found this:
image
ref: https://link.springer.com/article/10.1007/s11075-016-0173-0
I think it also explains why the initial case was performing badly too since |a| and |b| were much larger than |c|
and |x| is also pretty close to 1 so the power series was performing poorly.
In the range 0.85<x<0.9, it was true for every value as long as we keep a, b, c same.
( this is assuming wolfram as a better source of comparision )

Copy link
Copy Markdown
Member

@kgryte kgryte Apr 10, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for digging up the reference!

if ( isInteger( d ) === false ) {
// Try the power series first:
y = hys2f1( a, b, c, x, err );
Expand Down Expand Up @@ -223,7 +223,11 @@
y = -y;
}
q = pow( s, id );
y = ( id > 0.0 ) ? y*q : y1*q;
if ( id > 0.0 ) {
y *= q;
} else {
y1 *= q;
}
Comment thread
kgryte marked this conversation as resolved.
y += y1;
}
return {
Expand Down
40 changes: 40 additions & 0 deletions lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/main.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
var gamma = require( '@stdlib/math/base/special/gamma' );
var pow = require( '@stdlib/math/base/special/pow' );
var abs = require( '@stdlib/math/base/special/abs' );
var max = require( '@stdlib/math/base/special/max' );
var isNonPositiveInteger = require( './isnonpositiveinteger.js' );
var isInteger = require( './isinteger.js' );
var hys2f1 = require( './hys2f1.js' );
Expand All @@ -51,6 +52,41 @@
var ETHRESH = config.ETHRESH;


// FUNCTIONS //

/**
* Evaluates 2F1(a, b; b; x) when `b = c` is a negative integer using AMS55 #15.4.2.
*
* @private
* @param {number} a - first parameter
* @param {number} b - second parameter (equals c, a non-positive integer)
* @param {number} x - argument
* @returns {number} function value, or NaN if precision is insufficient
*/
function hyp2f1NegCEqualBC( a, b, x ) {
var collectorMax;
var collector;
var sum;
var k;

if ( !( abs( b ) < 1.0e5 ) ) {
return NaN;
}
collectorMax = 1.0;
collector = 1.0;
sum = 1.0;
for ( k = 1; k <= -b; k++ ) {
collector *= ( a + k - 1.0 ) * x / k;
collectorMax = max( abs( collector ), collectorMax );
sum += collector;
}
if ( 1.0e-16 * ( 1.0 + ( collectorMax / abs( sum ) ) ) > 1.0e-7 ) {
return NaN;
}
return sum;
}


Comment thread
kgryte marked this conversation as resolved.
Outdated
// MAIN //

/**
Expand Down Expand Up @@ -90,7 +126,7 @@
* var v = hyp2f1( 1.0, NaN, 2.0, 0.5 );
* // returns NaN
*/
function hyp2f1( a, b, c, x ) {

Check warning on line 129 in lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/main.js

View workflow job for this annotation

GitHub Actions / Lint Changed Files 

Function 'hyp2f1' has too many statements (122). Maximum allowed is 100
var negIntCaOrCb;
var negIntC;
var negIntB;
Expand Down Expand Up @@ -156,6 +192,10 @@
if ( ax < 1.0 || x === -1.0 ) {
if ( b === c ) {
// 2F1(a,b;b;x) = (1-x)**(-a):
if ( negIntB ) {
// For negative integer b=c use the finite polynomial (AMS55 #15.4.2):
return hyp2f1NegCEqualBC( a, b, x );
}
return pow( s, -a );
}
if ( a === c ) {
Expand Down
2 changes: 1 addition & 1 deletion lib/node_modules/@stdlib/math/base/special/hyp2f1/test/fixtures/python/edge_cases3.json
View file
Open in desktop

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions lib/node_modules/@stdlib/math/base/special/hyp2f1/test/fixtures/python/edge_cases5.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
{"a": [1.5, 1.5, 2.5, 0.5, 1.5, 0.5, 2.0, 1.5, 2.5, 1.5], "b": [2.5, 2.5, 1.5, 2.5, 2.5, 1.5, 3.0, 3.5, 2.5, 3.5], "c": [3.0, 3.0, 3.0, 2.0, 2.0, 1.0, 4.0, 4.0, 4.0, 3.0], "x": [0.9, 0.95, 0.9, 0.92, 0.91, 0.9, 0.9, 0.88, 0.93, 0.91], "expected": [14.663397526525134, 30.949523656432504, 14.663397526525134, 6.507774205759818, 106.96827020610523, 7.033214388515233, 21.78942310292968, 13.41253423166626, 36.54504530314331, 88.86455455103395]}
1 change: 1 addition & 0 deletions lib/node_modules/@stdlib/math/base/special/hyp2f1/test/fixtures/python/edge_cases6.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
{"a": [2.0, 2.0, 3.0, 0.5, 1.5, 2.0, 4.0, 0.5], "b": [-1.0, -2.0, -2.0, -3.0, -3.0, -4.0, -2.0, -1.0], "c": [-1.0, -2.0, -2.0, -3.0, -3.0, -4.0, -2.0, -1.0], "x": [0.5, 0.5, 0.3, 0.7, -0.5, 0.6, 0.8, -0.9], "expected": [2.0, 2.75, 2.44, 1.6409375, 0.4453125, 4.792, 10.600000000000001, 0.55]}
2 changes: 1 addition & 1 deletion lib/node_modules/@stdlib/math/base/special/hyp2f1/test/fixtures/python/outliers.json
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1 +1 @@
{"a": [-1.8, 1.8, 1.8, 1.8, 1.8, 5.0, 10.0, 10.0, 2.0, 2.0, -1.8, 1.8, 10.0, 10.0, 0.5, -8.0, -0.1156544430516313, -27.911722328444608, -10.184355055477743, -5.0083982772911, -0.7429095709749414, -14.55029516886238, -5.428179853908536, -1.8190858888165087, -0.9534205693212994, -31.70680212982947, -0.18350917503602782, -200.19574300099214, -8.837706832399117, -3.2052905962805456, -13.37614633859299], "b": [7.4, -2.5, 1.0, 7.4, 7.4, 7.4, 7.4, 7.4, 3.0, 2.5, -2.5, 7.4, 1.0, 7.4, -269.5, 18.016500331508873, 12.722140489351698, 14.018736008860946, 6.5450429328481174, 171.09919289684245, 11.725525989335528, 0.46553367087834685, 26.782963689614604, 21.025202609254364, 10.024920055648606, 0.11445767340979862, 14.696533286274349, 1.9444493115386567, 0.19096384651567044, 0.3048596186043564, 1.337786902832312], "c": [20.4, 20.4, 20.4, -1.8, 20.4, 20.4, 20.4, 20.4, 5.0, -3.25, 20.4, -1.8, -1.8, -1.8, 1.5, 10.805295997850628, -0.6030985338257255, -1.0168630024136034, -0.44051751620383994, -1.4923865283013589, -7.299203351011158, -0.5349901418098844, -4.162402828198655, -0.6960141045675334, -0.23845367675425955, -0.5570625216964069, -6.338990250643928, -0.6933677849479691, -2.3808990396202425, -0.7039276707364683, -0.05015281580946018], "x": [-10.0, -10.0, -10.0, -10.0, -10.0, -10.0, -10.0, 0.95, 0.99, 0.999, -10.0, -1.01, -0.99, -0.99, 0.998001, 0.90875647507, -0.0005448214561780684, 0.6228483452847298, 0.6555896492302491, 0.7073764056574465, 0.11534401417791407, 0.6086682776938632, 0.10047548647104909, 0.2570211738060222, 0.4419851803780206, 0.020574313007314604, -0.9003519736816863, 0.24492811068225362, -0.27624718180517105, 0.8496565195400692, -0.0022109253176725296], "expected": [16.334204105729413, 5.883122144128727, 0.5820646473844754, 0.010147094307759376, 0.07291374137469146, 0.0011628419587190388, 8.042653908304492e-06, 566.0092271657466, 27.699347904322664, 2.183739328012162e+26, -0.5786871007671348, 0.10921062752340924, -0.8752587026378827, 0.4234752106080307, 0.053963052503373715, -3.566216341436626e-09, 0.998681836639102, 5.33299803869569, 1.1177709865124041, 3.491332086988064e+82, 1.132062691110089, -0.0006927158077474787, 98.41698380414908, -221.57018712798285, 43.53554136829992, 1.0540592696289446, 0.5371016318543966, 0.0028040224553278668, -2.4067466563131665, 0.028751519134363426, 0.18419839456905535]}
{"a": [-1.8, 1.8, 1.8, 1.8, 1.8, 5.0, 10.0, 10.0, 2.0, 2.0, -1.8, 1.8, 10.0, 10.0, 0.5, -8.0, -0.1156544430516313, -27.911722328444608, -10.184355055477743, -5.0083982772911, -0.7429095709749414, -14.55029516886238, -5.428179853908536, -1.8190858888165087, -0.9534205693212994, -31.70680212982947, -0.18350917503602782, -200.19574300099214, -8.837706832399117, -3.2052905962805456, -13.37614633859299], "b": [7.4, -2.5, 1.0, 7.4, 7.4, 7.4, 7.4, 7.4, 3.0, 2.5, -2.5, 7.4, 1.0, 7.4, -269.5, 18.016500331508873, 12.722140489351698, 14.018736008860946, 6.5450429328481174, 171.09919289684245, 11.725525989335528, 0.46553367087834685, 26.782963689614604, 21.025202609254364, 10.024920055648606, 0.11445767340979862, 14.696533286274349, 1.9444493115386567, 0.19096384651567044, 0.3048596186043564, 1.337786902832312], "c": [20.4, 20.4, 20.4, -1.8, 20.4, 20.4, 20.4, 20.4, 5.0, -3.25, 20.4, -1.8, -1.8, -1.8, 1.5, 10.805295997850628, -0.6030985338257255, -1.0168630024136034, -0.44051751620383994, -1.4923865283013589, -7.299203351011158, -0.5349901418098844, -4.162402828198655, -0.6960141045675334, -0.23845367675425955, -0.5570625216964069, -6.338990250643928, -0.6933677849479691, -2.3808990396202425, -0.7039276707364683, -0.05015281580946018], "x": [-10.0, -10.0, -10.0, -10.0, -10.0, -10.0, -10.0, 0.95, 0.99, 0.999, -10.0, -1.01, -0.99, -0.99, 0.998001, 0.90875647507, -0.0005448214561780684, 0.6228483452847298, 0.6555896492302491, 0.7073764056574465, 0.11534401417791407, 0.6086682776938632, 0.10047548647104909, 0.2570211738060222, 0.4419851803780206, 0.020574313007314604, -0.9003519736816863, 0.24492811068225362, -0.27624718180517105, 0.8496565195400692, -0.0022109253176725296], "expected": [16.334204105729413, 5.883122144128727, 0.5820646473844754, 0.010147094307759376, 0.07291374137469146, 0.0011628419587190388, 8.042653908304492e-06, 566.0092271658314, 27.699347904322664, 2.1837393280121743e+26, -0.5786871007671348, 0.10921062752340924, -0.8752587026378827, 0.4234752106080307, 0.053963052503373715, -3.566216341436626e-09, 0.998681836639102, 5.33299803869569, 1.1177709865124041, 3.491332086988064e+82, 1.132062691110089, -0.0006927158077474787, 98.41698380414908, -221.57018712798285, 43.53554136829992, 1.0540592696289446, 0.5371016318543966, 0.0028040224553278668, -2.4067466563131665, 0.028751519134363426, 0.18419839456905535]}
40 changes: 40 additions & 0 deletions lib/node_modules/@stdlib/math/base/special/hyp2f1/test/fixtures/python/runner.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,6 +256,46 @@ def main():
# Outliers:
gen_outliers("outliers.json")

# Edge Cases #5
# d = c - a - b is a negative integer, x > 0.85, a and b are positive integers.
pts = [
(1.5, 2.5, 3.0, 0.90),
(1.5, 2.5, 3.0, 0.95),
(2.5, 1.5, 3.0, 0.90),
(0.5, 2.5, 2.0, 0.92),
(1.5, 2.5, 2.0, 0.91),
(0.5, 1.5, 1.0, 0.90),
(2.0, 3.0, 4.0, 0.90),
(1.5, 3.5, 4.0, 0.88),
(2.5, 2.5, 4.0, 0.93),
(1.5, 3.5, 3.0, 0.91),
]
a, b, c, x = zip(*pts)
a = np.array(a)
b = np.array(b)
c = np.array(c)
x = np.array(x)
gen(a, b, c, x, "edge_cases5.json")

# Edge Cases #6
# b = c is a negative integer: 2F1(a, b; b; x) -> finite polynomial sum.
pts = [
(2.0, -1.0, -1.0, 0.5),
(2.0, -2.0, -2.0, 0.5),
(3.0, -2.0, -2.0, 0.3),
(0.5, -3.0, -3.0, 0.7),
(1.5, -3.0, -3.0, -0.5),
(2.0, -4.0, -4.0, 0.6),
(4.0, -2.0, -2.0, 0.8),
(0.5, -1.0, -1.0, -0.9),
]
a, b, c, x = zip(*pts)
a = np.array(a)
b = np.array(b)
c = np.array(c)
x = np.array(x)
gen(a, b, c, x, "edge_cases6.json")


if __name__ == "__main__":
main()
54 changes: 54 additions & 0 deletions lib/node_modules/@stdlib/math/base/special/hyp2f1/test/test.js
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,8 @@ var edgeCases1 = require( './fixtures/python/edge_cases1.json' );
var edgeCases2 = require( './fixtures/python/edge_cases2.json' );
var edgeCases3 = require( './fixtures/python/edge_cases3.json' );
var edgeCases4 = require( './fixtures/python/edge_cases4.json' );
var edgeCases5 = require( './fixtures/python/edge_cases5.json' );
var edgeCases6 = require( './fixtures/python/edge_cases6.json' );
var outliers = require( './fixtures/python/outliers.json' );


Expand Down Expand Up @@ -259,6 +261,58 @@ tape( 'the function correctly evaluates the hypergeometric function', function t
t.end();
});

tape( 'the function correctly evaluates the hypergeometric function', function test( t ) {
var expected;
var delta;
var tol;
var a;
var b;
var c;
var x;
var v;
var i;

a = edgeCases5.a;
b = edgeCases5.b;
c = edgeCases5.c;
x = edgeCases5.x;
expected = edgeCases5.expected;

for ( i = 0; i < x.length; i++ ) {
v = hyp2f1( a[ i ], b[ i ], c[ i ], x[ i ] );
delta = abs( v - expected[ i ] );
tol = 8.0 * EPS * abs( expected[ i ] );
Comment thread
kgryte marked this conversation as resolved.
Outdated
t.ok( delta <= tol, 'within tolerance. a: ' + a[ i ] + ' b: ' + b[ i ] + ' c: ' + c[ i ] + ' x: ' + x[ i ] + '. Value: ' + v + '. Expected: ' + expected[ i ] + '. Delta: ' + delta + '. Tolerance: ' + tol + '.' );
}
t.end();
});

tape( 'the function correctly evaluates the hypergeometric function', function test( t ) {
var expected;
var delta;
var tol;
var a;
var b;
var c;
var x;
var v;
var i;

a = edgeCases6.a;
b = edgeCases6.b;
c = edgeCases6.c;
x = edgeCases6.x;
expected = edgeCases6.expected;

for ( i = 0; i < x.length; i++ ) {
v = hyp2f1( a[ i ], b[ i ], c[ i ], x[ i ] );
delta = abs( v - expected[ i ] );
tol = EPS * abs( expected[ i ] );
t.ok( delta <= tol, 'within tolerance. a: ' + a[ i ] + ' b: ' + b[ i ] + ' c: ' + c[ i ] + ' x: ' + x[ i ] + '. Value: ' + v + '. Expected: ' + expected[ i ] + '. Delta: ' + delta + '. Tolerance: ' + tol + '.' );
}
t.end();
});

tape( 'the function correctly evaluates the hypergeometric function', function test( t ) {
var expected;
var delta;
Expand Down