/
erf.go
339 lines (330 loc) · 11.4 KB
/
erf.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
// Copyright 2010 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package math
/*
Floating-point error function and complementary error function.
*/
// The original C code and the long comment below are
// from FreeBSD's /usr/src/lib/msun/src/s_erf.c and
// came with this notice. The go code is a simplified
// version of the original C.
//
// ====================================================
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
//
// Developed at SunPro, a Sun Microsystems, Inc. business.
// Permission to use, copy, modify, and distribute this
// software is freely granted, provided that this notice
// is preserved.
// ====================================================
//
//
// double erf(double x)
// double erfc(double x)
// x
// 2 |\
// erf(x) = --------- | exp(-t*t)dt
// sqrt(pi) \|
// 0
//
// erfc(x) = 1-erf(x)
// Note that
// erf(-x) = -erf(x)
// erfc(-x) = 2 - erfc(x)
//
// Method:
// 1. For |x| in [0, 0.84375]
// erf(x) = x + x*R(x**2)
// erfc(x) = 1 - erf(x) if x in [-.84375,0.25]
// = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375]
// where R = P/Q where P is an odd poly of degree 8 and
// Q is an odd poly of degree 10.
// -57.90
// | R - (erf(x)-x)/x | <= 2
//
//
// Remark. The formula is derived by noting
// erf(x) = (2/sqrt(pi))*(x - x**3/3 + x**5/10 - x**7/42 + ....)
// and that
// 2/sqrt(pi) = 1.128379167095512573896158903121545171688
// is close to one. The interval is chosen because the fix
// point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
// near 0.6174), and by some experiment, 0.84375 is chosen to
// guarantee the error is less than one ulp for erf.
//
// 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
// c = 0.84506291151 rounded to single (24 bits)
// erf(x) = sign(x) * (c + P1(s)/Q1(s))
// erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
// 1+(c+P1(s)/Q1(s)) if x < 0
// |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
// Remark: here we use the taylor series expansion at x=1.
// erf(1+s) = erf(1) + s*Poly(s)
// = 0.845.. + P1(s)/Q1(s)
// That is, we use rational approximation to approximate
// erf(1+s) - (c = (single)0.84506291151)
// Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
// where
// P1(s) = degree 6 poly in s
// Q1(s) = degree 6 poly in s
//
// 3. For x in [1.25,1/0.35(~2.857143)],
// erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
// erf(x) = 1 - erfc(x)
// where
// R1(z) = degree 7 poly in z, (z=1/x**2)
// S1(z) = degree 8 poly in z
//
// 4. For x in [1/0.35,28]
// erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
// = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
// = 2.0 - tiny (if x <= -6)
// erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
// erf(x) = sign(x)*(1.0 - tiny)
// where
// R2(z) = degree 6 poly in z, (z=1/x**2)
// S2(z) = degree 7 poly in z
//
// Note1:
// To compute exp(-x*x-0.5625+R/S), let s be a single
// precision number and s := x; then
// -x*x = -s*s + (s-x)*(s+x)
// exp(-x*x-0.5626+R/S) =
// exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
// Note2:
// Here 4 and 5 make use of the asymptotic series
// exp(-x*x)
// erfc(x) ~ ---------- * ( 1 + Poly(1/x**2) )
// x*sqrt(pi)
// We use rational approximation to approximate
// g(s)=f(1/x**2) = log(erfc(x)*x) - x*x + 0.5625
// Here is the error bound for R1/S1 and R2/S2
// |R1/S1 - f(x)| < 2**(-62.57)
// |R2/S2 - f(x)| < 2**(-61.52)
//
// 5. For inf > x >= 28
// erf(x) = sign(x) *(1 - tiny) (raise inexact)
// erfc(x) = tiny*tiny (raise underflow) if x > 0
// = 2 - tiny if x<0
//
// 7. Special case:
// erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
// erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
// erfc/erf(NaN) is NaN
const (
erx = 8.45062911510467529297e-01 // 0x3FEB0AC160000000
// Coefficients for approximation to erf in [0, 0.84375]
efx = 1.28379167095512586316e-01 // 0x3FC06EBA8214DB69
efx8 = 1.02703333676410069053e+00 // 0x3FF06EBA8214DB69
pp0 = 1.28379167095512558561e-01 // 0x3FC06EBA8214DB68
pp1 = -3.25042107247001499370e-01 // 0xBFD4CD7D691CB913
pp2 = -2.84817495755985104766e-02 // 0xBF9D2A51DBD7194F
pp3 = -5.77027029648944159157e-03 // 0xBF77A291236668E4
pp4 = -2.37630166566501626084e-05 // 0xBEF8EAD6120016AC
qq1 = 3.97917223959155352819e-01 // 0x3FD97779CDDADC09
qq2 = 6.50222499887672944485e-02 // 0x3FB0A54C5536CEBA
qq3 = 5.08130628187576562776e-03 // 0x3F74D022C4D36B0F
qq4 = 1.32494738004321644526e-04 // 0x3F215DC9221C1A10
qq5 = -3.96022827877536812320e-06 // 0xBED09C4342A26120
// Coefficients for approximation to erf in [0.84375, 1.25]
pa0 = -2.36211856075265944077e-03 // 0xBF6359B8BEF77538
pa1 = 4.14856118683748331666e-01 // 0x3FDA8D00AD92B34D
pa2 = -3.72207876035701323847e-01 // 0xBFD7D240FBB8C3F1
pa3 = 3.18346619901161753674e-01 // 0x3FD45FCA805120E4
pa4 = -1.10894694282396677476e-01 // 0xBFBC63983D3E28EC
pa5 = 3.54783043256182359371e-02 // 0x3FA22A36599795EB
pa6 = -2.16637559486879084300e-03 // 0xBF61BF380A96073F
qa1 = 1.06420880400844228286e-01 // 0x3FBB3E6618EEE323
qa2 = 5.40397917702171048937e-01 // 0x3FE14AF092EB6F33
qa3 = 7.18286544141962662868e-02 // 0x3FB2635CD99FE9A7
qa4 = 1.26171219808761642112e-01 // 0x3FC02660E763351F
qa5 = 1.36370839120290507362e-02 // 0x3F8BEDC26B51DD1C
qa6 = 1.19844998467991074170e-02 // 0x3F888B545735151D
// Coefficients for approximation to erfc in [1.25, 1/0.35]
ra0 = -9.86494403484714822705e-03 // 0xBF843412600D6435
ra1 = -6.93858572707181764372e-01 // 0xBFE63416E4BA7360
ra2 = -1.05586262253232909814e+01 // 0xC0251E0441B0E726
ra3 = -6.23753324503260060396e+01 // 0xC04F300AE4CBA38D
ra4 = -1.62396669462573470355e+02 // 0xC0644CB184282266
ra5 = -1.84605092906711035994e+02 // 0xC067135CEBCCABB2
ra6 = -8.12874355063065934246e+01 // 0xC054526557E4D2F2
ra7 = -9.81432934416914548592e+00 // 0xC023A0EFC69AC25C
sa1 = 1.96512716674392571292e+01 // 0x4033A6B9BD707687
sa2 = 1.37657754143519042600e+02 // 0x4061350C526AE721
sa3 = 4.34565877475229228821e+02 // 0x407B290DD58A1A71
sa4 = 6.45387271733267880336e+02 // 0x40842B1921EC2868
sa5 = 4.29008140027567833386e+02 // 0x407AD02157700314
sa6 = 1.08635005541779435134e+02 // 0x405B28A3EE48AE2C
sa7 = 6.57024977031928170135e+00 // 0x401A47EF8E484A93
sa8 = -6.04244152148580987438e-02 // 0xBFAEEFF2EE749A62
// Coefficients for approximation to erfc in [1/.35, 28]
rb0 = -9.86494292470009928597e-03 // 0xBF84341239E86F4A
rb1 = -7.99283237680523006574e-01 // 0xBFE993BA70C285DE
rb2 = -1.77579549177547519889e+01 // 0xC031C209555F995A
rb3 = -1.60636384855821916062e+02 // 0xC064145D43C5ED98
rb4 = -6.37566443368389627722e+02 // 0xC083EC881375F228
rb5 = -1.02509513161107724954e+03 // 0xC09004616A2E5992
rb6 = -4.83519191608651397019e+02 // 0xC07E384E9BDC383F
sb1 = 3.03380607434824582924e+01 // 0x403E568B261D5190
sb2 = 3.25792512996573918826e+02 // 0x40745CAE221B9F0A
sb3 = 1.53672958608443695994e+03 // 0x409802EB189D5118
sb4 = 3.19985821950859553908e+03 // 0x40A8FFB7688C246A
sb5 = 2.55305040643316442583e+03 // 0x40A3F219CEDF3BE6
sb6 = 4.74528541206955367215e+02 // 0x407DA874E79FE763
sb7 = -2.24409524465858183362e+01 // 0xC03670E242712D62
)
// Erf returns the error function of x.
//
// Special cases are:
// Erf(+Inf) = 1
// Erf(-Inf) = -1
// Erf(NaN) = NaN
func Erf(x float64) float64
func erf(x float64) float64 {
const (
VeryTiny = 2.848094538889218e-306 // 0x0080000000000000
Small = 1.0 / (1 << 28) // 2**-28
)
// special cases
switch {
case IsNaN(x):
return NaN()
case IsInf(x, 1):
return 1
case IsInf(x, -1):
return -1
}
sign := false
if x < 0 {
x = -x
sign = true
}
if x < 0.84375 { // |x| < 0.84375
var temp float64
if x < Small { // |x| < 2**-28
if x < VeryTiny {
temp = 0.125 * (8.0*x + efx8*x) // avoid underflow
} else {
temp = x + efx*x
}
} else {
z := x * x
r := pp0 + z*(pp1+z*(pp2+z*(pp3+z*pp4)))
s := 1 + z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))))
y := r / s
temp = x + x*y
}
if sign {
return -temp
}
return temp
}
if x < 1.25 { // 0.84375 <= |x| < 1.25
s := x - 1
P := pa0 + s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))))
Q := 1 + s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))))
if sign {
return -erx - P/Q
}
return erx + P/Q
}
if x >= 6 { // inf > |x| >= 6
if sign {
return -1
}
return 1
}
s := 1 / (x * x)
var R, S float64
if x < 1/0.35 { // |x| < 1 / 0.35 ~ 2.857143
R = ra0 + s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))))
S = 1 + s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))))
} else { // |x| >= 1 / 0.35 ~ 2.857143
R = rb0 + s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))))
S = 1 + s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))))
}
z := Float64frombits(Float64bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precision x
r := Exp(-z*z-0.5625) * Exp((z-x)*(z+x)+R/S)
if sign {
return r/x - 1
}
return 1 - r/x
}
// Erfc returns the complementary error function of x.
//
// Special cases are:
// Erfc(+Inf) = 0
// Erfc(-Inf) = 2
// Erfc(NaN) = NaN
func Erfc(x float64) float64
func erfc(x float64) float64 {
const Tiny = 1.0 / (1 << 56) // 2**-56
// special cases
switch {
case IsNaN(x):
return NaN()
case IsInf(x, 1):
return 0
case IsInf(x, -1):
return 2
}
sign := false
if x < 0 {
x = -x
sign = true
}
if x < 0.84375 { // |x| < 0.84375
var temp float64
if x < Tiny { // |x| < 2**-56
temp = x
} else {
z := x * x
r := pp0 + z*(pp1+z*(pp2+z*(pp3+z*pp4)))
s := 1 + z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))))
y := r / s
if x < 0.25 { // |x| < 1/4
temp = x + x*y
} else {
temp = 0.5 + (x*y + (x - 0.5))
}
}
if sign {
return 1 + temp
}
return 1 - temp
}
if x < 1.25 { // 0.84375 <= |x| < 1.25
s := x - 1
P := pa0 + s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))))
Q := 1 + s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))))
if sign {
return 1 + erx + P/Q
}
return 1 - erx - P/Q
}
if x < 28 { // |x| < 28
s := 1 / (x * x)
var R, S float64
if x < 1/0.35 { // |x| < 1 / 0.35 ~ 2.857143
R = ra0 + s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(ra5+s*(ra6+s*ra7))))))
S = 1 + s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(sa5+s*(sa6+s*(sa7+s*sa8)))))))
} else { // |x| >= 1 / 0.35 ~ 2.857143
if sign && x > 6 {
return 2 // x < -6
}
R = rb0 + s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(rb5+s*rb6)))))
S = 1 + s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))))
}
z := Float64frombits(Float64bits(x) & 0xffffffff00000000) // pseudo-single (20-bit) precision x
r := Exp(-z*z-0.5625) * Exp((z-x)*(z+x)+R/S)
if sign {
return 2 - r/x
}
return r / x
}
if sign {
return 2
}
return 0
}