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 | |
86 | 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 | |
97 | 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). |
133 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
134 | pub fn erff(x: f32) -> f32 { |
135 | let r: f32; |
136 | let s: f32; |
137 | let z: f32; |
138 | let y: f32; |
139 | let mut ix: u32; |
140 | let sign: usize; |
141 | |
142 | ix = x.to_bits(); |
143 | sign = (ix >> 31) as usize; |
144 | ix &= 0x7fffffff; |
145 | if ix >= 0x7f800000 { |
146 | /* erf(nan)=nan, erf(+-inf)=+-1 */ |
147 | return 1.0 - 2.0 * (sign as f32) + 1.0 / x; |
148 | } |
149 | if ix < 0x3f580000 { |
150 | /* |x| < 0.84375 */ |
151 | if ix < 0x31800000 { |
152 | /* |x| < 2**-28 */ |
153 | /*avoid underflow */ |
154 | return 0.125 * (8.0 * x + EFX8 * x); |
155 | } |
156 | z = x * x; |
157 | r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4))); |
158 | s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5)))); |
159 | y = r / s; |
160 | return x + x * y; |
161 | } |
162 | if ix < 0x40c00000 { |
163 | /* |x| < 6 */ |
164 | y = 1.0 - erfc2(ix, x); |
165 | } else { |
166 | let x1p_120 = f32::from_bits(0x03800000); |
167 | y = 1.0 - x1p_120; |
168 | } |
169 | |
170 | if sign != 0 { |
171 | -y |
172 | } else { |
173 | y |
174 | } |
175 | } |
176 | |
177 | /// Complementary error function (f32) |
178 | /// |
179 | /// Calculates the complementary probability. |
180 | /// Is `1 - erf(x)`. Is computed directly, so that you can use it to avoid |
181 | /// the loss of precision that would result from subtracting |
182 | /// large probabilities (on large `x`) from 1. |
183 | pub fn erfcf(x: f32) -> f32 { |
184 | let r: f32; |
185 | let s: f32; |
186 | let z: f32; |
187 | let y: f32; |
188 | let mut ix: u32; |
189 | let sign: usize; |
190 | |
191 | ix = x.to_bits(); |
192 | sign = (ix >> 31) as usize; |
193 | ix &= 0x7fffffff; |
194 | if ix >= 0x7f800000 { |
195 | /* erfc(nan)=nan, erfc(+-inf)=0,2 */ |
196 | return 2.0 * (sign as f32) + 1.0 / x; |
197 | } |
198 | |
199 | if ix < 0x3f580000 { |
200 | /* |x| < 0.84375 */ |
201 | if ix < 0x23800000 { |
202 | /* |x| < 2**-56 */ |
203 | return 1.0 - x; |
204 | } |
205 | z = x * x; |
206 | r = PP0 + z * (PP1 + z * (PP2 + z * (PP3 + z * PP4))); |
207 | s = 1.0 + z * (QQ1 + z * (QQ2 + z * (QQ3 + z * (QQ4 + z * QQ5)))); |
208 | y = r / s; |
209 | if sign != 0 || ix < 0x3e800000 { |
210 | /* x < 1/4 */ |
211 | return 1.0 - (x + x * y); |
212 | } |
213 | return 0.5 - (x - 0.5 + x * y); |
214 | } |
215 | if ix < 0x41e00000 { |
216 | /* |x| < 28 */ |
217 | if sign != 0 { |
218 | return 2.0 - erfc2(ix, x); |
219 | } else { |
220 | return erfc2(ix, x); |
221 | } |
222 | } |
223 | |
224 | let x1p_120 = f32::from_bits(0x03800000); |
225 | if sign != 0 { |
226 | 2.0 - x1p_120 |
227 | } else { |
228 | x1p_120 * x1p_120 |
229 | } |
230 | } |
231 | |