22
22
23
23
var isNonPositiveInteger = require ( '@stdlib/math/base/assert/is-nonpositive-integer' ) ;
24
24
var PINF = require ( '@stdlib/constants/float64/pinf' ) ;
25
+ var round = require ( '@stdlib/math/base/special/round' ) ;
25
26
var abs = require ( '@stdlib/math/base/special/abs' ) ;
26
- var hyp2f1ra = require ( './hyp2f1ra.js' ) ;
27
27
28
28
29
29
// VARIABLES //
@@ -32,6 +32,96 @@ var EPS = 1.0e-13;
32
32
var MAX_ITERATIONS = 10000 ;
33
33
34
34
35
+ // FUNCTIONS //
36
+
37
+ /**
38
+ * Evaluates the Gaussian hypergeometric function by two-term recurrence in `a`.
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 hyp2f1ra ( a , b , c , x , loss ) {
50
+ var f2Val ;
51
+ var f1Val ;
52
+ var f0Val ;
53
+ var err ;
54
+ var da ;
55
+ var f1 ;
56
+ var f0 ;
57
+ var t ;
58
+ var n ;
59
+
60
+ loss = 0.0 ;
61
+ err = 0.0 ;
62
+
63
+ if ( ( c < 0 && a <= c ) || ( c >= 0 && a >= c ) ) {
64
+ da = round ( a - c ) ;
65
+ }
66
+ else {
67
+ da = round ( a ) ;
68
+ }
69
+
70
+ t = a - da ;
71
+ if ( abs ( da ) > MAX_ITERATIONS ) {
72
+ loss = 1.0 ;
73
+ return {
74
+ 'value' : NaN ,
75
+ 'error' : loss
76
+ } ;
77
+ }
78
+
79
+ if ( da < 0 ) {
80
+ f2Val = 0 ;
81
+ f1 = hys2f1 ( t , b , c , x , err ) ;
82
+ loss += f1 . error ;
83
+ err = f1 . error ;
84
+ f0 = hys2f1 ( t - 1 , b , c , x , err ) ;
85
+ loss += f0 . error ;
86
+ t -= 1 ;
87
+ f1Val = f1 . value ;
88
+ f0Val = f0 . value ;
89
+ for ( n = 1 ; n < - da ; ++ n ) {
90
+ f2Val = f1Val ;
91
+ f1Val = f0Val ;
92
+
93
+ // eslint-disable-next-line max-len
94
+ f0Val = - ( ( ( ( ( 2 * t ) - c ) - ( t * x ) + ( b * x ) ) * f1Val ) + ( ( t * ( x - 1 ) ) * f2Val ) ) / ( c - t ) ;
95
+ t -= 1 ;
96
+ }
97
+ }
98
+ else {
99
+ f2Val = 0 ;
100
+ f1 = hys2f1 ( t , b , c , x , err ) ;
101
+ loss += f1 . error ;
102
+ err = f1 . error ;
103
+ f0 = hys2f1 ( t + 1 , b , c , x , err ) ; // CHECK IF err is THE SAME OR NEW ONE
104
+ loss += f0 . error ;
105
+ t += 1 ;
106
+ f1Val = f1 . value ;
107
+ f0Val = f0 . value ;
108
+ for ( n = 1 ; n < da ; ++ n ) {
109
+ f2Val = f1Val ;
110
+ f1Val = f0Val ;
111
+
112
+ // eslint-disable-next-line max-len
113
+ f0Val = - ( ( ( ( ( 2 * t ) - c ) - ( t * x ) + ( b * x ) ) * f1Val ) + ( ( c - t ) * f2Val ) ) / ( t * ( x - 1 ) ) ;
114
+ t += 1 ;
115
+ }
116
+ }
117
+
118
+ return {
119
+ 'value' : f0Val ,
120
+ 'error' : loss
121
+ } ;
122
+ }
123
+
124
+
35
125
// MAIN //
36
126
37
127
/**
@@ -90,7 +180,10 @@ function hys2f1( a, b, c, x, loss ) {
90
180
do {
91
181
if ( abs ( h ) < EPS ) {
92
182
loss = 1.0 ;
93
- return PINF ;
183
+ return {
184
+ 'value' : PINF ,
185
+ 'error' : loss
186
+ } ;
94
187
}
95
188
m = k + 1.0 ;
96
189
u *= ( ( f + k ) * ( g + k ) * x / ( ( h + k ) * m ) ) ;
@@ -103,7 +196,10 @@ function hys2f1( a, b, c, x, loss ) {
103
196
i += 1 ;
104
197
if ( i > MAX_ITERATIONS ) {
105
198
loss = 1.0 ;
106
- return s ;
199
+ return {
200
+ 'value' : s ,
201
+ 'error' : loss
202
+ } ;
107
203
}
108
204
} while ( s === 0 || abs ( u / s ) > EPS ) ;
109
205
0 commit comments