core_math_rs/f16/
erf.rs

1/* Correctly-rounded error function for binary16 value.
2
3Copyright (c) 2022-2025 Alexei Sibidanov and Paul Zimmermann
4
5This file is ported from 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 erff.c:
33at input we convert the inputs to _Float16 (exactly),
34we then use the same code than for erff,
35and at output we round to _Float16.
36The changes with respect to erff.c are marked with a comment
37"specific erff16". */
38
39// specific erff16: input renamed to xf16
40pub fn cr_erff16(xf16: f16) -> f16 {
41    let x = f32::from(xf16);
42    /* for 7 <= i < 63, C[i-7] is a degree-7 polynomial approximation of
43    erf(i/16+1/32+x) for -1/32 <= x <= 1/32 */
44    const C: [[f64; 8]; 56] = [
45        [
46            h64!("0x1.f86faa9428f9cp-2"),
47            h64!("0x1.cfc41e36c7dfap-1"),
48            h64!("-0x1.b2c7dc53508b9p-2"),
49            h64!("-0x1.5a9de93fa556ep-3"),
50            h64!("0x1.731793dbb01b5p-3"),
51            h64!("0x1.133e06426cf18p-6"),
52            h64!("-0x1.a12a6289cafd8p-5"),
53            h64!("0x1.717d6f1d6f557p-9"),
54        ],
55        [
56            h64!("0x1.1855a5fd3dd50p-1"),
57            h64!("0x1.b3aafcc27502fp-1"),
58            h64!("-0x1.cee5ac8e92bb2p-2"),
59            h64!("-0x1.fa02983ca2d79p-4"),
60            h64!("0x1.77cd746cb1922p-3"),
61            h64!("-0x1.fa6f277886487p-10"),
62            h64!("-0x1.8de75458db416p-5"),
63            h64!("0x1.00899c98551c9p-7"),
64        ],
65        [
66            h64!("0x1.32a54cb8db67ap-1"),
67            h64!("0x1.96164fafd8de5p-1"),
68            h64!("-0x1.e23a7ea0c9ad3p-2"),
69            h64!("-0x1.3f5ee15671cf4p-4"),
70            h64!("0x1.70e468a3d72d9p-3"),
71            h64!("-0x1.3da68037cfc99p-6"),
72            h64!("-0x1.69ed9ba1f9839p-5"),
73            h64!("0x1.8cab9244a4ff4p-7"),
74        ],
75        [
76            h64!("0x1.4b13713ad3513p-1"),
77            h64!("0x1.7791b886e7405p-1"),
78            h64!("-0x1.ecef423109bf5p-2"),
79            h64!("-0x1.15c3c5cec6847p-5"),
80            h64!("0x1.5f688fc931ba6p-3"),
81            h64!("-0x1.1da63ed190037p-5"),
82            h64!("-0x1.38427ca63cca4p-5"),
83            h64!("0x1.fa00e52525e17p-7"),
84        ],
85        [
86            h64!("0x1.61955607dd15dp-1"),
87            h64!("0x1.58a445da7c74ep-1"),
88            h64!("-0x1.ef6c246a0f66cp-2"),
89            h64!("0x1.e83e0d9d61330p-8"),
90            h64!("0x1.44cc65535bc9fp-3"),
91            h64!("-0x1.87d3c4860435dp-5"),
92            h64!("-0x1.f90b10501169bp-6"),
93            h64!("0x1.22295856d427ap-6"),
94        ],
95        [
96            h64!("0x1.762870f720c6fp-1"),
97            h64!("0x1.39ccc1b136d5cp-1"),
98            h64!("-0x1.ea4feea4e4744p-2"),
99            h64!("0x1.715e5952ebfbap-5"),
100            h64!("0x1.22cdbd83c75c4p-3"),
101            h64!("-0x1.da50aa1d925b6p-5"),
102            h64!("-0x1.754dc0a29b4ddp-6"),
103            h64!("0x1.350b6bef9392cp-6"),
104        ],
105        [
106            h64!("0x1.88d1cd474a2e0p-1"),
107            h64!("0x1.1b7e98fe26219p-1"),
108            h64!("-0x1.de65a22ce1419p-2"),
109            h64!("0x1.40686a3f16400p-4"),
110            h64!("0x1.f6b0cbb216b2bp-4"),
111            h64!("-0x1.09c7c903edd57p-4"),
112            h64!("-0x1.da7529fde641p-7"),
113            h64!("0x1.362a7a0588eabp-6"),
114        ],
115        [
116            h64!("0x1.999d4192a5717p-1"),
117            h64!("0x1.fc3ee5d1524b3p-2"),
118            h64!("-0x1.cc990045b55c8p-2"),
119            h64!("0x1.b37338e68b37dp-4"),
120            h64!("0x1.a0d120c872ea7p-4"),
121            h64!("-0x1.19bb2b07ecff6p-4"),
122            h64!("-0x1.a110f5f593aafp-8"),
123            h64!("0x1.272c15a57720ep-6"),
124        ],
125        [
126            h64!("0x1.a89c850b7d54dp-1"),
127            h64!("0x1.c40b0729ed54ap-2"),
128            h64!("-0x1.b5eaaef0a2346p-2"),
129            h64!("0x1.0847c7dacbae1p-3"),
130            h64!("0x1.47de0ba6d18fbp-4"),
131            h64!("-0x1.1d9de77a4b648p-4"),
132            h64!("0x1.30ffbe56f0726p-10"),
133            h64!("0x1.0a9cb99feea01p-6"),
134        ],
135        [
136            h64!("0x1.b5e62fce16096p-1"),
137            h64!("0x1.8eed36b886d95p-2"),
138            h64!("-0x1.9b64a06e50705p-2"),
139            h64!("0x1.2bb6e2c744df5p-3"),
140            h64!("0x1.dee3261ca61bcp-5"),
141            h64!("-0x1.16996004f7da5p-4"),
142            h64!("0x1.fdff37bae983ep-8"),
143            h64!("0x1.c750083e65f9ap-7"),
144        ],
145        [
146            h64!("0x1.c194b1d49a184p-1"),
147            h64!("0x1.5d4fd33729015p-2"),
148            h64!("-0x1.7e0f4f045addbp-2"),
149            h64!("0x1.444bc66c31a1bp-3"),
150            h64!("0x1.356dbf16ec8f1p-5"),
151            h64!("-0x1.0643de0906cd8p-4"),
152            h64!("0x1.b281af7bd3a2cp-7"),
153            h64!("0x1.6b97eaa2c6abdp-7"),
154        ],
155        [
156            h64!("0x1.cbc54b476248ep-1"),
157            h64!("0x1.2f7cc3fe6f423p-2"),
158            h64!("-0x1.5ee8429e36de8p-2"),
159            h64!("0x1.52a8395f96177p-3"),
160            h64!("0x1.313761ba257dcp-6"),
161            h64!("-0x1.dcf844d5fed8fp-5"),
162            h64!("0x1.1e1420f475fa9p-6"),
163            h64!("0x1.091c7dc1e18b2p-7"),
164        ],
165        [
166            h64!("0x1.d4970f9ce00d9p-1"),
167            h64!("0x1.059f59af7a905p-2"),
168            h64!("-0x1.3eda354de36c3p-2"),
169            h64!("0x1.57b85ad439779p-3"),
170            h64!("0x1.8e913b9778136p-10"),
171            h64!("-0x1.a2893bd3435f4p-5"),
172            h64!("0x1.4d3a90e37164ap-6"),
173            h64!("0x1.4ce7f6e19a902p-8"),
174        ],
175        [
176            h64!("0x1.dc29fb60715b0p-1"),
177            h64!("0x1.bf8e1b1ca2277p-3"),
178            h64!("-0x1.1eb7095e5d6d2p-2"),
179            h64!("0x1.549ea6f7a64f4p-3"),
180            h64!("-0x1.b10f12f3877a3p-7"),
181            h64!("-0x1.61420c8f7156ap-5"),
182            h64!("0x1.674f1f92a8812p-6"),
183            h64!("0x1.25543ffd74d52p-9"),
184        ],
185        [
186            h64!("0x1.e29e22a89d767p-1"),
187            h64!("0x1.7bd5c7df3fe99p-3"),
188            h64!("-0x1.fe674494077bfp-3"),
189            h64!("0x1.4a9feacf86578p-3"),
190            h64!("-0x1.a008269076644p-6"),
191            h64!("-0x1.1cf0e8fb4f1cbp-5"),
192            h64!("0x1.6e0d2ef105fb3p-6"),
193            h64!("-0x1.367205fbd7876p-12"),
194        ],
195        [
196            h64!("0x1.e812fc64db36ap-1"),
197            h64!("0x1.3fda6bc016991p-3"),
198            h64!("-0x1.c1cb278627920p-3"),
199            h64!("0x1.3b10512314f1ep-3"),
200            h64!("-0x1.1e6457bb1b9a9p-5"),
201            h64!("-0x1.b1f6474e2388cp-6"),
202            h64!("0x1.640a5345f7ec7p-6"),
203            h64!("-0x1.3dae5a997fdbp-9"),
204        ],
205        [
206            h64!("0x1.eca6ccd709544p-1"),
207            h64!("0x1.0b3f52ce8c380p-3"),
208            h64!("-0x1.8885019f63c6dp-3"),
209            h64!("0x1.274275fc91a05p-3"),
210            h64!("-0x1.57f73699a8372p-5"),
211            h64!("-0x1.3076a305fc7cep-6"),
212            h64!("0x1.4c6ae04843a41p-6"),
213            h64!("-0x1.0be5fcf5ecc91p-8"),
214        ],
215        [
216            h64!("0x1.f0762fde45ee7p-1"),
217            h64!("0x1.bb1c972f23e4ap-4"),
218            h64!("-0x1.5341e3c01b58dp-3"),
219            h64!("0x1.107929f6f0b60p-3"),
220            h64!("-0x1.7e1b34f976c02p-5"),
221            h64!("-0x1.73b62589c234ap-7"),
222            h64!("0x1.2a97ee1876486p-6"),
223            h64!("-0x1.595f40a3150fep-8"),
224        ],
225        [
226            h64!("0x1.f39bc242e43e6p-1"),
227            h64!("0x1.6c7e64e7281c5p-4"),
228            h64!("-0x1.2274b86835fd3p-3"),
229            h64!("0x1.efb890e5c770dp-4"),
230            h64!("-0x1.92c7db16847e0p-5"),
231            h64!("-0x1.45477db5e2dd4p-8"),
232            h64!("0x1.01fc6165fc866p-6"),
233            h64!("-0x1.8845509030c2cp-8"),
234        ],
235        [
236            h64!("0x1.f62fe80272419p-1"),
237            h64!("0x1.297db960e4f5dp-4"),
238            h64!("-0x1.ecb83b087c04fp-4"),
239            h64!("0x1.bce18363ca3d1p-4"),
240            h64!("-0x1.985aaf776482cp-5"),
241            h64!("0x1.cd953efdae886p-12"),
242            h64!("0x1.ab9a0b89b54ffp-7"),
243            h64!("-0x1.9b5e576ccc31cp-8"),
244        ],
245        [
246            h64!("0x1.f848acb544e95p-1"),
247            h64!("0x1.e1d4cf1e24501p-5"),
248            h64!("-0x1.9e12e1fde5552p-4"),
249            h64!("0x1.8a27806df3d1bp-4"),
250            h64!("-0x1.91674e5eb3319p-5"),
251            h64!("0x1.3bc75595b2db8p-8"),
252            h64!("0x1.51bc537ac61afp-7"),
253            h64!("-0x1.96b23b19ea04dp-8"),
254        ],
255        [
256            h64!("0x1.f9f9ba8d3c733p-1"),
257            h64!("0x1.83298d7172108p-5"),
258            h64!("-0x1.58d101f905a75p-4"),
259            h64!("0x1.58f1456f8639bp-4"),
260            h64!("-0x1.808d1850b8231p-5"),
261            h64!("0x1.0c1bd99c348a7p-7"),
262            h64!("0x1.f61e9d7bc48cap-8"),
263            h64!("-0x1.7f07c13441774p-8"),
264        ],
265        [
266            h64!("0x1.fb54641aebbc9p-1"),
267            h64!("0x1.34ac36ad8dafap-5"),
268            h64!("-0x1.1c8ec267f9405p-4"),
269            h64!("0x1.2a52c5d841848p-4"),
270            h64!("-0x1.68541c02b3b6bp-5"),
271            h64!("0x1.5afe400196379p-7"),
272            h64!("0x1.565b2d6eda3d6p-8"),
273            h64!("-0x1.596aaff29e739p-8"),
274        ],
275        [
276            h64!("0x1.fc67bcf2d7b8fp-1"),
277            h64!("0x1.e85c449e377efp-6"),
278            h64!("-0x1.d177f166c07c6p-5"),
279            h64!("0x1.fe23b7584b504p-5"),
280            h64!("-0x1.4b12109613313p-5"),
281            h64!("0x1.8d9905c0acf7dp-7"),
282            h64!("0x1.9265032a669dap-9"),
283            h64!("-0x1.2ac4a6dbcbf3ep-8"),
284        ],
285        [
286            h64!("0x1.fd40bd6d7a785p-1"),
287            h64!("0x1.7f5188610ddc7p-6"),
288            h64!("-0x1.7954423f7c998p-5"),
289            h64!("0x1.af5baae33887fp-5"),
290            h64!("-0x1.2ad77c7cbc474p-5"),
291            h64!("0x1.a7b8c47ec2a51p-7"),
292            h64!("0x1.46646ee094bccp-10"),
293            h64!("-0x1.ef19d8db8673p-9"),
294        ],
295        [
296            h64!("0x1.fdea6e062d0c9p-1"),
297            h64!("0x1.2a875b5ffab58p-6"),
298            h64!("-0x1.2f3178cd6dcd5p-5"),
299            h64!("0x1.68d1c45b94182p-5"),
300            h64!("-0x1.09648ed3aeaefp-5"),
301            h64!("0x1.ad8b150d38164p-7"),
302            h64!("-0x1.e9a6023d9429fp-13"),
303            h64!("-0x1.8722d19ee2e8ep-9"),
304        ],
305        [
306            h64!("0x1.fe6e1742f7cf5p-1"),
307            h64!("0x1.cd5ec93c1243ap-7"),
308            h64!("-0x1.e2ff3aaacb386p-6"),
309            h64!("0x1.2aa4e5823cc89p-5"),
310            h64!("-0x1.d049842dbe399p-6"),
311            h64!("0x1.a34edb21ab302p-7"),
312            h64!("-0x1.676e5996c7f9bp-10"),
313            h64!("-0x1.23b01a35140bfp-9"),
314        ],
315        [
316            h64!("0x1.fed37386190fbp-1"),
317            h64!("0x1.61beae53b72c2p-7"),
318            h64!("-0x1.7d6193f22c3c1p-6"),
319            h64!("0x1.e947279e3bb7dp-6"),
320            h64!("-0x1.906031b97ca97p-6"),
321            h64!("0x1.8d14d62561755p-7"),
322            h64!("-0x1.1f245e7178882p-9"),
323            h64!("-0x1.9257d4eb47685p-10"),
324        ],
325        [
326            h64!("0x1.ff20e0a7ba8c2p-1"),
327            h64!("0x1.0d1d69569b839p-7"),
328            h64!("-0x1.2a8ca0dc02752p-6"),
329            h64!("0x1.8cc071b709751p-6"),
330            h64!("-0x1.54a149f1b070cp-6"),
331            h64!("0x1.6e9137b13412cp-7"),
332            h64!("-0x1.6577ed3d8e83bp-9"),
333            h64!("-0x1.e9c1a5178a289p-11"),
334        ],
335        [
336            h64!("0x1.ff5b8fb26f5f6p-1"),
337            h64!("0x1.9646f35a7663cp-8"),
338            h64!("-0x1.cf68ed9311b0bp-7"),
339            h64!("0x1.3e8735b5a694fp-6"),
340            h64!("-0x1.1e1612d026fdfp-6"),
341            h64!("0x1.4afd8e6ca636dp-7"),
342            h64!("-0x1.8c375170ccb22p-9"),
343            h64!("-0x1.c799443c4fd3bp-12"),
344        ],
345        [
346            h64!("0x1.ff87b1913e853p-1"),
347            h64!("0x1.30499b5039596p-8"),
348            h64!("-0x1.64964201ec8bap-7"),
349            h64!("0x1.fa73d7eafba98p-7"),
350            h64!("-0x1.daa3022141fbbp-7"),
351            h64!("0x1.2509444c063b7p-7"),
352            h64!("-0x1.99482a2f8a0a1p-9"),
353            h64!("-0x1.403d1f76c9454p-15"),
354        ],
355        [
356            h64!("0x1.ffa89fe5b3625p-1"),
357            h64!("0x1.c4412bf4b8f35p-9"),
358            h64!("-0x1.100f347126cf0p-7"),
359            h64!("0x1.8ebda07671d40p-7"),
360            h64!("-0x1.850c6a31c98c1p-7"),
361            h64!("0x1.fdac860c67d21p-8"),
362            h64!("-0x1.927d03d2ba12cp-9"),
363            h64!("0x1.0ff620b4190fep-12"),
364        ],
365        [
366            h64!("0x1.ffc10194fcb64p-1"),
367            h64!("0x1.4d78bba8ca621p-9"),
368            h64!("-0x1.9ba107a443e02p-8"),
369            h64!("0x1.36f273fbc04ccp-7"),
370            h64!("-0x1.3b38716ac7e6fp-7"),
371            h64!("0x1.b3fe0181914acp-8"),
372            h64!("-0x1.7d3fe7de98c5cp-9"),
373            h64!("0x1.ea31f8e5317f7p-12"),
374        ],
375        [
376            h64!("0x1.ffd2eae369a07p-1"),
377            h64!("0x1.e7f232d9e266cp-10"),
378            h64!("-0x1.34c7442dd48d9p-8"),
379            h64!("0x1.e066bed070a0bp-8"),
380            h64!("-0x1.f914f3c42fc0dp-8"),
381            h64!("0x1.6f4664ed2260fp-8"),
382            h64!("-0x1.5e59910761d24p-9"),
383            h64!("0x1.39cbb6e84c126p-11"),
384        ],
385        [
386            h64!("0x1.ffdff92db56e5p-1"),
387            h64!("0x1.6235fbd7a4373p-10"),
388            h64!("-0x1.cb5e029b9e56ap-9"),
389            h64!("0x1.6fa4c7ef274dap-8"),
390            h64!("-0x1.903a089a835f3p-8"),
391            h64!("0x1.30f12e0ca1901p-8"),
392            h64!("-0x1.39d21b6957f99p-9"),
393            h64!("0x1.5d3f8495a703cp-11"),
394        ],
395        [
396            h64!("0x1.ffe96a78a04a9p-1"),
397            h64!("0x1.fe41cd9bb4f2cp-11"),
398            h64!("-0x1.52d7b28966c0cp-9"),
399            h64!("0x1.16c192d86a1a7p-8"),
400            h64!("-0x1.39bfce951100cp-8"),
401            h64!("0x1.f376a7869f9e3p-9"),
402            h64!("-0x1.12e6cef999c4fp-9"),
403            h64!("0x1.66acd4d667b5p-11"),
404        ],
405        [
406            h64!("0x1.fff0312b010b5p-1"),
407            h64!("0x1.6caa0d3583018p-11"),
408            h64!("-0x1.efb729f4cf75bp-10"),
409            h64!("0x1.a2da7cebe12acp-9"),
410            h64!("-0x1.e6c27a24bc759p-9"),
411            h64!("0x1.93b1f4d8ea65p-9"),
412            h64!("-0x1.d82050aa94a08p-10"),
413            h64!("0x1.5cd7dc75d6cbap-11"),
414        ],
415        [
416            h64!("0x1.fff50456dab8cp-1"),
417            h64!("0x1.0295ef6591865p-11"),
418            h64!("-0x1.679880e95a4dap-10"),
419            h64!("0x1.37d38e3a5c8ebp-9"),
420            h64!("-0x1.75b3708aebb8fp-9"),
421            h64!("0x1.4231c4b4b0296p-9"),
422            h64!("-0x1.8e26476489318p-10"),
423            h64!("0x1.45c3b570dd924p-11"),
424        ],
425        [
426            h64!("0x1.fff86cfd3e657p-1"),
427            h64!("0x1.6be02102b353dp-12"),
428            h64!("-0x1.02b157780d6aep-10"),
429            h64!("0x1.cc1d886861133p-10"),
430            h64!("-0x1.1bff6f12ec9abp-9"),
431            h64!("0x1.fc0f77bd9c736p-10"),
432            h64!("-0x1.4a3320bd0959dp-10"),
433            h64!("0x1.267f8b4f95d2p-11"),
434        ],
435        [
436            h64!("0x1.fffad0b901755p-1"),
437            h64!("0x1.fc0d55470cf5ep-13"),
438            h64!("-0x1.7121aff5e820ep-11"),
439            h64!("0x1.506d6992f7de5p-10"),
440            h64!("-0x1.ab595d3ecd0d6p-10"),
441            h64!("0x1.8bdd79daaf754p-10"),
442            h64!("-0x1.0d9b090f997c1p-10"),
443            h64!("0x1.031ab9fd1c7dap-11"),
444        ],
445        [
446            h64!("0x1.fffc7a37857d2p-1"),
447            h64!("0x1.5feada379d8a5p-13"),
448            h64!("-0x1.05304df58f3aap-11"),
449            h64!("0x1.e79c081b8600fp-11"),
450            h64!("-0x1.3e5dbe33232e0p-10"),
451            h64!("0x1.30eb208200729p-10"),
452            h64!("-0x1.b1d493b147945p-11"),
453            h64!("0x1.bd587bbc071bep-12"),
454        ],
455        [
456            h64!("0x1.fffd9fdeabccep-1"),
457            h64!("0x1.e3bcf436a1a49p-14"),
458            h64!("-0x1.6e953111ef0a1p-12"),
459            h64!("0x1.5e3edf6768654p-11"),
460            h64!("-0x1.d5be67c0547a4p-11"),
461            h64!("0x1.d07d9ffa1d435p-11"),
462            h64!("-0x1.58328f5f358cap-11"),
463            h64!("0x1.76d42d95c42c4p-12"),
464        ],
465        [
466            h64!("0x1.fffe68f4fa777p-1"),
467            h64!("0x1.49e17724f4cddp-14"),
468            h64!("-0x1.fe48c44e229c1p-13"),
469            h64!("0x1.f2bd95d76f188p-12"),
470            h64!("-0x1.57388cb12d011p-11"),
471            h64!("0x1.5decc25c5c079p-11"),
472            h64!("-0x1.0d7499d1b0d2dp-11"),
473            h64!("0x1.359332c94ecdcp-12"),
474        ],
475        [
476            h64!("0x1.fffef1960d85dp-1"),
477            h64!("0x1.be6abbb10a4cdp-15"),
478            h64!("-0x1.6040381a8c313p-13"),
479            h64!("0x1.5fff1dde9ee9dp-12"),
480            h64!("-0x1.f0c933efa9971p-12"),
481            h64!("0x1.04cbf4a5cd760p-11"),
482            h64!("-0x1.a07f150af6dadp-12"),
483            h64!("0x1.f68dd183426bap-13"),
484        ],
485        [
486            h64!("0x1.ffff4db27f146p-1"),
487            h64!("0x1.2bb5cc22e5cd8p-15"),
488            h64!("-0x1.e25894899f526p-14"),
489            h64!("0x1.ec8a8e5a72757p-13"),
490            h64!("-0x1.64256ae0a3cf9p-12"),
491            h64!("0x1.80a836c18c46cp-12"),
492            h64!("-0x1.3dea401af6775p-12"),
493            h64!("0x1.915ddff3fe0d1p-13"),
494        ],
495        [
496            h64!("0x1.ffff8b500e77cp-1"),
497            h64!("0x1.8f4ccca7fc769p-16"),
498            h64!("-0x1.478cffe305946p-14"),
499            h64!("0x1.559f04adde504p-13"),
500            h64!("-0x1.f9e1577d6961dp-13"),
501            h64!("0x1.18bda53c14716p-12"),
502            h64!("-0x1.df8634c35541cp-13"),
503            h64!("0x1.3bb5c6b616337p-13"),
504        ],
505        [
506            h64!("0x1.ffffb43555b5fp-1"),
507            h64!("0x1.07ebd2a2d26c8p-16"),
508            h64!("-0x1.b93e442a37f2bp-15"),
509            h64!("0x1.d5cf15159ce28p-14"),
510            h64!("-0x1.63f5e1469c006p-13"),
511            h64!("0x1.95a03acebac18p-13"),
512            h64!("-0x1.656e5e2a1f8e2p-13"),
513            h64!("0x1.e98c437189bdep-14"),
514        ],
515        [
516            h64!("0x1.ffffcf23ff5fcp-1"),
517            h64!("0x1.5a2adfa0b492cp-17"),
518            h64!("-0x1.26c88270759f0p-15"),
519            h64!("0x1.40473572b99a8p-14"),
520            h64!("-0x1.f057cbde578a5p-14"),
521            h64!("0x1.22178d1c3c948p-13"),
522            h64!("-0x1.0765b61a0d859p-13"),
523            h64!("0x1.765b3ea03ddbep-14"),
524        ],
525        [
526            h64!("0x1.ffffe0bd3e852p-1"),
527            h64!("0x1.c282cd3957a72p-18"),
528            h64!("-0x1.86ad6dfa44faap-16"),
529            h64!("0x1.b0f313f03a029p-15"),
530            h64!("-0x1.56e44abecd255p-14"),
531            h64!("0x1.9ad1ecfe34a89p-14"),
532            h64!("-0x1.7fe4033478618p-14"),
533            h64!("0x1.1a8184e049fbfp-14"),
534        ],
535        [
536            h64!("0x1.ffffec2641a9ep-1"),
537            h64!("0x1.22df29821407ep-18"),
538            h64!("-0x1.00c902a6cfd98p-16"),
539            h64!("0x1.22234eb88671fp-15"),
540            h64!("-0x1.d57a181c9e6e1p-15"),
541            h64!("0x1.200c283b54a90p-14"),
542            h64!("-0x1.14b4c3295a7d0p-14"),
543            h64!("0x1.a4f966f713bdep-15"),
544        ],
545        [
546            h64!("0x1.fffff37d63a36p-1"),
547            h64!("0x1.74adc8f405eecp-19"),
548            h64!("-0x1.4ed4228e44858p-17"),
549            h64!("0x1.81918baea92bap-16"),
550            h64!("-0x1.3e81b17a0009cp-15"),
551            h64!("0x1.9004a36116436p-15"),
552            h64!("-0x1.8aa1ba400e076p-15"),
553            h64!("0x1.35cd4e2340a9ep-15"),
554        ],
555        [
556            h64!("0x1.fffff82cdcf1bp-1"),
557            h64!("0x1.d9c73698fa87dp-20"),
558            h64!("-0x1.b11017ec67115p-18"),
559            h64!("0x1.fc0dfadf653f8p-17"),
560            h64!("-0x1.ac4e03cd2dfc2p-16"),
561            h64!("0x1.131806b5abbc5p-15"),
562            h64!("-0x1.1672ef66fcaafp-15"),
563            h64!("0x1.c2882c7debed7p-16"),
564        ],
565        [
566            h64!("0x1.fffffb248c39dp-1"),
567            h64!("0x1.2acee2f5ec66ap-20"),
568            h64!("-0x1.15cc570408a36p-18"),
569            h64!("0x1.4be757bbb75a3p-17"),
570            h64!("-0x1.1d6aa5f8d2940p-16"),
571            h64!("0x1.76c5937d5105ep-16"),
572            h64!("-0x1.84dffc3ca9302p-16"),
573            h64!("0x1.43c8315f2c30ap-16"),
574        ],
575        [
576            h64!("0x1.fffffd01f36afp-1"),
577            h64!("0x1.75fa8dbc840bap-21"),
578            h64!("-0x1.6186da0133f5ap-19"),
579            h64!("0x1.ae023231e1af5p-18"),
580            h64!("-0x1.790812f7ca394p-17"),
581            h64!("0x1.f9c25656d0ef2p-17"),
582            h64!("-0x1.0cc66682e304cp-16"),
583            h64!("0x1.cc170a75d6f9cp-17"),
584        ],
585        [
586            h64!("0x1.fffffe2ba0ea5p-1"),
587            h64!("0x1.d06ad6ecde88ep-22"),
588            h64!("-0x1.be46aa8edc9a1p-20"),
589            h64!("0x1.143860c7840b8p-18"),
590            h64!("-0x1.edaba78fb1260p-18"),
591            h64!("0x1.52138a96ecee2p-17"),
592            h64!("-0x1.6fca538c4e2eep-17"),
593            h64!("0x1.434040640bcefp-17"),
594        ],
595        [
596            h64!("0x1.fffffee3cc32cp-1"),
597            h64!("0x1.1e1e857adb8ddp-22"),
598            h64!("-0x1.1769ce5f2a6e8p-20"),
599            h64!("0x1.5fe5d479b0543p-19"),
600            h64!("-0x1.405d865c94c2ap-18"),
601            h64!("0x1.bfc94feb96afcp-18"),
602            h64!("-0x1.f245d5f3e8358p-18"),
603            h64!("0x1.c142456acf443p-18"),
604        ],
605    ];
606    let ax = x.abs();
607    let t = ax;
608    let ux = t.to_bits();
609    let s = f64::from(x);
610    let mut z = f64::from(ax);
611    //for x > 0x1.f5a888p+1, erf(x) rounds to 1 to nearest
612    if unlikely(ux > 0x407ad444) {
613        // |x| > 0x1.f5a888p+1
614        let os = 1.0_f32.copysign(x);
615        if ux > (0xff << 23) {
616            return (x + x) as f16; // nan
617        }
618        if ux == (0xff << 23) {
619            return os as f16; // +-inf
620        }
621        return (os - h32!("0x1p-25") * os) as f16;
622    }
623    let v = (16.0 * z).floor();
624    let i = (16.0 * ax) as usize;
625    if unlikely(ux < 0x3ee00000) {
626        // |x| < 0x1.cp-2
627        const C: [f64; 8] = [
628            h64!("0x1.20dd750429b6dp+0"),
629            h64!("-0x1.812746b0375fbp-2"),
630            h64!("0x1.ce2f219fd6f45p-4"),
631            h64!("-0x1.b82ce2cbf0838p-6"),
632            h64!("0x1.565bb655adb85p-8"),
633            h64!("-0x1.c025bfc879c94p-11"),
634            h64!("0x1.f81718f61309cp-14"),
635            h64!("-0x1.cc67bd88f5867p-17"),
636        ];
637        let z2 = s * s;
638        let z4 = z2 * z2;
639        let z8 = z4 * z4;
640        let mut c0 = C[0] + z2 * C[1];
641        let c2 = C[2] + z2 * C[3];
642        let mut c4 = C[4] + z2 * C[5];
643        let c6 = C[6] + z2 * C[7];
644        c0 += z4 * c2;
645        c4 += z4 * c6;
646        c0 += z8 * c4;
647        return (s * c0) as f16;
648    }
649    z = (z - 0.03125) - 0.0625 * v;
650    let c = C[i - 7];
651    let z2 = z * z;
652    let z4 = z2 * z2;
653    let mut c0 = c[0] + z * c[1];
654    let c2 = c[2] + z * c[3];
655    let mut c4 = c[4] + z * c[5];
656    let c6 = c[6] + z * c[7];
657    c0 += z2 * c2;
658    c4 += z2 * c6;
659    c0 += z4 * c4;
660    c0.copysign(s) as f16
661}
662
663// TODO: replace with std::hint::unlikely
664const fn unlikely(b: bool) -> bool {
665    b
666}
667
668#[cfg(test)]
669mod tests {
670    #[test]
671    fn exhaustive() {
672        for b in 0..=u16::MAX {
673            let x = f16::from_bits(b);
674            let y1 = super::cr_erff16(x);
675            let y2 = core_math::erff16(x);
676            assert_eq!(
677                y1.to_bits(),
678                y2.to_bits(),
679                "erff16({x} @ {b:#04x}) = ({y1} @ {y1b:#04x}) vs ({y2} @ {y2b:#04x})",
680                y1b = y1.to_bits(),
681                y2b = y2.to_bits()
682            );
683        }
684    }
685
686    #[test]
687    fn edge() {
688        assert_eq!(super::cr_erff16(-f16::NAN).to_bits(), (-f16::NAN).to_bits());
689        assert_eq!(super::cr_erff16(-f16::INFINITY), -1.0);
690        assert_eq!(super::cr_erff16(-0.0).to_bits(), (-0.0_f16).to_bits());
691        assert_eq!(super::cr_erff16(0.0).to_bits(), (0.0_f16).to_bits());
692        assert_eq!(super::cr_erff16(f16::INFINITY), 1.0);
693        assert_eq!(super::cr_erff16(f16::NAN).to_bits(), (f16::NAN).to_bits());
694    }
695}
core_math_rs/f16/erf.rs

core_math_rs/f16/
erf.rs
