zenith_float_num/ieee_soft/
mod.rs

1//! Software IEEE-754 binary32 / binary64 (`Ieee32` / `Ieee64`).
2//! Bits live in `u32` / `u64`. Arithmetic is integer only.
3
4mod arith;
5mod array;
6mod convert;
7mod simd;
8
9use crate::defs::Error;
10use crate::RoundingMode;
11use crate::Sign;
12use arith::{
13    add_bits, cmp_bits, div_bits, fma_bits, frexp_bits, from_i32_bits, mul_bits, next_after,
14    next_down, next_up, sqrt_bits, sub_bits, unpack, Class, BIN32, BIN64,
15};
16use core::cmp::Ordering;
17use core::num::FpCategory;
18
19pub use array::{ExactNumArray, Ieee32Array, Ieee64Array};
20pub use simd::IEEE_SIMD_LANE_WIDTH;
21
22/// Software IEEE-754 binary32 (24-bit significand, 8-bit exponent). Stored as `u32` bits.
23#[derive(Clone, Copy, Debug)]
24pub struct Ieee32(u32);
25
26/// Software IEEE-754 binary64 (53-bit significand, 11-bit exponent). Stored as `u64` bits.
27#[derive(Clone, Copy, Debug)]
28pub struct Ieee64(u64);
29
30macro_rules! impl_ieee {
31    ($ty:ident, $bits:ty, $fmt:ident, $as_u64:ident, $from_u64:ident) => {
32        impl $ty {
33            /// \(+0\).
34            pub const ZERO: Self = Self(0);
35
36            /// Build from the IEEE bit pattern.
37            #[inline]
38            pub const fn from_bits(bits: $bits) -> Self {
39                Self(bits)
40            }
41
42            /// Return the IEEE bit pattern.
43            #[inline]
44            pub const fn to_bits(self) -> $bits {
45                self.0
46            }
47
48            fn raw(self) -> u64 {
49                $as_u64(self.0)
50            }
51
52            /// Integer \(n\) converted with IEEE rounding.
53            pub fn from_i32(n: i32) -> Self {
54                Self($from_u64(from_i32_bits(n, $fmt)))
55            }
56
57            /// Classify the value.
58            pub fn classify(self) -> FpCategory {
59                match unpack(self.raw(), $fmt).class {
60                    Class::Nan => FpCategory::Nan,
61                    Class::Inf => FpCategory::Infinite,
62                    Class::Zero => FpCategory::Zero,
63                    Class::Sub => FpCategory::Subnormal,
64                    Class::Norm => FpCategory::Normal,
65                }
66            }
67
68            /// True if the value is NaN.
69            pub fn is_nan(self) -> bool {
70                unpack(self.raw(), $fmt).class == Class::Nan
71            }
72
73            /// True if the value is \(\pm\infty\).
74            pub fn is_infinite(self) -> bool {
75                unpack(self.raw(), $fmt).class == Class::Inf
76            }
77
78            /// True if the value is finite (including zero and subnormals).
79            pub fn is_finite(self) -> bool {
80                !self.is_nan() && !self.is_infinite()
81            }
82
83            /// True if the value is \(\pm 0\).
84            pub fn is_zero(self) -> bool {
85                unpack(self.raw(), $fmt).class == Class::Zero
86            }
87
88            /// True if the value is subnormal.
89            pub fn is_subnormal(self) -> bool {
90                unpack(self.raw(), $fmt).class == Class::Sub
91            }
92
93            /// True if the sign bit is set.
94            pub fn is_sign_negative(self) -> bool {
95                unpack(self.raw(), $fmt).sign
96            }
97
98            /// Associated error on NaN (`InvalidArgument`); `None` otherwise.
99            pub fn err(self) -> Option<Error> {
100                if self.is_nan() {
101                    Some(Error::InvalidArgument)
102                } else {
103                    None
104                }
105            }
106
107            /// IEEE add (to-nearest, ties to even).
108            pub fn add(self, rhs: Self) -> Self {
109                Self($from_u64(add_bits(self.raw(), rhs.raw(), $fmt)))
110            }
111
112            /// IEEE sub.
113            pub fn sub(self, rhs: Self) -> Self {
114                Self($from_u64(sub_bits(self.raw(), rhs.raw(), $fmt)))
115            }
116
117            /// IEEE mul.
118            pub fn mul(self, rhs: Self) -> Self {
119                Self($from_u64(mul_bits(self.raw(), rhs.raw(), $fmt)))
120            }
121
122            /// IEEE div. \(x/0\) is \(\pm\infty\); \(0/0\) is NaN.
123            pub fn div(self, rhs: Self) -> Self {
124                Self($from_u64(div_bits(self.raw(), rhs.raw(), $fmt)))
125            }
126
127            /// IEEE sqrt. Negative finite → NaN.
128            pub fn sqrt(self) -> Self {
129                Self($from_u64(sqrt_bits(self.raw(), $fmt)))
130            }
131
132            /// IEEE fused multiply-add \(a\cdot b + c\).
133            pub fn mul_add(self, b: Self, c: Self) -> Self {
134                Self($from_u64(fma_bits(self.raw(), b.raw(), c.raw(), $fmt)))
135            }
136
137            /// Next representable value toward \(+\infty\).
138            pub fn next_up(self) -> Self {
139                Self($from_u64(next_up(self.raw(), $fmt)))
140            }
141
142            /// Next representable value toward \(-\infty\).
143            pub fn next_down(self) -> Self {
144                Self($from_u64(next_down(self.raw(), $fmt)))
145            }
146
147            /// Next representable value toward `other`.
148            pub fn next_after(self, other: Self) -> Self {
149                Self($from_u64(next_after(self.raw(), other.raw(), $fmt)))
150            }
151
152            /// Split as \(m\cdot 2^e\) with \(m\in[1/2,1)\) (or a special).
153            pub fn frexp(self) -> (Self, i32) {
154                let (m, e) = frexp_bits(self.raw(), $fmt);
155                (Self($from_u64(m)), e)
156            }
157
158            /// Negate (flip the sign bit).
159            pub fn neg(self) -> Self {
160                Self($from_u64(self.raw() ^ $fmt.sign_mask()))
161            }
162
163            /// Absolute value.
164            pub fn abs(self) -> Self {
165                Self($from_u64(self.raw() & !$fmt.sign_mask()))
166            }
167
168            /// Copy the sign of `sign` onto `self`.
169            pub fn copysign(self, sign: Self) -> Self {
170                let mag = self.raw() & !$fmt.sign_mask();
171                let s = sign.raw() & $fmt.sign_mask();
172                Self($from_u64(mag | s))
173            }
174
175            /// Sign, or `None` if NaN.
176            pub fn sign(self) -> Option<Sign> {
177                if self.is_nan() {
178                    None
179                } else if unpack(self.raw(), $fmt).sign {
180                    Some(Sign::Neg)
181                } else {
182                    Some(Sign::Pos)
183                }
184            }
185
186            /// Compare; `None` if either is NaN. \(+0 = -0\).
187            pub fn cmp(self, other: Self) -> Option<Ordering> {
188                cmp_bits(self.raw(), other.raw(), $fmt).map(|d| {
189                    if d < 0 {
190                        Ordering::Less
191                    } else if d > 0 {
192                        Ordering::Greater
193                    } else {
194                        Ordering::Equal
195                    }
196                })
197            }
198        }
199
200        impl PartialEq for $ty {
201            fn eq(&self, other: &Self) -> bool {
202                self.cmp(*other) == Some(Ordering::Equal)
203            }
204        }
205
206        impl core::ops::Add for $ty {
207            type Output = Self;
208            fn add(self, rhs: Self) -> Self {
209                $ty::add(self, rhs)
210            }
211        }
212        impl core::ops::Sub for $ty {
213            type Output = Self;
214            fn sub(self, rhs: Self) -> Self {
215                $ty::sub(self, rhs)
216            }
217        }
218        impl core::ops::Mul for $ty {
219            type Output = Self;
220            fn mul(self, rhs: Self) -> Self {
221                $ty::mul(self, rhs)
222            }
223        }
224        impl core::ops::Div for $ty {
225            type Output = Self;
226            fn div(self, rhs: Self) -> Self {
227                $ty::div(self, rhs)
228            }
229        }
230        impl core::ops::Neg for $ty {
231            type Output = Self;
232            fn neg(self) -> Self {
233                $ty::neg(self)
234            }
235        }
236    };
237}
238
239const fn u32_as_u64(x: u32) -> u64 {
240    x as u64
241}
242const fn u64_as_u64(x: u64) -> u64 {
243    x
244}
245const fn u64_as_u32(x: u64) -> u32 {
246    x as u32
247}
248const fn u64_id(x: u64) -> u64 {
249    x
250}
251
252impl_ieee!(Ieee32, u32, BIN32, u32_as_u64, u64_as_u32);
253impl_ieee!(Ieee64, u64, BIN64, u64_as_u64, u64_id);
254
255impl Ieee32 {
256    /// Canonical quiet NaN.
257    pub const NAN: Self = Self(0x7FC0_0000);
258    /// \(+\infty\).
259    pub const INFINITY: Self = Self(0x7F80_0000);
260    /// \(-\infty\).
261    pub const NEG_INFINITY: Self = Self(0xFF80_0000);
262    /// \(-0\).
263    pub const NEG_ZERO: Self = Self(0x8000_0000);
264}
265
266impl Ieee64 {
267    /// Canonical quiet NaN.
268    pub const NAN: Self = Self(0x7FF8_0000_0000_0000);
269    /// \(+\infty\).
270    pub const INFINITY: Self = Self(0x7FF0_0000_0000_0000);
271    /// \(-\infty\).
272    pub const NEG_INFINITY: Self = Self(0xFFF0_0000_0000_0000);
273    /// \(-0\).
274    pub const NEG_ZERO: Self = Self(0x8000_0000_0000_0000);
275}
276
277impl Ieee32 {
278    /// Widen to `ExactNum` at precision `p` (use at least 32 bits).
279    pub fn to_exact(self, p: usize) -> crate::ExactNum {
280        convert::to_exact(self.to_bits() as u64, BIN32, p)
281    }
282
283    /// Round an `ExactNum` to binary32 (to-nearest, ties to even).
284    pub fn from_exact(x: &crate::ExactNum) -> Self {
285        Self::from_bits(convert::from_exact(x, BIN32) as u32)
286    }
287
288    /// Widen then apply `op`, then IEEE-round back.
289    pub fn map_exact<F>(self, p: usize, _rm: RoundingMode, op: F) -> Self
290    where
291        F: FnOnce(&crate::ExactNum) -> crate::ExactNum,
292    {
293        Self::from_exact(&op(&self.to_exact(p)))
294    }
295}
296
297impl Ieee64 {
298    /// Widen to `ExactNum` at precision `p` (use at least 64 bits).
299    pub fn to_exact(self, p: usize) -> crate::ExactNum {
300        convert::to_exact(self.to_bits(), BIN64, p)
301    }
302
303    /// Round an `ExactNum` to binary64 (to-nearest, ties to even).
304    pub fn from_exact(x: &crate::ExactNum) -> Self {
305        Self::from_bits(convert::from_exact(x, BIN64))
306    }
307
308    /// Widen then apply `op`, then IEEE-round back.
309    pub fn map_exact<F>(self, p: usize, _rm: RoundingMode, op: F) -> Self
310    where
311        F: FnOnce(&crate::ExactNum) -> crate::ExactNum,
312    {
313        Self::from_exact(&op(&self.to_exact(p)))
314    }
315}
316
317#[cfg(test)]
318mod tests {
319    use super::*;
320
321    const ONE32: u32 = 0x3F80_0000;
322    const TWO32: u32 = 0x4000_0000;
323    const HALF32: u32 = 0x3F00_0000;
324    const THIRD32: u32 = 0x3EAA_AAAB;
325    const THREE32: u32 = 0x4040_0000;
326    const INF32: u32 = 0x7F80_0000;
327    const QNAN32: u32 = 0x7FC0_0000;
328
329    const ONE64: u64 = 0x3FF0_0000_0000_0000;
330    const TWO64: u64 = 0x4000_0000_0000_0000;
331    const HALF64: u64 = 0x3FE0_0000_0000_0000;
332    const THIRD64: u64 = 0x3FD5_5555_5555_5555;
333    const THREE64: u64 = 0x4008_0000_0000_0000;
334    const INF64: u64 = 0x7FF0_0000_0000_0000;
335    const QNAN64: u64 = 0x7FF8_0000_0000_0000;
336    const MONE64: u64 = 0xBFF0_0000_0000_0000;
337
338    #[test]
339    fn bin32_add_mul_div_bits() {
340        let one = Ieee32::from_bits(ONE32);
341        let two = Ieee32::from_bits(TWO32);
342        let half = Ieee32::from_bits(HALF32);
343        let three = Ieee32::from_bits(THREE32);
344        assert_eq!((one + one).to_bits(), TWO32);
345        assert_eq!((one + two).to_bits(), THREE32);
346        assert_eq!((two * half).to_bits(), ONE32);
347        assert_eq!((one / Ieee32::from_i32(3)).to_bits(), THIRD32);
348        assert_eq!((three - one).to_bits(), TWO32);
349        assert!(Ieee32::from_bits(0)
350            .add(Ieee32::from_bits(0x8000_0000))
351            .is_zero());
352        assert_eq!(Ieee32::from_i32(1).to_bits(), ONE32);
353        assert_eq!(Ieee32::from_i32(-1).to_bits(), 0xBF80_0000);
354        assert_eq!(one.sqrt().to_bits(), ONE32);
355        assert!((Ieee32::from_i32(-1).sqrt()).is_nan());
356        assert_eq!((one / Ieee32::ZERO).to_bits(), INF32);
357        assert!((Ieee32::ZERO / Ieee32::ZERO).is_nan());
358        assert_eq!(Ieee32::INFINITY.to_bits(), INF32);
359        assert_eq!(Ieee32::NAN.to_bits(), QNAN32);
360        let tiny = Ieee32::from_bits(1);
361        assert!(tiny.is_subnormal());
362        assert_eq!((tiny + tiny).to_bits(), 2);
363    }
364
365    #[test]
366    fn bin64_add_mul_div_bits() {
367        let one = Ieee64::from_bits(ONE64);
368        let two = Ieee64::from_bits(TWO64);
369        let half = Ieee64::from_bits(HALF64);
370        let three = Ieee64::from_bits(THREE64);
371        assert_eq!((one + one).to_bits(), TWO64);
372        assert_eq!((one + two).to_bits(), THREE64);
373        assert_eq!((two * half).to_bits(), ONE64);
374        assert_eq!((one / Ieee64::from_i32(3)).to_bits(), THIRD64);
375        assert_eq!((three - one).to_bits(), TWO64);
376        assert_eq!(Ieee64::from_i32(1).to_bits(), ONE64);
377        assert_eq!(Ieee64::from_i32(-1).to_bits(), MONE64);
378        assert_eq!(one.sqrt().to_bits(), ONE64);
379        assert_eq!(Ieee64::from_i32(4).sqrt().to_bits(), TWO64);
380        // IEEE sqrt(2) (to-nearest, ties to even).
381        assert_eq!(two.sqrt().to_bits(), 0x3FF6_A09E_667F_3BCD);
382        assert_eq!((one / Ieee64::ZERO).to_bits(), INF64);
383        assert_eq!(Ieee64::INFINITY.to_bits(), INF64);
384        assert_eq!(Ieee64::NAN.to_bits(), QNAN64);
385        let tiny = Ieee64::from_bits(1);
386        assert!(tiny.is_subnormal());
387        assert_eq!((tiny + tiny).to_bits(), 2);
388        assert_eq!(one.mul_add(two, one).to_bits(), THREE64);
389        assert_eq!(one.neg().to_bits(), MONE64);
390        assert_eq!(Ieee64::ZERO.cmp(Ieee64::NEG_ZERO), Some(Ordering::Equal));
391        assert_eq!(one.next_up().next_down().to_bits(), ONE64);
392        let (m, e) = two.frexp();
393        assert_eq!(e, 2);
394        assert_eq!(m.to_bits(), HALF64);
395    }
396
397    #[test]
398    fn widen_roundtrip_one() {
399        let x = Ieee64::from_bits(ONE64);
400        let e = x.to_exact(64);
401        assert_eq!(Ieee64::from_exact(&e).to_bits(), ONE64);
402        let y = Ieee32::from_bits(ONE32);
403        assert_eq!(Ieee32::from_exact(&y.to_exact(64)).to_bits(), ONE32);
404    }
405}
zenith_float_num/ieee_soft/mod.rs

zenith_float_num/ieee_soft/
mod.rs
