Skip to main content

c/
math_lrint.rs

1// Port of musl lrint/llrint/lrintf/llrintf/lrintl/llrintl
2//
3// On x86_64: c_long = i64 (LONG_MAX >= 2^53), c_longlong = i64, long double = f64.
4// All six functions reduce to: round-to-integer-using-current-FP-mode, then cast.
5//
6// The rint/rintf helpers use musl's `x + toint - toint` trick which rounds via
7// the current FP rounding mode (MXCSR on x86_64), unlike Rust's `as` which
8// truncates toward zero.
9
10// musl rint: rounds f64 to integer using current FP rounding mode
11#[inline]
12fn rint(x: f64) -> f64 {
13    let u = asuint64(x);
14    let e = ((u >> 52) & 0x7ff) as i32;
15    let s = u >> 63;
16
17    // Already an integer (exponent >= 52 means no fractional bits)
18    if e >= 0x3ff + 52 {
19        return x;
20    }
21
22    // toint = 1/DBL_EPSILON = 2^52 = 4503599627370496.0
23    let toint = 4503599627370496.0f64;
24    let y = if s != 0 {
25        x - toint + toint
26    } else {
27        x + toint - toint
28    };
29
30    if y == 0.0 {
31        if s != 0 { -0.0 } else { 0.0 }
32    } else {
33        y
34    }
35}
36
37// musl rintf: rounds f32 to integer using current FP rounding mode
38#[inline]
39fn rintf(x: f32) -> f32 {
40    let u = asuint(x);
41    let e = ((u >> 23) & 0xff) as i32;
42    let s = u >> 31;
43
44    // Already an integer (exponent >= 23 means no fractional bits)
45    if e >= 0x7f + 23 {
46        return x;
47    }
48
49    // toint = 1/FLT_EPSILON = 2^23 = 8388608.0
50    let toint = 8388608.0f32;
51    let y = if s != 0 {
52        x - toint + toint
53    } else {
54        x + toint - toint
55    };
56
57    if y == 0.0 {
58        if s != 0 { -0.0f32 } else { 0.0f32 }
59    } else {
60        y
61    }
62}
63
64// On x86_64 with c_long = i64 (LONG_MAX >= 2^53), the musl slow path
65// (lrint_slow with FE_INEXACT handling) is not needed. Simple cast suffices
66// because: (1) rint returns an exact integer as f64, (2) the largest integer
67// representable as f64 is 2^53 which fits in i64, (3) for |x| >= 2^52 rint
68// returns x unchanged and the saturating `as` cast handles overflow.
69
70#[no_mangle]
71pub extern "C" fn lrint(x: f64) -> c_long {
72    rint(x) as c_long
73}
74
75#[no_mangle]
76pub extern "C" fn llrint(x: f64) -> c_longlong {
77    rint(x) as c_longlong
78}
79
80#[no_mangle]
81pub extern "C" fn lrintf(x: f32) -> c_long {
82    rintf(x) as c_long
83}
84
85#[no_mangle]
86pub extern "C" fn llrintf(x: f32) -> c_longlong {
87    rintf(x) as c_longlong
88}
89
90// On x86_64 with -mlong-double-64: long double = f64, delegate to lrint/llrint.
91// On aarch64: long double = f128 (IEEE quad), cast to f64 then round.
92#[cfg(target_arch = "x86_64")]
93#[no_mangle]
94pub extern "C" fn lrintl(x: f64) -> c_long {
95    lrint(x)
96}
97
98#[cfg(target_arch = "x86_64")]
99#[no_mangle]
100pub extern "C" fn llrintl(x: f64) -> c_longlong {
101    llrint(x)
102}
103
104#[cfg(target_arch = "aarch64")]
105#[no_mangle]
106pub extern "C" fn lrintl(x: f128) -> c_long {
107    lrint(x as f64)
108}
109
110#[cfg(target_arch = "aarch64")]
111#[no_mangle]
112pub extern "C" fn llrintl(x: f128) -> c_longlong {
113    llrint(x as f64)
114}