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, fetestexcept, FE_INEXACT, FE_TONEAREST, FE_UNDERFLOW, |
33 | }; |
34 | |
35 | /* |
36 | * Fused multiply-add: Compute x * y + z with a single rounding error. |
37 | * |
38 | * A double has more than twice as much precision than a float, so |
39 | * direct double-precision arithmetic suffices, except where double |
40 | * rounding occurs. |
41 | */ |
42 | |
43 | /// Floating multiply add (f32) |
44 | /// |
45 | /// Computes `(x*y)+z`, rounded as one ternary operation: |
46 | /// Computes the value (as if) to infinite precision and rounds once to the result format, |
47 | /// according to the rounding mode characterized by the value of FLT_ROUNDS. |
48 | #[cfg_attr (all(test, assert_no_panic), no_panic::no_panic)] |
49 | pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 { |
50 | let xy: f64; |
51 | let mut result: f64; |
52 | let mut ui: u64; |
53 | let e: i32; |
54 | |
55 | xy = x as f64 * y as f64; |
56 | result = xy + z as f64; |
57 | ui = result.to_bits(); |
58 | e = (ui >> 52) as i32 & 0x7ff; |
59 | /* Common case: The double precision result is fine. */ |
60 | if ( |
61 | /* not a halfway case */ |
62 | ui & 0x1fffffff) != 0x10000000 || |
63 | /* NaN */ |
64 | e == 0x7ff || |
65 | /* exact */ |
66 | (result - xy == z as f64 && result - z as f64 == xy) || |
67 | /* not round-to-nearest */ |
68 | fegetround() != FE_TONEAREST |
69 | { |
70 | /* |
71 | underflow may not be raised correctly, example: |
72 | fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f) |
73 | */ |
74 | if e < 0x3ff - 126 && e >= 0x3ff - 149 && fetestexcept(FE_INEXACT) != 0 { |
75 | feclearexcept(FE_INEXACT); |
76 | // prevent `xy + vz` from being CSE'd with `xy + z` above |
77 | let vz: f32 = unsafe { read_volatile(&z) }; |
78 | result = xy + vz as f64; |
79 | if fetestexcept(FE_INEXACT) != 0 { |
80 | feraiseexcept(FE_UNDERFLOW); |
81 | } else { |
82 | feraiseexcept(FE_INEXACT); |
83 | } |
84 | } |
85 | z = result as f32; |
86 | return z; |
87 | } |
88 | |
89 | /* |
90 | * If result is inexact, and exactly halfway between two float values, |
91 | * we need to adjust the low-order bit in the direction of the error. |
92 | */ |
93 | let neg = ui >> 63 != 0; |
94 | let err = if neg == (z as f64 > xy) { |
95 | xy - result + z as f64 |
96 | } else { |
97 | z as f64 - result + xy |
98 | }; |
99 | if neg == (err < 0.0) { |
100 | ui += 1; |
101 | } else { |
102 | ui -= 1; |
103 | } |
104 | f64::from_bits(ui) as f32 |
105 | } |
106 | |
107 | #[cfg (test)] |
108 | mod tests { |
109 | #[test ] |
110 | fn issue_263() { |
111 | let a = f32::from_bits(1266679807); |
112 | let b = f32::from_bits(1300234242); |
113 | let c = f32::from_bits(1115553792); |
114 | let expected = f32::from_bits(1501560833); |
115 | assert_eq!(super::fmaf(a, b, c), expected); |
116 | } |
117 | } |
118 | |