1 /* origin: FreeBSD /usr/src/lib/msun/src/s_erff.c */
2 /*
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
4 */
5 /*
6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice
12 * is preserved.
13 * ====================================================
14 */
15
16 use super::{expf, fabsf};
17
18 const ERX: f32 = 8.4506291151e-01; /* 0x3f58560b */
19 /*
20 * Coefficients for approximation to erf on [0,0.84375]
21 */
22 const EFX8: f32 = 1.0270333290e+00; /* 0x3f8375d4 */
23 const PP0: f32 = 1.2837916613e-01; /* 0x3e0375d4 */
24 const PP1: f32 = -3.2504209876e-01; /* 0xbea66beb */
25 const PP2: f32 = -2.8481749818e-02; /* 0xbce9528f */
26 const PP3: f32 = -5.7702702470e-03; /* 0xbbbd1489 */
27 const PP4: f32 = -2.3763017452e-05; /* 0xb7c756b1 */
28 const QQ1: f32 = 3.9791721106e-01; /* 0x3ecbbbce */
29 const QQ2: f32 = 6.5022252500e-02; /* 0x3d852a63 */
30 const QQ3: f32 = 5.0813062117e-03; /* 0x3ba68116 */
31 const QQ4: f32 = 1.3249473704e-04; /* 0x390aee49 */
32 const QQ5: f32 = -3.9602282413e-06; /* 0xb684e21a */
33 /*
34 * Coefficients for approximation to erf in [0.84375,1.25]
35 */
36 const PA0: f32 = -2.3621185683e-03; /* 0xbb1acdc6 */
37 const PA1: f32 = 4.1485610604e-01; /* 0x3ed46805 */
38 const PA2: f32 = -3.7220788002e-01; /* 0xbebe9208 */
39 const PA3: f32 = 3.1834661961e-01; /* 0x3ea2fe54 */
40 const PA4: f32 = -1.1089469492e-01; /* 0xbde31cc2 */
41 const PA5: f32 = 3.5478305072e-02; /* 0x3d1151b3 */
42 const PA6: f32 = -2.1663755178e-03; /* 0xbb0df9c0 */
43 const QA1: f32 = 1.0642088205e-01; /* 0x3dd9f331 */
44 const QA2: f32 = 5.4039794207e-01; /* 0x3f0a5785 */
45 const QA3: f32 = 7.1828655899e-02; /* 0x3d931ae7 */
46 const QA4: f32 = 1.2617121637e-01; /* 0x3e013307 */
47 const QA5: f32 = 1.3637083583e-02; /* 0x3c5f6e13 */
48 const QA6: f32 = 1.1984500103e-02; /* 0x3c445aa3 */
49 /*
50 * Coefficients for approximation to erfc in [1.25,1/0.35]
51 */
52 const RA0: f32 = -9.8649440333e-03; /* 0xbc21a093 */
53 const RA1: f32 = -6.9385856390e-01; /* 0xbf31a0b7 */
54 const RA2: f32 = -1.0558626175e+01; /* 0xc128f022 */
55 const RA3: f32 = -6.2375331879e+01; /* 0xc2798057 */
56 const RA4: f32 = -1.6239666748e+02; /* 0xc322658c */
57 const RA5: f32 = -1.8460508728e+02; /* 0xc3389ae7 */
58 const RA6: f32 = -8.1287437439e+01; /* 0xc2a2932b */
59 const RA7: f32 = -9.8143291473e+00; /* 0xc11d077e */
60 const SA1: f32 = 1.9651271820e+01; /* 0x419d35ce */
61 const SA2: f32 = 1.3765776062e+02; /* 0x4309a863 */
62 const SA3: f32 = 4.3456588745e+02; /* 0x43d9486f */
63 const SA4: f32 = 6.4538726807e+02; /* 0x442158c9 */
64 const SA5: f32 = 4.2900814819e+02; /* 0x43d6810b */
65 const SA6: f32 = 1.0863500214e+02; /* 0x42d9451f */
66 const SA7: f32 = 6.5702495575e+00; /* 0x40d23f7c */
67 const SA8: f32 = -6.0424413532e-02; /* 0xbd777f97 */
68 /*
69 * Coefficients for approximation to erfc in [1/.35,28]
70 */
71 const RB0: f32 = -9.8649431020e-03; /* 0xbc21a092 */
72 const RB1: f32 = -7.9928326607e-01; /* 0xbf4c9dd4 */
73 const RB2: f32 = -1.7757955551e+01; /* 0xc18e104b */
74 const RB3: f32 = -1.6063638306e+02; /* 0xc320a2ea */
75 const RB4: f32 = -6.3756646729e+02; /* 0xc41f6441 */
76 const RB5: f32 = -1.0250950928e+03; /* 0xc480230b */
77 const RB6: f32 = -4.8351919556e+02; /* 0xc3f1c275 */
78 const SB1: f32 = 3.0338060379e+01; /* 0x41f2b459 */
79 const SB2: f32 = 3.2579251099e+02; /* 0x43a2e571 */
80 const SB3: f32 = 1.5367296143e+03; /* 0x44c01759 */
81 const SB4: f32 = 3.1998581543e+03; /* 0x4547fdbb */
82 const SB5: f32 = 2.5530502930e+03; /* 0x451f90ce */
83 const SB6: f32 = 4.7452853394e+02; /* 0x43ed43a7 */
84 const SB7: f32 = -2.2440952301e+01; /* 0xc1b38712 */
85
erfc1(x: f32) -> f3286 fn erfc1(x: f32) -> f32 {
87 let s: f32;
88 let p: f32;
89 let q: f32;
90
91 s = fabsf(x) - 1.0;
92 p = PA0 + s * (PA1 + s * (PA2 + s * (PA3 + s * (PA4 + s * (PA5 + s * PA6)))));
93 q = 1.0 + s * (QA1 + s * (QA2 + s * (QA3 + s * (QA4 + s * (QA5 + s * QA6)))));
94 return 1.0 - ERX - p / q;
95 }
96
erfc2(mut ix: u32, mut x: f32) -> f3297 fn erfc2(mut ix: u32, mut x: f32) -> f32 {
98 let s: f32;
99 let r: f32;
100 let big_s: f32;
101 let z: f32;
102
103 if ix < 0x3fa00000 {
104 /* |x| < 1.25 */
105 return erfc1(x);
106 }
107
108 x = fabsf(x);
109 s = 1.0 / (x * x);
110 if ix < 0x4036db6d {
111 /* |x| < 1/0.35 */
112 r = RA0 + s * (RA1 + s * (RA2 + s * (RA3 + s * (RA4 + s * (RA5 + s * (RA6 + s * RA7))))));
113 big_s = 1.0
114 + s * (SA1
115 + s * (SA2 + s * (SA3 + s * (SA4 + s * (SA5 + s * (SA6 + s * (SA7 + s * SA8)))))));
116 } else {
117 /* |x| >= 1/0.35 */
118 r = RB0 + s * (RB1 + s * (RB2 + s * (RB3 + s * (RB4 + s * (RB5 + s * RB6)))));
119 big_s =
120 1.0 + s * (SB1 + s * (SB2 + s * (SB3 + s * (SB4 + s * (SB5 + s * (SB6 + s * SB7))))));
121 }
122 ix = x.to_bits();
123 z = f32::from_bits(ix & 0xffffe000);
124
125 expf(-z * z - 0.5625) * expf((z - x) * (z + x) + r / big_s) / x
126 }
127
128 /// Error function (f32)
129 ///
130 /// Calculates an approximation to the “error function”, which estimates
131 /// the probability that an observation will fall within x standard
132 /// deviations of the mean (assuming a normal distribution).
erff(x: f32) -> f32133 pub fn erff(x: f32) -> f32 {
134 let r: f32;
135 let s: f32;
136 let z: f32;
137 let y: f32;
138 let mut ix: u32;
139 let sign: usize;
140
141 ix = x.to_bits();
142 sign = (ix >> 31) as usize;
143 ix &= 0x7fffffff;
144 if ix >= 0x7f800000 {
145 /* erf(nan)=nan, erf(+-inf)=+-1 */
146 return 1.0 - 2.0 * (sign as f32) + 1.0 / x;
147 }
148 if ix < 0x3f580000 {
149 /* |x| < 0.84375 */
150 if ix < 0x31800000 {
151 /* |x| < 2**-28 */
152 /*avoid underflow */
153 return 0.125 * (8.0 * x + EFX8 * x);
154 }
155 z = x * x;
156 r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4)));
157 s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5))));
158 y = r / s;
159 return x + x * y;
160 }
161 if ix < 0x40c00000 {
162 /* |x| < 6 */
163 y = 1.0 - erfc2(ix, x);
164 } else {
165 let x1p_120 = f32::from_bits(0x03800000);
166 y = 1.0 - x1p_120;
167 }
168
169 if sign != 0 {
170 -y
171 } else {
172 y
173 }
174 }
175
176 /// Error function (f32)
177 ///
178 /// Calculates the complementary probability.
179 /// Is `1 - erf(x)`. Is computed directly, so that you can use it to avoid
180 /// the loss of precision that would result from subtracting
181 /// large probabilities (on large `x`) from 1.
erfcf(x: f32) -> f32182 pub fn erfcf(x: f32) -> f32 {
183 let r: f32;
184 let s: f32;
185 let z: f32;
186 let y: f32;
187 let mut ix: u32;
188 let sign: usize;
189
190 ix = x.to_bits();
191 sign = (ix >> 31) as usize;
192 ix &= 0x7fffffff;
193 if ix >= 0x7f800000 {
194 /* erfc(nan)=nan, erfc(+-inf)=0,2 */
195 return 2.0 * (sign as f32) + 1.0 / x;
196 }
197
198 if ix < 0x3f580000 {
199 /* |x| < 0.84375 */
200 if ix < 0x23800000 {
201 /* |x| < 2**-56 */
202 return 1.0 - x;
203 }
204 z = x * x;
205 r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4)));
206 s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5))));
207 y = r / s;
208 if sign != 0 || ix < 0x3e800000 {
209 /* x < 1/4 */
210 return 1.0 - (x + x * y);
211 }
212 return 0.5 - (x - 0.5 + x * y);
213 }
214 if ix < 0x41e00000 {
215 /* |x| < 28 */
216 if sign != 0 {
217 return 2.0 - erfc2(ix, x);
218 } else {
219 return erfc2(ix, x);
220 }
221 }
222
223 let x1p_120 = f32::from_bits(0x03800000);
224 if sign != 0 {
225 2.0 - x1p_120
226 } else {
227 x1p_120 * x1p_120
228 }
229 }
230