pxfm/
round.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 6/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::bits::{get_exponent_f32, get_exponent_f64, mantissa_f32, mantissa_f64};
30
31#[inline]
32pub const fn roundf(x: f32) -> f32 {
33    // If x is infinity NaN or zero, return it.
34    if !x.is_normal() && !x.is_subnormal() {
35        return x;
36    }
37
38    let exponent = get_exponent_f32(x);
39
40    const FRACTION_LENGTH: u32 = 23;
41
42    // If the exponent is greater than the most negative mantissa
43    // exponent, then x is already an integer.
44    if exponent >= FRACTION_LENGTH as i32 {
45        return x;
46    }
47
48    if exponent == -1 {
49        // Absolute value of x is greater than equal to 0.5 but less than 1.
50        return if x.is_sign_negative() { -1.0 } else { 1.0 };
51    }
52
53    if exponent <= -2 {
54        // Absolute value of x is less than 0.5.
55        return if x.is_sign_negative() { -0.0 } else { 0.0 };
56    }
57
58    let trim_size = (FRACTION_LENGTH as i32).wrapping_sub(exponent);
59    let half_bit_set = mantissa_f32(x) & (1u32 << (trim_size - 1)) != 0;
60    let x_u = x.to_bits();
61    let trunc_u: u32 = (x_u >> trim_size).wrapping_shl(trim_size as u32);
62
63    // If x is already an integer, return it.
64    if trunc_u == x_u {
65        return x;
66    }
67
68    let trunc_value = f32::from_bits(trunc_u);
69
70    if !half_bit_set {
71        // Fractional part is less than 0.5 so round value is the
72        // same as the trunc value.
73        trunc_value
74    } else if x.is_sign_negative() {
75        trunc_value - 1.0
76    } else {
77        trunc_value + 1.0
78    }
79}
80
81// infinity, NaNs are assumed already handled somewhere
82#[inline]
83pub(crate) fn froundf_finite(x: f32) -> f32 {
84    #[cfg(any(
85        all(
86            any(target_arch = "x86", target_arch = "x86_64"),
87            target_feature = "sse4.1"
88        ),
89        target_arch = "aarch64"
90    ))]
91    {
92        x.round()
93    }
94    #[cfg(not(any(
95        all(
96            any(target_arch = "x86", target_arch = "x86_64"),
97            target_feature = "sse4.1"
98        ),
99        target_arch = "aarch64"
100    )))]
101    {
102        let exponent = get_exponent_f32(x);
103
104        const FRACTION_LENGTH: u32 = 23;
105
106        // If the exponent is greater than the most negative mantissa
107        // exponent, then x is already an integer.
108        if exponent >= FRACTION_LENGTH as i32 {
109            return x;
110        }
111
112        if exponent == -1 {
113            // Absolute value of x is greater than equal to 0.5 but less than 1.
114            return if x.is_sign_negative() { -1.0 } else { 1.0 };
115        }
116
117        if exponent <= -2 {
118            // Absolute value of x is less than 0.5.
119            return if x.is_sign_negative() { -0.0 } else { 0.0 };
120        }
121
122        let trim_size = (FRACTION_LENGTH as i32).wrapping_sub(exponent);
123        let half_bit_set = mantissa_f32(x) & (1u32 << (trim_size - 1)) != 0;
124        let x_u = x.to_bits();
125        let trunc_u: u32 = (x_u >> trim_size).wrapping_shl(trim_size as u32);
126
127        // If x is already an integer, return it.
128        if trunc_u == x_u {
129            return x;
130        }
131
132        let trunc_value = f32::from_bits(trunc_u);
133
134        if !half_bit_set {
135            // Fractional part is less than 0.5 so round value is the
136            // same as the trunc value.
137            trunc_value
138        } else if x.is_sign_negative() {
139            trunc_value - 1.0
140        } else {
141            trunc_value + 1.0
142        }
143    }
144}
145
146#[inline]
147pub const fn round(x: f64) -> f64 {
148    // If x is infinity NaN or zero, return it.
149    if !x.is_normal() && !x.is_subnormal() {
150        return x;
151    }
152
153    let exponent = get_exponent_f64(x);
154
155    const FRACTION_LENGTH: u64 = 52;
156
157    // If the exponent is greater than the most negative mantissa
158    // exponent, then x is already an integer.
159    if exponent >= FRACTION_LENGTH as i64 {
160        return x;
161    }
162
163    if exponent == -1 {
164        // Absolute value of x is greater than equal to 0.5 but less than 1.
165        return if x.is_sign_negative() { -1.0 } else { 1.0 };
166    }
167
168    if exponent <= -2 {
169        // Absolute value of x is less than 0.5.
170        return if x.is_sign_negative() { -0.0 } else { 0.0 };
171    }
172
173    let trim_size = (FRACTION_LENGTH as i64).wrapping_sub(exponent);
174    let half_bit_set = mantissa_f64(x) & (1u64 << (trim_size.wrapping_sub(1))) != 0;
175    let x_u = x.to_bits();
176    let trunc_u: u64 = (x_u >> trim_size).wrapping_shl(trim_size as u32);
177
178    // If x is already an integer, return it.
179    if trunc_u == x_u {
180        return x;
181    }
182
183    let trunc_value = f64::from_bits(trunc_u);
184
185    if !half_bit_set {
186        // Fractional part is less than 0.5 so round value is the
187        // same as the trunc value.
188        trunc_value
189    } else if x.is_sign_negative() {
190        trunc_value - 1.0
191    } else {
192        trunc_value + 1.0
193    }
194}
195
196// infinity, NaNs are assumed already handled somewhere
197#[inline]
198pub(crate) fn fround_finite(x: f64) -> f64 {
199    #[cfg(any(
200        all(
201            any(target_arch = "x86", target_arch = "x86_64"),
202            target_feature = "sse4.1"
203        ),
204        target_arch = "aarch64"
205    ))]
206    {
207        x.round()
208    }
209    #[cfg(not(any(
210        all(
211            any(target_arch = "x86", target_arch = "x86_64"),
212            target_feature = "sse4.1"
213        ),
214        target_arch = "aarch64"
215    )))]
216    {
217        let exponent = get_exponent_f64(x);
218
219        const FRACTION_LENGTH: u64 = 52;
220
221        // If the exponent is greater than the most negative mantissa
222        // exponent, then x is already an integer.
223        if exponent >= FRACTION_LENGTH as i64 {
224            return x;
225        }
226
227        if exponent == -1 {
228            // Absolute value of x is greater than equal to 0.5 but less than 1.
229            return if x.is_sign_negative() { -1.0 } else { 1.0 };
230        }
231
232        if exponent <= -2 {
233            // Absolute value of x is less than 0.5.
234            return if x.is_sign_negative() { -0.0 } else { 0.0 };
235        }
236
237        let trim_size = (FRACTION_LENGTH as i64).wrapping_sub(exponent);
238        let half_bit_set = mantissa_f64(x) & (1u64 << (trim_size.wrapping_sub(1))) != 0;
239        let x_u = x.to_bits();
240        let trunc_u: u64 = (x_u >> trim_size).wrapping_shl(trim_size as u32);
241
242        // If x is already an integer, return it.
243        if trunc_u == x_u {
244            return x;
245        }
246
247        let trunc_value = f64::from_bits(trunc_u);
248
249        if !half_bit_set {
250            // Fractional part is less than 0.5 so round value is the
251            // same as the trunc value.
252            trunc_value
253        } else if x.is_sign_negative() {
254            trunc_value - 1.0
255        } else {
256            trunc_value + 1.0
257        }
258    }
259}
260
261pub(crate) trait RoundFinite {
262    fn round_finite(self) -> Self;
263}
264
265impl RoundFinite for f32 {
266    #[inline]
267    fn round_finite(self) -> Self {
268        froundf_finite(self)
269    }
270}
271
272impl RoundFinite for f64 {
273    #[inline]
274    fn round_finite(self) -> Self {
275        fround_finite(self)
276    }
277}
278
279#[cfg(test)]
280mod tests {
281    use super::*;
282
283    #[test]
284    fn test_roundf() {
285        assert_eq!(roundf(0f32), 0.0f32.round());
286        assert_eq!(roundf(1f32), 1.0f32.round());
287        assert_eq!(roundf(1.2f32), 1.2f32.round());
288        assert_eq!(roundf(-1.2f32), (-1.2f32).round());
289        assert_eq!(roundf(-1.6f32), (-1.6f32).round());
290        assert_eq!(roundf(-1.5f32), (-1.5f32).round());
291        assert_eq!(roundf(1.6f32), 1.6f32.round());
292        assert_eq!(roundf(1.5f32), 1.5f32.round());
293        assert_eq!(roundf(2.5f32), 2.5f32.round());
294    }
295
296    #[test]
297    fn test_round() {
298        assert_eq!(round(0.), 0.0f64.round());
299        assert_eq!(round(1.), 1.0f64.round());
300        assert_eq!(round(1.2), 1.2f64.round());
301        assert_eq!(round(-1.2), (-1.2f64).round());
302        assert_eq!(round(-1.6), (-1.6f64).round());
303        assert_eq!(round(-1.5), (-1.5f64).round());
304        assert_eq!(round(1.6), 1.6f64.round());
305        assert_eq!(round(1.5), 1.5f64.round());
306        assert_eq!(round(2.5), 2.5f64.round());
307    }
308}