Skip to content

Commit 2244f8e

Browse files
officiallyaneekgrytestdlib-bot
authored
refactor: improve the implementation of math/base/special/hyp2f1
PR-URL: #11325 Co-authored-by: Athan Reines <kgryte@gmail.com> Reviewed-by: Athan Reines <kgryte@gmail.com> Co-authored-by: stdlib-bot <noreply@stdlib.io>
1 parent 73ba549 commit 2244f8e

File tree

9 files changed

+179
-4
lines changed

9 files changed

+179
-4
lines changed
Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2026 The Stdlib Authors.
5+
*
6+
* Licensed under the Apache License, Version 2.0 (the "License");
7+
* you may not use this file except in compliance with the License.
8+
* You may obtain a copy of the License at
9+
*
10+
* http://www.apache.org/licenses/LICENSE-2.0
11+
*
12+
* Unless required by applicable law or agreed to in writing, software
13+
* distributed under the License is distributed on an "AS IS" BASIS,
14+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15+
* See the License for the specific language governing permissions and
16+
* limitations under the License.
17+
*
18+
*
19+
* ## Notice
20+
*
21+
* The original C code and copyright notice are from the [Cephes Mathematical Library]{@link https://www.netlib.org/cephes/}. The implementation has been modified for JavaScript.
22+
*
23+
* ```text
24+
* (C) Copyright Stephen L. Moshier 1984, 1987, 1992, 2000.
25+
*
26+
* Use, modification and distribution are subject to the
27+
* Cephes Mathematical Library License. (See accompanying file
28+
* LICENSE or copy at https://smath.com/en-US/view/CephesMathLibrary/license)
29+
* ```
30+
*/
31+
32+
'use strict';
33+
34+
// MODULES //
35+
36+
var abs = require( '@stdlib/math/base/special/abs' );
37+
var max = require( '@stdlib/math/base/special/max' );
38+
39+
40+
// MAIN //
41+
42+
/**
43+
* Evaluates 2F1(a, b; b; x) when `b = c` is a negative integer using AMS55 #15.4.2.
44+
*
45+
* @private
46+
* @param {number} a - first parameter
47+
* @param {number} b - second parameter (equals c, a non-positive integer)
48+
* @param {number} x - argument
49+
* @returns {number} function value, or NaN if precision is insufficient
50+
*/
51+
function hyp2f1NegCEqualBC( a, b, x ) {
52+
var collectorMax;
53+
var collector;
54+
var sum;
55+
var k;
56+
57+
collectorMax = 1.0;
58+
collector = 1.0;
59+
sum = 1.0;
60+
61+
if ( abs( b ) >= 1.0e5 ) {
62+
return NaN;
63+
}
64+
for ( k = 1; k <= -b; k++ ) {
65+
collector *= ( ( a + k - 1.0 ) * x ) / k;
66+
collectorMax = max( abs( collector ), collectorMax );
67+
sum += collector;
68+
}
69+
if ( 1.0e-16 * ( 1.0 + ( collectorMax / abs( sum ) ) ) > 1.0e-7 ) {
70+
return NaN;
71+
}
72+
return sum;
73+
}
74+
75+
76+
// EXPORTS //
77+
78+
module.exports = hyp2f1NegCEqualBC;

lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/hyt2f1.js

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ function hyt2f1( a, b, c, x, loss ) {
122122
d = c - a - b;
123123
id = round( d );
124124

125-
if ( x > 0.9 && !negIntA && !negIntB ) {
125+
if ( x > 0.85 && !negIntA && !negIntB ) {
126126
if ( isInteger( d ) === false ) {
127127
// Try the power series first:
128128
y = hys2f1( a, b, c, x, err );
@@ -223,7 +223,11 @@ function hyt2f1( a, b, c, x, loss ) {
223223
y = -y;
224224
}
225225
q = pow( s, id );
226-
y = ( id > 0.0 ) ? y*q : y1*q;
226+
if ( id > 0.0 ) {
227+
y *= q;
228+
} else {
229+
y1 *= q;
230+
}
227231
y += y1;
228232
}
229233
return {

lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/main.js

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ var gamma = require( '@stdlib/math/base/special/gamma' );
4040
var pow = require( '@stdlib/math/base/special/pow' );
4141
var abs = require( '@stdlib/math/base/special/abs' );
4242
var isNonPositiveInteger = require( './isnonpositiveinteger.js' );
43+
var hyp2f1NegCEqualBC = require( './hyp2f1negcequalbc.js' );
4344
var isInteger = require( './isinteger.js' );
4445
var hys2f1 = require( './hys2f1.js' );
4546
var hyt2f1 = require( './hyt2f1.js' );
@@ -156,6 +157,10 @@ function hyp2f1( a, b, c, x ) {
156157
if ( ax < 1.0 || x === -1.0 ) {
157158
if ( b === c ) {
158159
// 2F1(a,b;b;x) = (1-x)**(-a):
160+
if ( negIntB ) {
161+
// For negative integer b=c use the finite polynomial (AMS55 #15.4.2):
162+
return hyp2f1NegCEqualBC( a, b, x );
163+
}
159164
return pow( s, -a );
160165
}
161166
if ( a === c ) {

lib/node_modules/@stdlib/math/base/special/hyp2f1/test/fixtures/python/edge_cases3.json

Lines changed: 1 addition & 1 deletion
Large diffs are not rendered by default.
Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
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]}
Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
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]}
Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +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]}
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.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]}

