]>
Commit | Line | Data |
---|---|---|
8faf50e0 XL |
1 | /* origin: FreeBSD /usr/src/lib/msun/src/s_fmaf.c */ |
2 | /*- | |
3 | * Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG> | |
4 | * All rights reserved. | |
5 | * | |
6 | * Redistribution and use in source and binary forms, with or without | |
7 | * modification, are permitted provided that the following conditions | |
8 | * are met: | |
9 | * 1. Redistributions of source code must retain the above copyright | |
10 | * notice, this list of conditions and the following disclaimer. | |
11 | * 2. Redistributions in binary form must reproduce the above copyright | |
12 | * notice, this list of conditions and the following disclaimer in the | |
13 | * documentation and/or other materials provided with the distribution. | |
14 | * | |
15 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND | |
16 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE | |
17 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE | |
18 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE | |
19 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL | |
20 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS | |
21 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) | |
22 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT | |
23 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY | |
24 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF | |
25 | * SUCH DAMAGE. | |
26 | */ | |
27 | ||
28 | use core::f32; | |
29 | use core::ptr::read_volatile; | |
30 | ||
31 | use super::fenv::{ | |
32 | feclearexcept, fegetround, feraiseexcept, fesetround, fetestexcept, FE_INEXACT, FE_TONEAREST, | |
33 | FE_TOWARDZERO, FE_UNDERFLOW, | |
34 | }; | |
35 | ||
36 | /* | |
37 | * Fused multiply-add: Compute x * y + z with a single rounding error. | |
38 | * | |
39 | * A double has more than twice as much precision than a float, so | |
40 | * direct double-precision arithmetic suffices, except where double | |
41 | * rounding occurs. | |
42 | */ | |
dc9dc135 XL |
43 | |
44 | /// Floating multiply add (f32) | |
45 | /// | |
46 | /// Computes `(x*y)+z`, rounded as one ternary operation: | |
47 | /// Computes the value (as if) to infinite precision and rounds once to the result format, | |
48 | /// according to the rounding mode characterized by the value of FLT_ROUNDS. | |
48663c56 | 49 | #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] |
8faf50e0 XL |
50 | pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 { |
51 | let xy: f64; | |
52 | let mut result: f64; | |
53 | let mut ui: u64; | |
54 | let e: i32; | |
55 | ||
56 | xy = x as f64 * y as f64; | |
57 | result = xy + z as f64; | |
58 | ui = result.to_bits(); | |
59 | e = (ui >> 52) as i32 & 0x7ff; | |
60 | /* Common case: The double precision result is fine. */ | |
61 | if ( | |
62 | /* not a halfway case */ | |
63 | ui & 0x1fffffff) != 0x10000000 || | |
64 | /* NaN */ | |
65 | e == 0x7ff || | |
66 | /* exact */ | |
67 | (result - xy == z as f64 && result - z as f64 == xy) || | |
68 | /* not round-to-nearest */ | |
69 | fegetround() != FE_TONEAREST | |
70 | { | |
71 | /* | |
72 | underflow may not be raised correctly, example: | |
73 | fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f) | |
74 | */ | |
75 | if e < 0x3ff - 126 && e >= 0x3ff - 149 && fetestexcept(FE_INEXACT) != 0 { | |
76 | feclearexcept(FE_INEXACT); | |
77 | // prevent `xy + vz` from being CSE'd with `xy + z` above | |
78 | let vz: f32 = unsafe { read_volatile(&z) }; | |
79 | result = xy + vz as f64; | |
80 | if fetestexcept(FE_INEXACT) != 0 { | |
81 | feraiseexcept(FE_UNDERFLOW); | |
82 | } else { | |
83 | feraiseexcept(FE_INEXACT); | |
84 | } | |
85 | } | |
86 | z = result as f32; | |
87 | return z; | |
88 | } | |
89 | ||
90 | /* | |
91 | * If result is inexact, and exactly halfway between two float values, | |
92 | * we need to adjust the low-order bit in the direction of the error. | |
93 | */ | |
94 | fesetround(FE_TOWARDZERO); | |
95 | // prevent `vxy + z` from being CSE'd with `xy + z` above | |
96 | let vxy: f64 = unsafe { read_volatile(&xy) }; | |
97 | let mut adjusted_result: f64 = vxy + z as f64; | |
98 | fesetround(FE_TONEAREST); | |
99 | if result == adjusted_result { | |
100 | ui = adjusted_result.to_bits(); | |
101 | ui += 1; | |
102 | adjusted_result = f64::from_bits(ui); | |
103 | } | |
104 | z = adjusted_result as f32; | |
105 | z | |
106 | } |