Skip to content

Commit 539d6f8

Browse files
devmotiondevmotion
authored
Use @evalpoly instead of Base.Math.@horner (#517)
1 parent 198fd80 commit 539d6f8

File tree

6 files changed

+150
-161
lines changed

6 files changed

+150
-161
lines changed

‎src/beta_inc.jl‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,3 @@
1-
using Base.Math: @horner
21
using Base.MPFR: ROUNDING_MODE
32
const exparg_n = log(nextfloat(floatmin(Float64)))
43
const exparg_p = log(prevfloat(floatmax(Float64)))
@@ -36,7 +35,7 @@ function _loggammadiv(a::Float64, b::Float64)
3635

3736
# SET W = stirling(b) - stirling(a+b)
3837
t = inv(b)^2
39-
w = @horner(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁)
38+
w = @evalpoly(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁)
4039
w *= c/b
4140

4241
#COMBINING
@@ -66,11 +65,11 @@ function stirling_corr(a0::Float64, b0::Float64)
6665
s₉ = 1.0 + (x + x²*s₇)
6766
s₁₁ = 1.0 + (x + x²*s₉)
6867
t = inv(b)^2
69-
w = @horner(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁)
68+
w = @evalpoly(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁)
7069
w *= c/b
7170
# COMPUTE stirling(a) + w
7271
t = inv(a)^2
73-
return @horner(t, .833333333333333E-01, -.277777777760991E-02, .793650666825390E-03, -.595202931351870E-03, .837308034031215E-03, -.165322962780713E-02)/a + w
72+
return @evalpoly(t, .833333333333333E-01, -.277777777760991E-02, .793650666825390E-03, -.595202931351870E-03, .837308034031215E-03, -.165322962780713E-02)/a + w
7473
end
7574

7675
@doc raw"""
@@ -948,7 +947,7 @@ function _beta_inc_inv(a::Float64, b::Float64, p::Float64, q::Float64=1-p)
948947
#Initial approx
949948
x = p
950949
r = sqrt(-2*log(p))
951-
p_approx = r - @horner(r, 2.30753e+00, 0.27061e+00) / @horner(r, 1.0, .99229e+00, .04481e+00)
950+
p_approx = r - @evalpoly(r, 2.30753e+00, 0.27061e+00) / @evalpoly(r, 1.0, .99229e+00, .04481e+00)
952951
fpu = floatmin(Float64)
953952
lb = logbeta(a, b)
954953

‎src/ellip.jl‎

Lines changed: 24 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,3 @@
1-
using Base.Math: @horner
2-
3-
4-
5-
61
@doc raw"""
72
ellipk(m)
83
@@ -57,31 +52,31 @@ function _ellipk(m::Float64)
5752

5853
elseif 0.0 < x < 0.1 #Table 2 from paper
5954
t = x-0.05
60-
t = @horner(t,
55+
t = @evalpoly(t,
6156
1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 ,
6257
0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 ,
6358
0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 ,
6459
0.085842591595413900 , 0.081541118718303215)
6560

6661
elseif 0.1 <= x < 0.2 #Table 3
6762
t = x-0.15
68-
t = @horner(t ,
63+
t = @evalpoly(t ,
6964
1.635256732264579992 , 0.471190626148732291 , 0.309728410831499587 ,
7065
0.252208311773135699 , 0.226725623219684650 , 0.215774446729585976 ,
7166
0.213108771877348910 , 0.216029124605188282 , 0.223255831633057896 ,
7267
0.234180501294209925 , 0.248557682972264071 , 0.266363809892617521)
7368

7469
elseif 0.2 <= x < 0.3 #Table 4
7570
t = x-0.25
76-
t = @horner(t ,
71+
t = @evalpoly(t ,
7772
1.685750354812596043 , 0.541731848613280329 , 0.401524438390690257 ,
7873
0.369642473420889090 , 0.376060715354583645 , 0.405235887085125919 ,
7974
0.453294381753999079 , 0.520518947651184205 , 0.609426039204995055 ,
8075
0.724263522282908870 , 0.871013847709812357 , 1.057652872753547036)
8176

8277
elseif 0.3 <= x < 0.4 #Table 5
8378
t = x-0.35
84-
t = @horner(t ,
79+
t = @evalpoly(t ,
8580
1.744350597225613243 , 0.634864275371935304 , 0.539842564164445538 ,
8681
0.571892705193787391 , 0.670295136265406100 , 0.832586590010977199 ,
8782
1.073857448247933265 , 1.422091460675497751 , 1.920387183402304829 ,
@@ -90,7 +85,7 @@ function _ellipk(m::Float64)
9085

9186
elseif 0.4 <= x < 0.5 #Table 6
9287
t = x-0.45
93-
t = @horner(t,
88+
t = @evalpoly(t,
9489
1.813883936816982644 , 0.763163245700557246 , 0.761928605321595831 ,
9590
0.951074653668427927 , 1.315180671703161215 , 1.928560693477410941 ,
9691
2.937509342531378755 , 4.594894405442878062 , 7.330071221881720772 ,
@@ -99,15 +94,15 @@ function _ellipk(m::Float64)
9994

10095
elseif 0.5 <= x < 0.6 #Table 7
10196
t = x-0.55
102-
t = @horner(t ,
97+
t = @evalpoly(t ,
10398
1.898924910271553526 , 0.950521794618244435 , 1.151077589959015808 ,
10499
1.750239106986300540 , 2.952676812636875180 , 5.285800396121450889 ,
105100
9.832485716659979747 , 18.78714868327559562 , 36.61468615273698145 ,
106101
72.45292395127771801 , 145.1079577347069102 , 293.4786396308497026 ,
107102
598.3851815055010179 , 1228.420013075863451 , 2536.529755382764488)
108103
elseif 0.6 <= x < 0.7 #Table 8
109104
t = x-0.65
110-
t = @horner(t ,
105+
t = @evalpoly(t ,
111106
2.007598398424376302 , 1.248457231212347337 , 1.926234657076479729 ,
112107
3.751289640087587680 , 8.119944554932045802 , 18.66572130873555361 ,
113108
44.60392484291437063 , 109.5092054309498377 , 274.2779548232413480 ,
@@ -117,7 +112,7 @@ function _ellipk(m::Float64)
117112

118113
elseif 0.7 <= x < 0.8 #Table 9
119114
t = x-0.75
120-
t = @horner(t,
115+
t = @evalpoly(t,
121116
2.156515647499643235 , 1.791805641849463243 , 3.826751287465713147 ,
122117
10.38672468363797208 , 31.40331405468070290 , 100.9237039498695416 ,
123118
337.3268282632272897 , 1158.707930567827917 , 4060.990742193632092 ,
@@ -128,7 +123,7 @@ function _ellipk(m::Float64)
128123

129124
elseif 0.8 <= x < 0.85 #Table 10
130125
t = x-0.825
131-
t = @horner(t ,
126+
t = @evalpoly(t ,
132127
2.318122621712510589 , 2.616920150291232841 , 7.897935075731355823 ,
133128
30.50239715446672327 , 131.4869365523528456 , 602.9847637356491617 ,
134129
2877.024617809972641 , 14110.51991915180325 , 70621.44088156540229 ,
@@ -138,7 +133,7 @@ function _ellipk(m::Float64)
138133

139134
elseif 0.85 <= x < 0.9 #Table 11
140135
t = x-0.875
141-
t = @horner(t,
136+
t = @evalpoly(t,
142137
2.473596173751343912 , 3.727624244118099310 , 15.60739303554930496 ,
143138
84.12850842805887747 , 506.9818197040613935 , 3252.277058145123644 ,
144139
21713.24241957434256 , 149037.0451890932766 , 1043999.331089990839 ,
@@ -151,14 +146,14 @@ function _ellipk(m::Float64)
151146
@assert 0.9 <= x < 1
152147
td = 1-x
153148
td1 = td-0.05
154-
qd = @horner(td,
149+
qd = @evalpoly(td,
155150
0.0, (1.0/16.0), (1.0/32.0), (21.0/1024.0), (31.0/2048.0), (6257.0/524288.0),
156151
(10293.0/1048576.0), (279025.0/33554432.0), (483127.0/67108864.0),
157152
(435506703.0/68719476736.0), (776957575.0/137438953472.0) ,
158153
(22417045555.0/4398046511104.0) , (40784671953.0/8796093022208.0) ,
159154
(9569130097211.0/2251799813685248.0) , (17652604545791.0/4503599627370496.0))
160155

161-
kdm = @horner(td1 ,
156+
kdm = @evalpoly(td1 ,
162157
1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 ,
163158
0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 ,
164159
0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 ,
@@ -229,39 +224,39 @@ function _ellipe(m::Float64)
229224

230225
elseif 0.0 < x < 0.1 #Table 2 from paper
231226
t = x-0.05
232-
t = @horner(t ,
227+
t = @evalpoly(t ,
233228
+1.550973351780472328 , -0.400301020103198524 , -0.078498619442941939 ,
234229
-0.034318853117591992 , -0.019718043317365499 , -0.013059507731993309 ,
235230
-0.009442372874146547 , -0.007246728512402157 , -0.005807424012956090 ,
236231
-0.004809187786009338)
237232

238233
elseif 0.1 <= x < 0.2 #Table 3
239234
t = x-0.15
240-
t = @horner(t ,
235+
t = @evalpoly(t ,
241236
+1.510121832092819728 , -0.417116333905867549 , -0.090123820404774569 ,
242237
-0.043729944019084312 , -0.027965493064761785 , -0.020644781177568105 ,
243238
-0.016650786739707238 , -0.014261960828842520 , -0.012759847429264803 ,
244239
-0.011799303775587354 , -0.011197445703074968)
245240

246241
elseif 0.2 <= x < 0.3 #Table 4
247242
t = x-0.25
248-
t = @horner(t ,
243+
t = @evalpoly(t ,
249244
+1.467462209339427155 , -0.436576290946337775 , -0.105155557666942554 ,
250245
-0.057371843593241730 , -0.041391627727340220 , -0.034527728505280841 ,
251246
-0.031495443512532783 , -0.030527000890325277 , -0.030916984019238900 ,
252247
-0.032371395314758122 , -0.034789960386404158)
253248

254249
elseif 0.3 <= x < 0.4 #Table 5
255250
t = x-0.35
256-
t = @horner(t,
251+
t = @evalpoly(t,
257252
+1.422691133490879171 , -0.459513519621048674 , -0.125250539822061878,
258253
-0.078138545094409477 , -0.064714278472050002 , -0.062084339131730311,
259254
-0.065197032815572477 , -0.072793895362578779 , -0.084959075171781003,
260255
-0.102539850131045997 , -0.127053585157696036 , -0.160791120691274606)
261256

262257
elseif 0.4 <= x < 0.5 #Table 6
263258
t = x-0.45
264-
t = @horner(t ,
259+
t = @evalpoly(t ,
265260
+1.375401971871116291 , -0.487202183273184837 , -0.153311701348540228 ,
266261
-0.111849444917027833 , -0.108840952523135768 , -0.122954223120269076 ,
267262
-0.152217163962035047 , -0.200495323642697339 , -0.276174333067751758 ,
@@ -270,15 +265,15 @@ function _ellipe(m::Float64)
270265

271266
elseif 0.5 <= x < 0.6 #Table 7
272267
t = x-0.55
273-
t = @horner(t ,
268+
t = @evalpoly(t ,
274269
+1.325024497958230082 , -0.521727647557566767 , -0.194906430482126213 ,
275270
-0.171623726822011264 , -0.202754652926419141 , -0.278798953118534762 ,
276271
-0.420698457281005762 , -0.675948400853106021 , -1.136343121839229244 ,
277272
-1.976721143954398261 , -3.531696773095722506 , -6.446753640156048150 ,
278273
-11.97703130208884026)
279274
elseif 0.6 <= x < 0.7 #Table 8
280275
t = x-0.65
281-
t = @horner(t,
276+
t = @evalpoly(t,
282277
+1.270707479650149744 , -0.566839168287866583 , -0.262160793432492598 ,
283278
-0.292244173533077419 , -0.440397840850423189 , -0.774947641381397458 ,
284279
-1.498870837987561088 , -3.089708310445186667 , -6.667595903381001064 ,
@@ -287,7 +282,7 @@ function _ellipe(m::Float64)
287282

288283
elseif 0.7 <= x < 0.8 #Table 9
289284
t = x-0.75
290-
t = @horner(t,
285+
t = @evalpoly(t,
291286
+1.211056027568459525 , -0.630306413287455807 , -0.387166409520669145 ,
292287
-0.592278235311934603 , -1.237555584513049844 , -3.032056661745247199 ,
293288
-8.181688221573590762 , -23.55507217389693250 , -71.04099935893064956 ,
@@ -297,7 +292,7 @@ function _ellipe(m::Float64)
297292

298293
elseif 0.8 <= x < 0.85 #Table 10
299294
t = x-0.825
300-
t = @horner(t,
295+
t = @evalpoly(t,
301296
+1.161307152196282836 , -0.701100284555289548 , -0.580551474465437362 ,
302297
-1.243693061077786614 , -3.679383613496634879 , -12.81590924337895775 ,
303298
-49.25672530759985272 , -202.1818735434090269 , -869.8602699308701437 ,
@@ -306,7 +301,7 @@ function _ellipe(m::Float64)
306301

307302
elseif 0.85 <= x < 0.9 #Table 11
308303
t = x-0.875
309-
t = @horner(t,
304+
t = @evalpoly(t,
310305
+1.124617325119752213 , -0.770845056360909542 , -0.844794053644911362 ,
311306
-2.490097309450394453 , -10.23971741154384360 , -49.74900546551479866 ,
312307
-267.0986675195705196 , -1532.665883825229947 , -9222.313478526091951 ,
@@ -319,12 +314,12 @@ function _ellipe(m::Float64)
319314
td = 1-x
320315
td1 = td-0.05
321316

322-
kdm = @horner(td1 ,
317+
kdm = @evalpoly(td1 ,
323318
1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 ,
324319
0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 ,
325320
0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 ,
326321
0.085842591595413900 , 0.081541118718303215)
327-
edm = @horner(td1 ,
322+
edm = @evalpoly(td1 ,
328323
+1.550973351780472328 , -0.400301020103198524 , -0.078498619442941939 ,
329324
-0.034318853117591992 , -0.019718043317365499 , -0.013059507731993309 ,
330325
-0.009442372874146547 , -0.007246728512402157 , -0.005807424012956090 ,

0 commit comments

Comments
 (0)