lib/node_modules/@stdlib/math/base/special/hyp2f1/test/fixtures/python/runner.py

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -256,6 +256,46 @@ def main():
256256
# Outliers:
257257
gen_outliers("outliers.json")
258258

259+
# Edge Cases #5
260+
# d = c - a - b is a negative integer, x > 0.85, a and b are positive integers.
261+
pts = [
262+
(1.5, 2.5, 3.0, 0.90),
263+
(1.5, 2.5, 3.0, 0.95),
264+
(2.5, 1.5, 3.0, 0.90),
265+
(0.5, 2.5, 2.0, 0.92),
266+
(1.5, 2.5, 2.0, 0.91),
267+
(0.5, 1.5, 1.0, 0.90),
268+
(2.0, 3.0, 4.0, 0.90),
269+
(1.5, 3.5, 4.0, 0.88),
270+
(2.5, 2.5, 4.0, 0.93),
271+
(1.5, 3.5, 3.0, 0.91),
272+
]
273+
a, b, c, x = zip(*pts)
274+
a = np.array(a)
275+
b = np.array(b)
276+
c = np.array(c)
277+
x = np.array(x)
278+
gen(a, b, c, x, "edge_cases5.json")
279+
280+
# Edge Cases #6
281+
# b = c is a negative integer: 2F1(a, b; b; x) -> finite polynomial sum.
282+
pts = [
283+
(2.0, -1.0, -1.0, 0.5),
284+
(2.0, -2.0, -2.0, 0.5),
285+
(3.0, -2.0, -2.0, 0.3),
286+
(0.5, -3.0, -3.0, 0.7),
287+
(1.5, -3.0, -3.0, -0.5),
288+
(2.0, -4.0, -4.0, 0.6),
289+
(4.0, -2.0, -2.0, 0.8),
290+
(0.5, -1.0, -1.0, -0.9),
291+
]
292+
a, b, c, x = zip(*pts)
293+
a = np.array(a)
294+
b = np.array(b)
295+
c = np.array(c)
296+
x = np.array(x)
297+
gen(a, b, c, x, "edge_cases6.json")
298+
259299

260300
if __name__ == "__main__":
261301
main()

lib/node_modules/@stdlib/math/base/special/hyp2f1/test/test.js

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,8 @@ var edgeCases1 = require( './fixtures/python/edge_cases1.json' );
3434
var edgeCases2 = require( './fixtures/python/edge_cases2.json' );
3535
var edgeCases3 = require( './fixtures/python/edge_cases3.json' );
3636
var edgeCases4 = require( './fixtures/python/edge_cases4.json' );
37+
var edgeCases5 = require( './fixtures/python/edge_cases5.json' );
38+
var edgeCases6 = require( './fixtures/python/edge_cases6.json' );
3739
var outliers = require( './fixtures/python/outliers.json' );
3840

3941

@@ -238,6 +240,50 @@ tape( 'the function correctly evaluates the hypergeometric function', function t
238240
t.end();
239241
});
240242

243+
tape( 'the function correctly evaluates the hypergeometric function', function test( t ) {
244+
var expected;
245+
var a;
246+
var b;
247+
var c;
248+
var x;
249+
var v;
250+
var i;
251+
252+
a = edgeCases5.a;
253+
b = edgeCases5.b;
254+
c = edgeCases5.c;
255+
x = edgeCases5.x;
256+
expected = edgeCases5.expected;
257+
258+
for ( i = 0; i < x.length; i++ ) {
259+
v = hyp2f1( a[ i ], b[ i ], c[ i ], x[ i ] );
260+
t.strictEqual( isAlmostSameValue( v, expected[ i ], 12 ), true, 'returns expected value.' );
261+
}
262+
t.end();
263+
});
264+
265+
tape( 'the function correctly evaluates the hypergeometric function', function test( t ) {
266+
var expected;
267+
var a;
268+
var b;
269+
var c;
270+
var x;
271+
var v;
272+
var i;
273+
274+
a = edgeCases6.a;
275+
b = edgeCases6.b;
276+
c = edgeCases6.c;
277+
x = edgeCases6.x;
278+
expected = edgeCases6.expected;
279+
280+
for ( i = 0; i < x.length; i++ ) {
281+
v = hyp2f1( a[ i ], b[ i ], c[ i ], x[ i ] );
282+
t.strictEqual( isAlmostSameValue( v, expected[ i ], 1 ), true, 'returns expected value.' );
283+
}
284+
t.end();
285+
});
286+
241287
tape( 'the function correctly evaluates the hypergeometric function', function test( t ) {
242288
var expected;
243289
var a;

0 commit comments

Comments
 (0)