core_math_rs/f16/
erfc.rs

1/* Correctly-rounded complementary error function for the binary16 format
2
3Copyright (c) 2023-2025 Alexei Sibidanov and Paul Zimmermann
4
5This file is part of the CORE-MATH project
6(https://core-math.gitlabpages.inria.fr/).
7
8Permission is hereby granted, free of charge, to any person obtaining a copy
9of this software and associated documentation files (the "Software"), to deal
10in the Software without restriction, including without limitation the rights
11to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12copies of the Software, and to permit persons to whom the Software is
13furnished to do so, subject to the following conditions:
14
15The above copyright notice and this permission notice shall be included in all
16copies or substantial portions of the Software.
17
18THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24SOFTWARE.
25*/
26
27#![expect(clippy::approx_constant)]
28
29use hexf::hexf32 as h32;
30use hexf::hexf64 as h64;
31
32/* This code is based on the binary32 code erfcf.c:
33at input we convert the inputs to _Float16 (exactly),
34we then use the same code than for erfcf,
35and at output we round to _Float16.
36The changes with respect to erfcf.c are marked with a comment
37"specific erfcf16". */
38
39const E: [f64; 128] = [
40    h64!("0x1p+0"),
41    h64!("0x1.0163da9fb3335p+0"),
42    h64!("0x1.02c9a3e778061p+0"),
43    h64!("0x1.04315e86e7f85p+0"),
44    h64!("0x1.059b0d3158574p+0"),
45    h64!("0x1.0706b29ddf6dep+0"),
46    h64!("0x1.0874518759bc8p+0"),
47    h64!("0x1.09e3ecac6f383p+0"),
48    h64!("0x1.0b5586cf9890fp+0"),
49    h64!("0x1.0cc922b7247f7p+0"),
50    h64!("0x1.0e3ec32d3d1a2p+0"),
51    h64!("0x1.0fb66affed31bp+0"),
52    h64!("0x1.11301d0125b51p+0"),
53    h64!("0x1.12abdc06c31ccp+0"),
54    h64!("0x1.1429aaea92dep+0"),
55    h64!("0x1.15a98c8a58e51p+0"),
56    h64!("0x1.172b83c7d517bp+0"),
57    h64!("0x1.18af9388c8deap+0"),
58    h64!("0x1.1a35beb6fcb75p+0"),
59    h64!("0x1.1bbe084045cd4p+0"),
60    h64!("0x1.1d4873168b9aap+0"),
61    h64!("0x1.1ed5022fcd91dp+0"),
62    h64!("0x1.2063b88628cd6p+0"),
63    h64!("0x1.21f49917ddc96p+0"),
64    h64!("0x1.2387a6e756238p+0"),
65    h64!("0x1.251ce4fb2a63fp+0"),
66    h64!("0x1.26b4565e27cddp+0"),
67    h64!("0x1.284dfe1f56381p+0"),
68    h64!("0x1.29e9df51fdee1p+0"),
69    h64!("0x1.2b87fd0dad99p+0"),
70    h64!("0x1.2d285a6e4030bp+0"),
71    h64!("0x1.2ecafa93e2f56p+0"),
72    h64!("0x1.306fe0a31b715p+0"),
73    h64!("0x1.32170fc4cd831p+0"),
74    h64!("0x1.33c08b26416ffp+0"),
75    h64!("0x1.356c55f929ff1p+0"),
76    h64!("0x1.371a7373aa9cbp+0"),
77    h64!("0x1.38cae6d05d866p+0"),
78    h64!("0x1.3a7db34e59ff7p+0"),
79    h64!("0x1.3c32dc313a8e5p+0"),
80    h64!("0x1.3dea64c123422p+0"),
81    h64!("0x1.3fa4504ac801cp+0"),
82    h64!("0x1.4160a21f72e2ap+0"),
83    h64!("0x1.431f5d950a897p+0"),
84    h64!("0x1.44e086061892dp+0"),
85    h64!("0x1.46a41ed1d0057p+0"),
86    h64!("0x1.486a2b5c13cdp+0"),
87    h64!("0x1.4a32af0d7d3dep+0"),
88    h64!("0x1.4bfdad5362a27p+0"),
89    h64!("0x1.4dcb299fddd0dp+0"),
90    h64!("0x1.4f9b2769d2ca7p+0"),
91    h64!("0x1.516daa2cf6642p+0"),
92    h64!("0x1.5342b569d4f82p+0"),
93    h64!("0x1.551a4ca5d920fp+0"),
94    h64!("0x1.56f4736b527dap+0"),
95    h64!("0x1.58d12d497c7fdp+0"),
96    h64!("0x1.5ab07dd485429p+0"),
97    h64!("0x1.5c9268a5946b7p+0"),
98    h64!("0x1.5e76f15ad2148p+0"),
99    h64!("0x1.605e1b976dc09p+0"),
100    h64!("0x1.6247eb03a5585p+0"),
101    h64!("0x1.6434634ccc32p+0"),
102    h64!("0x1.6623882552225p+0"),
103    h64!("0x1.68155d44ca973p+0"),
104    h64!("0x1.6a09e667f3bcdp+0"),
105    h64!("0x1.6c012750bdabfp+0"),
106    h64!("0x1.6dfb23c651a2fp+0"),
107    h64!("0x1.6ff7df9519484p+0"),
108    h64!("0x1.71f75e8ec5f74p+0"),
109    h64!("0x1.73f9a48a58174p+0"),
110    h64!("0x1.75feb564267c9p+0"),
111    h64!("0x1.780694fde5d3fp+0"),
112    h64!("0x1.7a11473eb0187p+0"),
113    h64!("0x1.7c1ed0130c132p+0"),
114    h64!("0x1.7e2f336cf4e62p+0"),
115    h64!("0x1.80427543e1a12p+0"),
116    h64!("0x1.82589994cce13p+0"),
117    h64!("0x1.8471a4623c7adp+0"),
118    h64!("0x1.868d99b4492edp+0"),
119    h64!("0x1.88ac7d98a6699p+0"),
120    h64!("0x1.8ace5422aa0dbp+0"),
121    h64!("0x1.8cf3216b5448cp+0"),
122    h64!("0x1.8f1ae99157736p+0"),
123    h64!("0x1.9145b0b91ffc6p+0"),
124    h64!("0x1.93737b0cdc5e5p+0"),
125    h64!("0x1.95a44cbc8520fp+0"),
126    h64!("0x1.97d829fde4e5p+0"),
127    h64!("0x1.9a0f170ca07bap+0"),
128    h64!("0x1.9c49182a3f09p+0"),
129    h64!("0x1.9e86319e32323p+0"),
130    h64!("0x1.a0c667b5de565p+0"),
131    h64!("0x1.a309bec4a2d33p+0"),
132    h64!("0x1.a5503b23e255dp+0"),
133    h64!("0x1.a799e1330b358p+0"),
134    h64!("0x1.a9e6b5579fdbfp+0"),
135    h64!("0x1.ac36bbfd3f37ap+0"),
136    h64!("0x1.ae89f995ad3adp+0"),
137    h64!("0x1.b0e07298db666p+0"),
138    h64!("0x1.b33a2b84f15fbp+0"),
139    h64!("0x1.b59728de5593ap+0"),
140    h64!("0x1.b7f76f2fb5e47p+0"),
141    h64!("0x1.ba5b030a1064ap+0"),
142    h64!("0x1.bcc1e904bc1d2p+0"),
143    h64!("0x1.bf2c25bd71e09p+0"),
144    h64!("0x1.c199bdd85529cp+0"),
145    h64!("0x1.c40ab5fffd07ap+0"),
146    h64!("0x1.c67f12e57d14bp+0"),
147    h64!("0x1.c8f6d9406e7b5p+0"),
148    h64!("0x1.cb720dcef9069p+0"),
149    h64!("0x1.cdf0b555dc3fap+0"),
150    h64!("0x1.d072d4a07897cp+0"),
151    h64!("0x1.d2f87080d89f2p+0"),
152    h64!("0x1.d5818dcfba487p+0"),
153    h64!("0x1.d80e316c98398p+0"),
154    h64!("0x1.da9e603db3285p+0"),
155    h64!("0x1.dd321f301b46p+0"),
156    h64!("0x1.dfc97337b9b5fp+0"),
157    h64!("0x1.e264614f5a129p+0"),
158    h64!("0x1.e502ee78b3ff6p+0"),
159    h64!("0x1.e7a51fbc74c83p+0"),
160    h64!("0x1.ea4afa2a490dap+0"),
161    h64!("0x1.ecf482d8e67f1p+0"),
162    h64!("0x1.efa1bee615a27p+0"),
163    h64!("0x1.f252b376bba97p+0"),
164    h64!("0x1.f50765b6e454p+0"),
165    h64!("0x1.f7bfdad9cbe14p+0"),
166    h64!("0x1.fa7c1819e90d8p+0"),
167    h64!("0x1.fd3c22b8f71f1p+0"),
168];
169
170// specific erfcf16: input renamed to xf16
171pub fn cr_erfcf16(xf16: f16) -> f16 {
172    let x = f32::from(xf16);
173    let ax = x.abs();
174    let axd = f64::from(ax);
175    let x2 = axd * axd;
176    let t = x;
177    let at = t.to_bits() & ((!0) >> 1);
178    let sgn = t.to_bits() >> 31;
179    let i = at > 0x40051000;
180    /* for x < -0x1.ea8f94p+1, erfc(x) rounds to 2 (to nearest) */
181    if unlikely(t.to_bits() > 0xc07547ca) {
182        // x < -0x1.ea8f94p+1
183        if unlikely(t.to_bits() >= 0xff800000) {
184            // -Inf or NaN
185            if t.to_bits() == 0xff800000 {
186                return 2.0; // -Inf
187            }
188            return (x + x) as f16; // NaN
189        }
190        return (2.0 - h32!("0x1p-25")) as f16; // rounds to 2 or nextbelow(2)
191    }
192    /* at is the absolute value of x
193    for x >= 0x1.41bbf8p+3, erfc(x) < 2^-150, thus rounds to 0 or to 2^-149
194    depending on the rounding mode */
195    if unlikely(at >= 0x4120ddfc) {
196        // |x| >= 0x1.41bbf8p+3
197        if unlikely(at >= 0x7f800000) {
198            // +Inf or NaN
199            if at == 0x7f800000 {
200                return 0.0; // +Inf
201            }
202            return (x + x) as f16; // NaN
203        }
204        // 0x1p-149f * 0.25f rounds to 0 or 2^-149 depending on rounding
205        return (h32!("0x1p-149") * 0.25) as f16;
206    }
207    if unlikely(at <= 0x3db80000) {
208        // |x| <= 0x1.7p-4
209        if unlikely(t.to_bits() == 0xb76c9f62) {
210            // x = -0x1.d93ec4p-17
211            return (h32!("0x1.00010ap+0") + h32!("0x1p-25")) as f16; // exceptional case
212        }
213        /* for |x| <= 0x1.c5bf88p-26. erfc(x) rounds to 1 (to nearest) */
214        if unlikely(at <= 0x32e2dfc4) {
215            // |x| <= 0x1.c5bf88p-26
216            if unlikely(at == 0) {
217                return 1.0;
218            };
219            const D: [f32; 2] = [h32!("-0x1p-26"), h32!("0x1p-25")];
220            return (1.0 + D[sgn as usize]) as f16;
221        }
222        /* around 0, erfc(x) behaves as 1 - (odd polynomial) */
223        const C: [f64; 5] = [
224            h64!("0x1.20dd750429b6dp+0"),
225            h64!("-0x1.812746b03610bp-2"),
226            h64!("0x1.ce2f218831d2fp-4"),
227            h64!("-0x1.b82c609607dcbp-6"),
228            h64!("0x1.553af09b8008ep-8"),
229        ];
230        let f0 = f64::from(x) * (C[0] + x2 * (C[1] + x2 * (C[2] + x2 * (C[3] + x2 * (C[4])))));
231        return (1.0 - f0) as f16;
232    }
233
234    /* now -0x1.ea8f94p+1 <= x <= 0x1.41bbf8p+3, with |x| > 0x1.7p-4 */
235    const ILN2: f64 = h64!("0x1.71547652b82fep+0");
236    const LN2H: f64 = h64!("0x1.62e42fefap-8");
237    const LN2L: f64 = h64!("0x1.cf79abd6f5dc8p-47");
238    let jt = f64::mul_add(x2, ILN2, -(1024.0 + h64!("0x1p-8")));
239    let j = ((jt.to_bits() << 12) as i64) >> 48;
240    let s_ = (((j >> 7) + (0x3ff | i64::from(sgn) << 11)) as u64) << 52;
241    const CH: [f64; 4] = [
242        h64!("-0x1.ffffffffff333p-2"),
243        h64!("0x1.5555555556a14p-3"),
244        h64!("-0x1.55556666659b4p-5"),
245        h64!("0x1.1111074cc7b22p-7"),
246    ];
247    let d = (x2 + LN2H * (j as f64)) + LN2L * (j as f64);
248    let d2 = d * d;
249    let e0 = E[(j & 127) as usize];
250    let f = d + d2 * ((CH[0] + d * CH[1]) + d2 * (CH[2] + d * CH[3]));
251    const CT: [[f64; 16]; 2] = [
252        [
253            h64!("0x1.c162355429b28p-1"),
254            h64!("0x1.d99999999999ap+1"),
255            h64!("0x1.da951cece2b85p-2"),
256            h64!("-0x1.70ef6cff4bcc4p+0"),
257            h64!("0x1.3d7f7b3d617dep+1"),
258            h64!("-0x1.9d0aa47537c51p+1"),
259            h64!("0x1.9754ea9a3fcb1p+1"),
260            h64!("-0x1.27a5453fcc015p+1"),
261            h64!("0x1.1ef2e0531aebap+0"),
262            h64!("-0x1.eca090f5a1c06p-3"),
263            h64!("-0x1.7a3cd173a063cp-4"),
264            h64!("0x1.30fa68a68fdddp-4"),
265            h64!("0x1.55ad9a326993ap-10"),
266            h64!("-0x1.07e7b0bb39fbfp-6"),
267            h64!("0x1.2328706c0e95p-10"),
268            h64!("0x1.d6aa0b7b19cfep-9"),
269        ],
270        [
271            h64!("0x1.137c8983f8516p+2"),
272            h64!("0x1.799999999999ap+1"),
273            h64!("0x1.05b53aa241333p-3"),
274            h64!("-0x1.a3f53872bf87p-3"),
275            h64!("0x1.de4c30742c9d5p-4"),
276            h64!("-0x1.cb24bfa591986p-5"),
277            h64!("0x1.666aec059ca5fp-6"),
278            h64!("-0x1.a61250eb26b0bp-8"),
279            h64!("0x1.2b28b7924b34dp-10"),
280            h64!("0x1.41b13a9d45013p-15"),
281            h64!("-0x1.6dd5e8a273613p-14"),
282            h64!("0x1.09ce8ea5e8da5p-16"),
283            h64!("0x1.33923b4102981p-18"),
284            h64!("-0x1.1dfd161e3f984p-19"),
285            h64!("-0x1.c87618fcae3b3p-23"),
286            h64!("0x1.e8a6ffa0ba2c7p-23"),
287        ],
288    ];
289    let z = (axd - CT[usize::from(i)][0]) / (axd + CT[usize::from(i)][1]);
290    let z2 = z * z;
291    let z4 = z2 * z2;
292    let z8 = z4 * z4;
293    let c = &CT[usize::from(i)][3..];
294    let mut s = (((c[0] + z * c[1]) + z2 * (c[2] + z * c[3]))
295        + z4 * ((c[4] + z * c[5]) + z2 * (c[6] + z * c[7])))
296        + z8 * (((c[8] + z * c[9]) + z2 * (c[10] + z * c[11])) + z4 * (c[12]));
297    s = CT[usize::from(i)][2] + z * s;
298    const OFF: [f64; 2] = [0.0, 2.0];
299    let r = (f64::from_bits(s_) * (e0 - f * e0)) * s;
300    let y = OFF[sgn as usize] + r;
301    y as f16
302}
303
304// TODO: replace with std::hint::unlikely
305const fn unlikely(b: bool) -> bool {
306    b
307}
308
309#[cfg(test)]
310mod tests {
311    #[test]
312    fn exhaustive() {
313        for b in 0..=u16::MAX {
314            let x = f16::from_bits(b);
315            let y1 = super::cr_erfcf16(x);
316            let y2 = core_math::erfcf16(x);
317            assert_eq!(
318                y1.to_bits(),
319                y2.to_bits(),
320                "erfcf16({x} @ {b:#04x}) = ({y1} @ {y1b:#04x}) vs ({y2} @ {y2b:#04x})",
321                y1b = y1.to_bits(),
322                y2b = y2.to_bits()
323            );
324        }
325    }
326
327    #[test]
328    fn edge() {
329        assert_eq!(
330            super::cr_erfcf16(-f16::NAN).to_bits(),
331            (-f16::NAN).to_bits()
332        );
333        assert_eq!(super::cr_erfcf16(-f16::INFINITY), 2.0);
334        assert_eq!(super::cr_erfcf16(-0.0), 1.0);
335        assert_eq!(super::cr_erfcf16(0.0), 1.0);
336        assert_eq!(
337            super::cr_erfcf16(f16::INFINITY).to_bits(),
338            (0.0_f16).to_bits()
339        );
340        assert_eq!(super::cr_erfcf16(f16::NAN).to_bits(), (f16::NAN).to_bits());
341    }
342}
core_math_rs/f16/erfc.rs

core_math_rs/f16/
erfc.rs
