@@ -127,23 +127,23 @@ function _lgamma_r(x::Float64)
127
127
if f ≥ 0.7315998077392578
128
128
y = 1.0 + c - x
129
129
z = y * y
130
- p1 = evalpoly (z, ( 7.72156649015328655494e-02 , 6.73523010531292681824e-02 , 7.38555086081402883957e-03 , 1.19270763183362067845e-03 , 2.20862790713908385557e-04 , 2.52144565451257326939e-05 ) )
131
- p2 = z * evalpoly (z, ( 3.22467033424113591611e-01 , 2.05808084325167332806e-02 , 2.89051383673415629091e-03 , 5.10069792153511336608e-04 , 1.08011567247583939954e-04 , 4.48640949618915160150e-05 ) )
130
+ p1 = @ evalpoly (z, 7.72156649015328655494e-02 , 6.73523010531292681824e-02 , 7.38555086081402883957e-03 , 1.19270763183362067845e-03 , 2.20862790713908385557e-04 , 2.52144565451257326939e-05 )
131
+ p2 = z * @ evalpoly (z, 3.22467033424113591611e-01 , 2.05808084325167332806e-02 , 2.89051383673415629091e-03 , 5.10069792153511336608e-04 , 1.08011567247583939954e-04 , 4.48640949618915160150e-05 )
132
132
p = muladd (p1, y, p2)
133
133
r += muladd (y, - 0.5 , p)
134
134
elseif f ≥ 0.2316399812698364 # or, the lb? 0.2316322326660156
135
135
y = x - 0.46163214496836225 - c
136
136
z = y * y
137
137
w = z * y
138
- p1 = evalpoly (w, ( 4.83836122723810047042e-01 , - 3.27885410759859649565e-02 , 6.10053870246291332635e-03 , - 1.40346469989232843813e-03 , 3.15632070903625950361e-04 ) )
139
- p2 = evalpoly (w, ( - 1.47587722994593911752e-01 , 1.79706750811820387126e-02 , - 3.68452016781138256760e-03 , 8.81081882437654011382e-04 , - 3.12754168375120860518e-04 ) )
140
- p3 = evalpoly (w, ( 6.46249402391333854778e-02 , - 1.03142241298341437450e-02 , 2.25964780900612472250e-03 , - 5.38595305356740546715e-04 , 3.35529192635519073543e-04 ) )
138
+ p1 = @ evalpoly (w, 4.83836122723810047042e-01 , - 3.27885410759859649565e-02 , 6.10053870246291332635e-03 , - 1.40346469989232843813e-03 , 3.15632070903625950361e-04 )
139
+ p2 = @ evalpoly (w, - 1.47587722994593911752e-01 , 1.79706750811820387126e-02 , - 3.68452016781138256760e-03 , 8.81081882437654011382e-04 , - 3.12754168375120860518e-04 )
140
+ p3 = @ evalpoly (w, 6.46249402391333854778e-02 , - 1.03142241298341437450e-02 , 2.25964780900612472250e-03 , - 5.38595305356740546715e-04 , 3.35529192635519073543e-04 )
141
141
p = muladd (z, p1, - muladd (w, - muladd (p3, y, p2), - 3.63867699703950536541e-18 ))
142
142
r += p - 1.21486290535849611461e-1
143
143
else
144
144
y = x - c
145
- p1 = y * evalpoly (y, ( - 7.72156649015328655494e-02 , 6.32827064025093366517e-01 , 1.45492250137234768737 , 9.77717527963372745603e-01 , 2.28963728064692451092e-01 , 1.33810918536787660377e-02 ) )
146
- p2 = evalpoly (y, ( 1.0 , 2.45597793713041134822 , 2.12848976379893395361 , 7.69285150456672783825e-01 , 1.04222645593369134254e-01 , 3.21709242282423911810e-03 ) )
145
+ p1 = y * @ evalpoly (y, - 7.72156649015328655494e-02 , 6.32827064025093366517e-01 , 1.45492250137234768737 , 9.77717527963372745603e-01 , 2.28963728064692451092e-01 , 1.33810918536787660377e-02 )
146
+ p2 = @ evalpoly (y, 1.0 , 2.45597793713041134822 , 2.12848976379893395361 , 7.69285150456672783825e-01 , 1.04222645593369134254e-01 , 3.21709242282423911810e-03 )
147
147
r += muladd (y, - 0.5 , p1 / p2)
148
148
end
149
149
elseif ix < 0x40200000 #= x < 8.0 =#
@@ -157,13 +157,13 @@ function _lgamma_r(x::Float64)
157
157
u = x + p
158
158
z *= u
159
159
end
160
- p = y * evalpoly (y, ( - 7.72156649015328655494e-2 , 2.14982415960608852501e-1 , 3.25778796408930981787e-1 , 1.46350472652464452805e-1 , 2.66422703033638609560e-2 , 1.84028451407337715652e-3 , 3.19475326584100867617e-5 ) )
161
- q = evalpoly (y, ( 1.0 , 1.39200533467621045958 , 7.21935547567138069525e-1 , 1.71933865632803078993e-1 , 1.86459191715652901344e-2 , 7.77942496381893596434e-4 , 7.32668430744625636189e-6 ) )
160
+ p = y * @ evalpoly (y, - 7.72156649015328655494e-2 , 2.14982415960608852501e-1 , 3.25778796408930981787e-1 , 1.46350472652464452805e-1 , 2.66422703033638609560e-2 , 1.84028451407337715652e-3 , 3.19475326584100867617e-5 )
161
+ q = @ evalpoly (y, 1.0 , 1.39200533467621045958 , 7.21935547567138069525e-1 , 1.71933865632803078993e-1 , 1.86459191715652901344e-2 , 7.77942496381893596434e-4 , 7.32668430744625636189e-6 )
162
162
r = log (z) + muladd (0.5 , y, p / q)
163
163
elseif ix < 0x43900000 #= 8.0 ≤ x < 2^58 =#
164
164
z = 1.0 / x
165
165
y = z * z
166
- w = muladd (z, evalpoly (y, ( 8.33333333333329678849e-2 , - 2.77777777728775536470e-3 , 7.93650558643019558500e-4 , - 5.95187557450339963135e-4 , 8.36339918996282139126e-4 , - 1.63092934096575273989e-3 ) ), 4.18938533204672725052e-1 )
166
+ w = muladd (z, @ evalpoly (y, 8.33333333333329678849e-2 , - 2.77777777728775536470e-3 , 7.93650558643019558500e-4 , - 5.95187557450339963135e-4 , 8.36339918996282139126e-4 , - 1.63092934096575273989e-3 ), 4.18938533204672725052e-1 )
167
167
r = muladd (x - 0.5 , log (x) - 1.0 , w)
168
168
else #= 2^58 ≤ x ≤ Inf =#
169
169
r = muladd (x, log (x), - x)
0 commit comments