Skip to content

Commit 2ebec87

Browse files
committed
feat: initial files
--- 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 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: missing_dependencies - task: lint_r status: na - task: lint_c_src status: na - 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: na - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- type: pre_push_report description: Results of running various checks prior to pushing changes. report: - task: run_javascript_examples status: na - task: run_c_examples status: na - task: run_cpp_examples status: na - task: run_javascript_readme_examples status: na - task: run_c_benchmarks status: na - task: run_cpp_benchmarks status: na - task: run_fortran_benchmarks status: na - task: run_javascript_benchmarks status: na - task: run_julia_benchmarks status: na - task: run_python_benchmarks status: na - task: run_r_benchmarks status: na - task: run_javascript_tests status: na ---
1 parent 438f8e0 commit 2ebec87

File tree

11 files changed

+714
-338
lines changed

11 files changed

+714
-338
lines changed
Lines changed: 85 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/**
22
* @license Apache-2.0
33
*
4-
* Copyright (c) 2018 The Stdlib Authors.
4+
* Copyright (c) 2025 The Stdlib Authors.
55
*
66
* Licensed under the Apache License, Version 2.0 (the "License");
77
* you may not use this file except in compliance with the License.
@@ -20,93 +20,106 @@
2020

2121
// MODULES //
2222

23-
var abs = require( '@stdlib/math/base/special/abs' );
2423
var round = require( '@stdlib/math/base/special/round' );
25-
var isInteger = require( '@stdlib/math/base/assert/is-integer' );
26-
var gamma = require( '@stdlib/math/base/special/gamma' );
27-
var isNegativeInteger = require( '@stdlib/math/base/assert/is-negative-integer' );
28-
var isNonPositiveInteger = require( '@stdlib/math/base/assert/is-nonpositive-integer' );
29-
var PINF = require( '@stdlib/constants/float64/pinf' );
30-
var UINT32_MAX = require( '@stdlib/constants/uint32/max' );
31-
// hyt2f1(double a, double b, double c, double x, double *loss);
32-
// hys2f1(double a, double b, double c, double x, double *loss);
33-
// hyp2f1ra(double a, double b, double c, double x, double *loss);
24+
var hys2f1 = require( './hys2f1.js' );
25+
3426

3527
// VARIABLES //
3628

37-
var EPS = 1.0e-13;
38-
var ETHRESH = 1.0e-12;
3929
var MAX_ITERATIONS = 10000;
40-
// const MAX_ITERATIONS = 1000; // Define a reasonable iteration limit
41-
42-
/*
43-
* Evaluate hypergeometric function by two-term recurrence in `a`.
44-
*
45-
* This avoids some of the loss of precision in the strongly alternating
46-
* hypergeometric series, and can be used to reduce the `a` and `b` parameters
47-
* to smaller values.
48-
*
49-
* AMS55 #15.2.10
50-
*/
51-
function hyp2f1ra(a, b, c, x, loss) {
52-
let f2, f1, f0;
53-
let n;
54-
let t, err, da;
55-
56-
// Don't cross c or zero
57-
if ((c < 0 && a <= c) || (c >= 0 && a >= c)) {
58-
da = round(a - c);
59-
}
60-
else {
61-
da = round(a);
62-
}
63-
t = a - da;
6430

65-
loss = 0.0; // Using an object to store loss
31+
32+
// MAIN //
33+
34+
/**
35+
* Evaluates the Gaussian hypergeometric function by two-term recurrence in `a`.
36+
*
37+
* @private
38+
* @param {number} a - input value
39+
* @param {number} b - input value
40+
* @param {number} c - input value
41+
* @param {number} x - input value
42+
* @param {number} loss - starting loss of significance
43+
* @returns {Object} the function value and error
44+
*
45+
*/
46+
function hyp2f1ra( a, b, c, x, loss ) {
47+
var f2Val;
48+
var f1Val;
49+
var f0Val;
50+
var err;
51+
var abs;
52+
var da;
53+
var f1;
54+
var f0;
55+
var t;
56+
var n;
57+
58+
loss = 0.0;
6659
err = 0.0;
6760

68-
if (abs(da) > MAX_ITERATIONS) {
69-
loss = 1.0;
70-
return { value: NaN, error: loss};
71-
}
61+
if ( ( c < 0 && a <= c ) || ( c >= 0 && a >= c ) ) {
62+
da = round( a - c );
63+
}
64+
else {
65+
da = round( a );
66+
}
67+
68+
t = a - da;
69+
if ( abs( da ) > MAX_ITERATIONS ) {
70+
loss = 1.0;
71+
return {
72+
'value': NaN,
73+
'error': loss
74+
};
75+
}
7276

73-
if (da < 0) {
74-
// Recurse down
75-
f2Val = 0;
76-
f1 = hys2f1(t, b, c, x, err);
77-
loss += f1.error;
77+
if ( da < 0 ) {
78+
f2Val = 0;
79+
f1 = hys2f1( t, b, c, x, err );
80+
loss += f1.error;
7881
err = f1.error;
79-
f0 = hys2f1(t - 1, b, c, x, err); // CHECK IF err is THE SAME OR NEW ONE
80-
loss += f0.error;
81-
t -= 1;
82+
f0 = hys2f1( t - 1, b, c, x, err );
83+
loss += f0.error;
84+
t -= 1;
8285
f1Val = f1.value;
8386
f0Val = f0.value;
84-
for (n = 1; n < -da; ++n) {
85-
f2Val = f1Val;
86-
f1Val = f0Val;
87-
f0Val = -((2 * t - c - t * x + b * x) * f1Val - t * (x - 1) * f2Val) / (c - t);
88-
t -= 1;
89-
}
90-
} else {
91-
// Recurse up
92-
f2Val = 0;
93-
f1 = hys2f1(t, b, c, x, err);
94-
loss += f1.error;
87+
for ( n = 1; n < -da; ++n ) {
88+
f2Val = f1Val;
89+
f1Val = f0Val;
90+
91+
// eslint-disable-next-line max-len
92+
f0Val = -(((((2 * t) - c) - (t * x) + (b * x)) * f1Val) - ((t * (x - 1)) * f2Val)) / (c - t);
93+
t -= 1;
94+
}
95+
}
96+
else {
97+
f2Val = 0;
98+
f1 = hys2f1( t, b, c, x, err );
99+
loss += f1.error;
95100
err = f1.error;
96-
f0 = hys2f1(t + 1, b, c, x, err); // CHECK IF err is THE SAME OR NEW ONE
101+
f0 = hys2f1( t + 1, b, c, x, err ); // CHECK IF err is THE SAME OR NEW ONE
97102
loss += f0.error;
98-
t += 1;
103+
t += 1;
99104
f1Val = f1.value;
100105
f0Val = f0.value;
101-
for (n = 1; n < da; ++n) {
102-
f2Val = f1Val;
103-
f1Val = f0Val;
104-
f0Val = -((2 * t - c - t * x + b * x) * f1Val + (c - t) * f2Val) / (t * (x - 1));
105-
t += 1;
106-
}
107-
}
108-
109-
return { value: f0, error: loss };
106+
for ( n = 1; n < da; ++n ) {
107+
f2Val = f1Val;
108+
f1Val = f0Val;
109+
110+
// eslint-disable-next-line max-len
111+
f0Val = -(((((2 * t) - c) - (t * x) + (b * x)) * f1Val) + ((c - t) * f2Val)) / (t * (x - 1));
112+
t += 1;
113+
}
114+
}
115+
116+
return {
117+
'value': f0,
118+
'error': loss
119+
};
110120
}
111121

122+
123+
// EXPORTS //
124+
112125
module.exports = hyp2f1ra;
Lines changed: 79 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
/**
22
* @license Apache-2.0
33
*
4-
* Copyright (c) 2018 The Stdlib Authors.
4+
* Copyright (c) 2025 The Stdlib Authors.
55
*
66
* Licensed under the Apache License, Version 2.0 (the "License");
77
* you may not use this file except in compliance with the License.
@@ -20,83 +20,101 @@
2020

2121
// MODULES //
2222

23-
var abs = require( '@stdlib/math/base/special/abs' );
24-
var round = require( '@stdlib/math/base/special/round' );
25-
var isInteger = require( '@stdlib/math/base/assert/is-integer' );
26-
var gamma = require( '@stdlib/math/base/special/gamma' );
27-
var isNegativeInteger = require( '@stdlib/math/base/assert/is-negative-integer' );
2823
var isNonPositiveInteger = require( '@stdlib/math/base/assert/is-nonpositive-integer' );
2924
var PINF = require( '@stdlib/constants/float64/pinf' );
30-
var UINT32_MAX = require( '@stdlib/constants/uint32/max' );
31-
var EPS = require( '@stdlib/constants/float64/eps' );
32-
const { error } = require('console');
33-
// hyt2f1(double a, double b, double c, double x, double *loss);
34-
// hys2f1(double a, double b, double c, double x, double *loss);
35-
// hyp2f1ra(double a, double b, double c, double x, double *loss);
25+
var abs = require( '@stdlib/math/base/special/abs' );
26+
var hyp2f1ra = require( './hyp2f1ra.js' );
27+
3628

3729
// VARIABLES //
3830

3931
var EPS = 1.0e-13;
40-
var ETHRESH = 1.0e-12;
4132
var MAX_ITERATIONS = 10000;
4233

43-
function hys2f1(a, b, c, x, loss) {
44-
let f, g, h, k, m, s, u, umax;
45-
let i;
46-
let ib, intflag = 0;
4734

48-
if (abs(b) > abs(a)) {
49-
f = b;
35+
// MAIN //
36+
37+
/**
38+
* Evaluates the power series expansion of Gaussian hypergeometric function.
39+
*
40+
* @private
41+
* @param {number} a - input value
42+
* @param {number} b - input value
43+
* @param {number} c - input value
44+
* @param {number} x - input value
45+
* @param {number} loss - starting loss of significance
46+
* @returns {Object} the function value and error
47+
*
48+
*/
49+
function hys2f1( a, b, c, x, loss ) {
50+
var intFlag;
51+
var umax;
52+
var f;
53+
var g;
54+
var h;
55+
var k;
56+
var m;
57+
var s;
58+
var u;
59+
var i;
60+
61+
intFlag = 0;
62+
63+
if ( abs( b ) > abs( a ) ) {
64+
f = b;
5065
b = a;
5166
a = f;
52-
}
67+
}
5368

54-
if (isNonPositiveInteger(b) && abs(b) < abs(a)) {
55-
// .. except when `b` is a smaller negative integer
56-
f = b;
69+
if ( isNonPositiveInteger( b ) && abs( b ) < abs( a ) ) {
70+
f = b;
5771
b = a;
5872
a = f;
59-
intflag = 1;
60-
}
61-
62-
if ((abs(a) > abs(c) + 1 || intflag) && abs(c - a) > 2 && abs(a) > 2) {
63-
// |a| >> |c| implies that large cancellation error is to be expected.
64-
// We try to reduce it with the recurrence relations
65-
return hyp2f1ra(a, b, c, x, loss);
66-
}
67-
68-
i = 0;
69-
umax = 0.0;
70-
f = a;
71-
g = b;
72-
h = c;
73-
s = 1.0;
74-
u = 1.0;
75-
k = 0.0;
76-
77-
do {
78-
if (abs(h) < EPS) {
79-
loss = 1.0;
80-
return PINF;
81-
}
82-
m = k + 1.0;
83-
u = u * ((f + k) * (g + k) * x / ((h + k) * m));
84-
s += u;
85-
k = abs(u); // remember largest term summed
86-
if (k > umax) {
73+
intFlag = 1;
74+
}
75+
76+
// eslint-disable-next-line max-len
77+
if ( ( ( abs( a ) > abs( c ) + 1 ) || intFlag ) && ( abs( c - a ) > 2 ) && ( abs( a ) > 2 ) ) {
78+
return hyp2f1ra( a, b, c, x, loss );
79+
}
80+
81+
i = 0;
82+
umax = 0.0;
83+
f = a;
84+
g = b;
85+
h = c;
86+
s = 1.0;
87+
u = 1.0;
88+
k = 0.0;
89+
90+
do {
91+
if ( abs( h ) < EPS ) {
92+
loss = 1.0;
93+
return PINF;
94+
}
95+
m = k + 1.0;
96+
u *= ( ( f + k ) * ( g + k ) * x / ( ( h + k ) * m ) );
97+
s += u;
98+
k = abs( u );
99+
if ( k > umax ) {
87100
umax = k;
88101
}
89-
k = m;
90-
if (++i > MAX_ITERATIONS) {
91-
loss = 1.0;
92-
return s;
93-
}
94-
} while (s === 0 || abs(u / s) > EPS);
95-
96-
// return estimated relative error
97-
loss = (EPS * umax) / abs(s) + (EPS * i);
102+
k = m;
103+
i += 1;
104+
if ( i > MAX_ITERATIONS ) {
105+
loss = 1.0;
106+
return s;
107+
}
108+
} while ( s === 0 || abs( u / s ) > EPS );
98109

99-
return { value: s, error: loss };
110+
loss = ( ( EPS * umax ) / abs( s ) ) + ( EPS * i );
111+
return {
112+
'value': s,
113+
'error': loss
114+
};
100115
}
101116

117+
118+
// EXPORTS //
119+
102120
module.exports = hys2f1;

0 commit comments

Comments
 (0)