Skip to main content

needle_core/
math.rs

1//! Thin wrappers around libm for no_std float math.
2
3#[inline(always)]
4pub fn exp(x: f32) -> f32 {
5    libm::expf(x)
6}
7
8#[inline(always)]
9pub fn ln(x: f32) -> f32 {
10    libm::logf(x)
11}
12
13#[inline(always)]
14pub fn sin(x: f32) -> f32 {
15    libm::sinf(x)
16}
17
18#[inline(always)]
19pub fn cos(x: f32) -> f32 {
20    libm::cosf(x)
21}
22
23#[inline(always)]
24pub fn sqrt(x: f32) -> f32 {
25    libm::sqrtf(x)
26}
27
28#[inline(always)]
29pub fn tanh(x: f32) -> f32 {
30    libm::tanhf(x)
31}
32
33#[inline(always)]
34pub fn powf(base: f32, exp: f32) -> f32 {
35    libm::powf(base, exp)
36}
37
38#[inline(always)]
39pub fn round(x: f32) -> f32 {
40    libm::roundf(x)
41}
42
43#[inline(always)]
44pub fn abs(x: f32) -> f32 {
45    libm::fabsf(x)
46}
47
48/// Convert IEEE 754 half-precision (binary16) bits to `f32`.
49///
50/// Exact for every input, including subnormals, infinities and NaN payloads.
51/// Lives here so the `no_std` v2 loader and `cq` can share it.
52#[inline]
53pub fn f16_to_f32(bits: u16) -> f32 {
54    let sign = ((bits >> 15) as u32) << 31;
55    let exp = ((bits >> 10) & 0x1F) as u32;
56    let mant = (bits & 0x3FF) as u32;
57
58    if exp == 0 {
59        // Subnormal (or zero): value = mant * 2^-24.
60        let val = mant as f32 / (1u32 << 24) as f32;
61        return if sign != 0 { -val } else { val };
62    }
63    if exp == 31 {
64        // Inf or NaN — preserve the mantissa payload.
65        return f32::from_bits(sign | 0x7F80_0000 | (mant << 13));
66    }
67    f32::from_bits(sign | ((exp + 127 - 15) << 23) | (mant << 13))
68}
69
70#[cfg(test)]
71mod f16_tests {
72    use super::f16_to_f32;
73
74    #[test]
75    fn known_values() {
76        assert_eq!(f16_to_f32(0x0000), 0.0);
77        assert_eq!(f16_to_f32(0x8000), -0.0);
78        assert_eq!(f16_to_f32(0x3C00), 1.0);
79        assert_eq!(f16_to_f32(0xBC00), -1.0);
80        assert_eq!(f16_to_f32(0x4000), 2.0);
81        assert_eq!(f16_to_f32(0x3555), 0.33325195); // nearest half to 1/3
82        assert_eq!(f16_to_f32(0x7BFF), 65504.0); // largest finite half
83        assert_eq!(f16_to_f32(0x0001), 5.9604645e-8); // smallest subnormal
84        assert_eq!(f16_to_f32(0x03FF), 6.0975552e-5); // largest subnormal
85        assert_eq!(f16_to_f32(0x0400), 6.1035156e-5); // smallest normal
86        assert!(f16_to_f32(0x7C00).is_infinite() && f16_to_f32(0x7C00) > 0.0);
87        assert!(f16_to_f32(0xFC00).is_infinite() && f16_to_f32(0xFC00) < 0.0);
88        assert!(f16_to_f32(0x7E00).is_nan());
89    }
90
91    /// Round-trip against the f32 -> f16 direction for every representable half.
92    #[test]
93    fn every_finite_half_round_trips() {
94        for bits in 0u32..=0xFFFF {
95            let b = bits as u16;
96            let v = f16_to_f32(b);
97            if v.is_nan() {
98                continue;
99            }
100            // f16 -> f32 is lossless, so re-narrowing must return the same bits.
101            let back = f32_to_f16_bits(v);
102            assert_eq!(back, b, "bits={b:#06x} value={v}");
103        }
104    }
105
106    /// Round-to-nearest-even narrowing, test-only reference.
107    fn f32_to_f16_bits(v: f32) -> u16 {
108        let x = v.to_bits();
109        let sign = ((x >> 31) as u16) << 15;
110        let exp = ((x >> 23) & 0xFF) as i32;
111        let mant = x & 0x7F_FFFF;
112        if exp == 0xFF {
113            return sign | 0x7C00 | ((mant >> 13) as u16);
114        }
115        let unbiased = exp - 127;
116        if unbiased < -24 {
117            return sign;
118        }
119        if unbiased < -14 {
120            // Subnormal half: shift the implicit 1 back in.
121            let shift = (-14 - unbiased) as u32;
122            let m = (mant | 0x80_0000) >> (shift + 13);
123            return sign | m as u16;
124        }
125        if unbiased > 15 {
126            return sign | 0x7C00;
127        }
128        sign | (((unbiased + 15) as u16) << 10) | ((mant >> 13) as u16)
129    }
130}