Skip to content

Commit 438f8e0

Browse files
anandkaranubcanandkaranubc
committed
feat: temp commit
1 parent fd45831 commit 438f8e0

File tree

4 files changed

+598
-0
lines changed

4 files changed

+598
-0
lines changed

‎lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/hyp2f1ra.js‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2018 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+
'use strict';
20+
21+
// MODULES //
22+
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' );
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);
34+
35+
// VARIABLES //
36+
37+
var EPS = 1.0e-13;
38+
var ETHRESH = 1.0e-12;
39+
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;
64+
65+
loss = 0.0; // Using an object to store loss
66+
err = 0.0;
67+
68+
if (abs(da) > MAX_ITERATIONS) {
69+
loss = 1.0;
70+
return { value: NaN, error: loss};
71+
}
72+
73+
if (da < 0) {
74+
// Recurse down
75+
f2Val = 0;
76+
f1 = hys2f1(t, b, c, x, err);
77+
loss += f1.error;
78+
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+
f1Val = f1.value;
83+
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;
95+
err = f1.error;
96+
f0 = hys2f1(t + 1, b, c, x, err); // CHECK IF err is THE SAME OR NEW ONE
97+
loss += f0.error;
98+
t += 1;
99+
f1Val = f1.value;
100+
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 };
110+
}
111+
112+
module.exports = hyp2f1ra;

‎lib/node_modules/@stdlib/math/base/special/hyp2f1/lib/hys2f1.js‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2018 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+
'use strict';
20+
21+
// MODULES //
22+
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' );
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+
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);
36+
37+
// VARIABLES //
38+
39+
var EPS = 1.0e-13;
40+
var ETHRESH = 1.0e-12;
41+
var MAX_ITERATIONS = 10000;
42+
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;
47+
48+
if (abs(b) > abs(a)) {
49+
f = b;
50+
b = a;
51+
a = f;
52+
}
53+
54+
if (isNonPositiveInteger(b) && abs(b) < abs(a)) {
55+
// .. except when `b` is a smaller negative integer
56+
f = b;
57+
b = a;
58+
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) {
87+
umax = k;
88+
}
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);
98+
99+
return { value: s, error: loss };
100+
}
101+
102+
module.exports = hys2f1;

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

Lines changed: 180 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2018 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+
'use strict';
20+
21+
// MODULES //
22+
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' );
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);
34+
35+
// VARIABLES //
36+
37+
var EPS = 1.0e-13;
38+
var ETHRESH = 1.0e-12;
39+
var MAX_ITERATIONS = 10000;
40+
41+
function hyt2f1(a, b, c, x, loss) {
42+
let p, q, r, s, t, y, w, d, err, err1;
43+
let ax, id, d1, d2, e, y1;
44+
let i, aid, sign;
45+
var val;
46+
47+
let neg_int_a = isNonPositiveInteger(a);
48+
let neg_int_b = isNonPositiveInteger(b);
49+
50+
err = 0.0;
51+
err1 = 0.0;
52+
s = 1.0 - x;
53+
54+
if (x < -0.5 && !(neg_int_a || neg_int_b)) {
55+
if (b > a) {
56+
y = hys2f1(a, c - b, c, -x / s, err);
57+
val = pow(s, -a) * y.value;
58+
}
59+
else {
60+
y = hys2f1(a, c - b, c, -x / s, err);
61+
val = pow(s, -b) * y.value;
62+
}
63+
loss = y.error;
64+
return { value: val, error: loss };
65+
}
66+
67+
d = c - a - b;
68+
id = round(d)
69+
70+
if (x > 0.9 && !(neg_int_a || neg_int_b)) {
71+
if (!isInteger(d)) {
72+
let sgngam;
73+
74+
y = hys2f1(a, b, c, x, err);
75+
if (y.error < ETHRESH) {
76+
return y;
77+
}
78+
79+
err = y.error;
80+
q = hys2f1(a, b, 1.0 - d, s, err);
81+
qVal = q.value;
82+
err = q.error;
83+
sign = 1;
84+
w = gammaln(d);
85+
sign *= gammasgn(d);
86+
w -= gammaln(c - a);
87+
sign *= gammasgn(c-a);
88+
w -= gammaln(c - b);
89+
sign *= gammasgn(c-b);
90+
qVal *= sign * exp(w);
91+
92+
r = hys2f1(c - a, c - b, d + 1.0, s, err1)
93+
err1 = r.error;
94+
rVal = pow(s, d) * r.value;
95+
sign = 1;
96+
w = gammaln(-d);
97+
sign *= gammasgn(-d);
98+
w -= gammaln(a);
99+
sign *= gammasgn(a);
100+
w -= gammaln(b);
101+
sign *= gammasgn(b);
102+
rVal *= sign * exp(w);
103+
104+
y = qVal + rVal;
105+
err += err1 + (EPS * max(abs(qVal), abs(rVal))) / y;
106+
y *= gamma(c);
107+
return { value: y, error: err };
108+
}
109+
else {
110+
if(id >= 0.0) {
111+
e = d;
112+
d1 = d;
113+
d2 = 0.0;
114+
aid = id;
115+
}
116+
else {
117+
e = -d;
118+
d1 = 0.0;
119+
d2 = d;
120+
aid = -id;
121+
}
122+
123+
ax = log(s);
124+
y = digamma(1.0) + digamma(1.0 + e) - digamma(a + d1) - digamma(b + d1) - ax;
125+
y /= gamma(e + 1.0);
126+
127+
p = (a + d1) * (b + d1) * s / gamma(e + 2.0);
128+
t = 1.0;
129+
130+
do {
131+
r = digamma(1.0 + t) + digamma(1.0 + t + e) - digamma(a + t + d1) - digamma(b + t + d1) - ax;
132+
q = p * r;
133+
y += q;
134+
p *= s * (a + t + d1) / (t + 1.0);
135+
p *= (b + t + d1) / (t + 1.0 + e);
136+
t += 1.0;
137+
138+
if (t > MAX_ITERATIONS) {
139+
loss = 1.0;
140+
return { value: NaN, error: loss };
141+
}
142+
} while (y === 0 || abs(q / y) > EPS);
143+
144+
if (id === 0.0) {
145+
y *= gamma(c) / (gamma(a) * gamma(b));
146+
return { value: y, loss: err };
147+
}
148+
149+
y1 = 1.0;
150+
151+
if (aid !== 1) {
152+
t = 0.0;
153+
p = 1.0;
154+
for (i = 1; i < aid; i++) {
155+
r = 1.0 - e + t;
156+
p *= s * (a + t + d2) * (b + t + d2) / r;
157+
t += 1.0;
158+
p /= t;
159+
y1 += p;
160+
}
161+
}
162+
163+
p = gamma(c);
164+
y1 *= gamma(e) * p / (gamma(a + d1) * gamma(b + d1));
165+
y *= p / (gamma(a + d2) * gamma(b + d2));
166+
if (aid % 2 !== 0) {
167+
y = -y;
168+
}
169+
q = pow(s, id);
170+
y = (id > 0.0) ? y * q : y1 * q;
171+
y += y1;
172+
return { value: y, loss: err };
173+
}
174+
}
175+
176+
y = hys2f1(a, b, c, x, err);
177+
return y;
178+
}
179+
180+
module.exports = hyt2f1;

0 commit comments

Comments
 (0)