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// long double = f64 on this platform (LDBL_MANT_DIG == 53)
91// musl: lrintl(x) -> lrint(x), llrintl(x) -> llrint(x)
92#[no_mangle]
93pub extern "C" fn lrintl(x: f64) -> c_long {
94    lrint(x)
95}
96
97#[no_mangle]
98pub extern "C" fn llrintl(x: f64) -> c_longlong {
99    llrint(x)
100}