nanvm_lib/vm/number/
binary64.rs

1//! The IEEE 754 binary64 layout of a [`Number`], in one place: a sign bit,
2//! an 11-bit biased exponent and a 52-bit fraction. [`Number::decompose`]
3//! reads a finite value as `±significand · 2^exponent`, and
4//! [`Number::compose`] rounds such a value, of any width, back to the
5//! nearest `Number`.
6
7use super::Number;
8use crate::sign::Sign;
9
10/// The width of the fraction field: the significand bits stored under the
11/// implicit leading one.
12const FRACTION_BITS: u32 = 52;
13
14/// The fraction field.
15const FRACTION_MASK: u64 = 0xf_ffff_ffff_ffff;
16
17/// The biased exponent field, once the fraction is shifted away. All ones
18/// is an infinity or a `NaN`.
19const EXPONENT_MASK: u64 = 0x7ff;
20
21/// The sign bit.
22const SIGN_BIT: u64 = 0x8000_0000_0000_0000;
23
24/// The leading bit of a normal value's significand, which the format
25/// implies rather than stores.
26const IMPLICIT_BIT: u64 = 0x10_0000_0000_0000;
27
28/// The number of significand bits a `Number` holds: the fraction and the
29/// implicit bit.
30const PRECISION: u32 = FRACTION_BITS + 1;
31
32/// What the biased exponent field exceeds the exponent of the significand's
33/// lowest bit by: the format's bias of 1023 plus the 52 fraction bits, as
34/// the significand is read as an integer.
35const BIAS: i32 = 1_075;
36
37/// The exponent of the significand's lowest bit for a biased exponent of
38/// zero or one: a subnormal's, and the smallest normal's. `2^-1074` is the
39/// smallest positive `Number`.
40const MIN_EXPONENT: i32 = 1 - BIAS;
41
42impl Number {
43    /// `self` as `±significand · 2^exponent`, exactly: a normal value's
44    /// `significand` has the implicit leading bit folded in, a subnormal's
45    /// is the fraction as it is, and both zeros are a `0` significand with
46    /// their own sign. `None` for a `NaN` and the infinities.
47    ///
48    /// Not reduced to lowest terms: `significand` may end in zero bits.
49    pub fn decompose(self) -> Option<(Sign, u64, i32)> {
50        let bits = self.0.to_bits();
51        let sign = if bits & SIGN_BIT == 0 {
52            Sign::Positive
53        } else {
54            Sign::Negative
55        };
56        let biased = (bits >> FRACTION_BITS) & EXPONENT_MASK;
57        let fraction = bits & FRACTION_MASK;
58        match biased {
59            EXPONENT_MASK => None,
60            0 => Some((sign, fraction, MIN_EXPONENT)),
61            _ => Some((sign, fraction | IMPLICIT_BIT, biased as i32 - BIAS)),
62        }
63    }
64
65    /// `±significand · 2^exponent` rounded to the nearest `Number`, a tie
66    /// going to the even significand: subnormal under the normal range, an
67    /// infinity past it, and a zero of `sign`'s sign for a `0` significand.
68    ///
69    /// `sticky` says the exact value is a little more than that: some bit
70    /// under `significand`'s lowest is set. It can only break a tie upward,
71    /// so it needs a rounded-away bit of `significand` above it to say where
72    /// the half is. With 53 significant bits or fewer the extra could be as
73    /// much as a whole unit of the result, so that call is refused with
74    /// `None`, in every build: the types can't carry the precondition.
75    pub(crate) fn compose(
76        sign: Sign,
77        significand: u128,
78        exponent: i32,
79        sticky: bool,
80    ) -> Option<Self> {
81        if sticky && significand >> PRECISION == 0 {
82            return None;
83        }
84        let sign_bit = match sign {
85            Sign::Positive => 0,
86            Sign::Negative => SIGN_BIT,
87        };
88        if significand == 0 {
89            return Some(Number(f64::from_bits(sign_bit)));
90        }
91        let width = 128 - significand.leading_zeros();
92        // The exponent of the result's lowest significand bit: `PRECISION`
93        // bits under the top of `significand`, but never under a
94        // subnormal's.
95        let unit = (i64::from(exponent) + i64::from(width) - i64::from(PRECISION))
96            .max(i64::from(MIN_EXPONENT));
97        // How many of `significand`'s low bits the result has no room for.
98        let dropped = unit - i64::from(exponent);
99        let kept = if dropped <= 0 {
100            // Exact: `significand` fits, at most `PRECISION` bits once
101            // shifted up to `unit`.
102            (significand << dropped.unsigned_abs()) as u64
103        } else {
104            let dropped = u32::try_from(dropped).unwrap_or(u32::MAX);
105            let kept = significand.checked_shr(dropped).unwrap_or(0) as u64;
106            let remainder = if dropped < 128 {
107                significand & ((1 << dropped) - 1)
108            } else {
109                significand
110            };
111            // `half` is the weight of the highest dropped bit; past 128 bits
112            // it exceeds every `remainder`.
113            let round_up = 1u128.checked_shl(dropped - 1).is_some_and(|half| {
114                remainder > half || (remainder == half && (sticky || kept & 1 != 0))
115            });
116            kept + u64::from(round_up)
117        };
118        // Rounding up all ones carries out of the significand: one bit
119        // fewer, one exponent higher, the same value.
120        let (kept, unit) = if kept >> PRECISION != 0 {
121            (kept >> 1, unit + 1)
122        } else {
123            (kept, unit)
124        };
125        let magnitude = if kept < IMPLICIT_BIT {
126            // Subnormal, or rounded down to zero: `unit` is `MIN_EXPONENT`,
127            // and the biased exponent is `0`.
128            kept
129        } else {
130            let biased = unit + i64::from(BIAS);
131            if biased >= EXPONENT_MASK as i64 {
132                EXPONENT_MASK << FRACTION_BITS
133            } else {
134                ((biased as u64) << FRACTION_BITS) | (kept & FRACTION_MASK)
135            }
136        };
137        Some(Number(f64::from_bits(sign_bit | magnitude)))
138    }
139}
140
141#[cfg(test)]
142mod tests {
143    use crate::{sign::Sign, vm::Number};
144
145    fn decompose(v: f64) -> Option<(Sign, u64, i32)> {
146        Number::from(v).decompose()
147    }
148
149    fn compose(significand: u128, exponent: i32, sticky: bool) -> f64 {
150        Number::compose(Sign::Positive, significand, exponent, sticky)
151            .unwrap()
152            .into()
153    }
154
155    /// Every value `decompose` reads, `compose` rebuilds bit for bit.
156    fn round_trip(v: f64) {
157        let (sign, significand, exponent) = decompose(v).unwrap();
158        let back =
159            f64::from(Number::compose(sign, u128::from(significand), exponent, false).unwrap());
160        assert_eq!(back.to_bits(), v.to_bits(), "{v:e}");
161    }
162
163    #[test]
164    fn decompose_normal() {
165        assert_eq!(
166            decompose(1.0),
167            Some((Sign::Positive, 0x10_0000_0000_0000, -52))
168        );
169        assert_eq!(
170            decompose(-1.5),
171            Some((Sign::Negative, 0x18_0000_0000_0000, -52))
172        );
173        // 2^53: the first integer whose successor is not a `Number`.
174        assert_eq!(
175            decompose(9_007_199_254_740_992.0),
176            Some((Sign::Positive, 0x10_0000_0000_0000, 1))
177        );
178        assert_eq!(
179            decompose(f64::MAX),
180            Some((Sign::Positive, 0x1f_ffff_ffff_ffff, 971))
181        );
182    }
183
184    #[test]
185    fn decompose_subnormal_boundary() {
186        // The smallest normal: the implicit bit, at the subnormals' exponent.
187        assert_eq!(
188            decompose(f64::MIN_POSITIVE),
189            Some((Sign::Positive, 0x10_0000_0000_0000, -1_074))
190        );
191        // The largest subnormal: no implicit bit, the same exponent.
192        assert_eq!(
193            decompose(f64::from_bits(0xf_ffff_ffff_ffff)),
194            Some((Sign::Positive, 0xf_ffff_ffff_ffff, -1_074))
195        );
196        assert_eq!(
197            decompose(f64::from_bits(1)),
198            Some((Sign::Positive, 1, -1_074))
199        );
200        assert_eq!(decompose(0.0), Some((Sign::Positive, 0, -1_074)));
201        assert_eq!(decompose(-0.0), Some((Sign::Negative, 0, -1_074)));
202    }
203
204    #[test]
205    fn decompose_not_finite() {
206        assert_eq!(decompose(f64::INFINITY), None);
207        assert_eq!(decompose(f64::NEG_INFINITY), None);
208        assert_eq!(decompose(f64::NAN), None);
209    }
210
211    #[test]
212    fn round_trips() {
213        for v in [
214            0.0,
215            -0.0,
216            1.0,
217            -1.5,
218            0.1,
219            9_007_199_254_740_992.0,
220            f64::MAX,
221            -f64::MAX,
222            f64::MIN_POSITIVE,
223            f64::from_bits(0xf_ffff_ffff_ffff),
224            f64::from_bits(1),
225            -f64::from_bits(1),
226        ] {
227            round_trip(v);
228        }
229    }
230
231    #[test]
232    fn exact() {
233        assert_eq!(compose(1, 0, false), 1.0);
234        assert_eq!(compose(3, -1, false), 1.5);
235        assert_eq!(compose(1 << 53, 0, false), 9_007_199_254_740_992.0);
236        assert_eq!(compose(1 << 100, -100, false), 1.0);
237        let zero = f64::from(Number::compose(Sign::Negative, 0, 5, false).unwrap());
238        assert_eq!(zero.to_bits(), 0x8000_0000_0000_0000);
239    }
240
241    #[test]
242    fn two_to_the_53_ties() {
243        let two_53 = 9_007_199_254_740_992.0;
244        // 2^53 + 1 is half way between 2^53 and 2^53 + 2: the even one.
245        assert_eq!(compose((1 << 53) + 1, 0, false), two_53);
246        // A sticky bit under it breaks the tie upward.
247        assert_eq!(compose(((1 << 53) + 1) << 8, -8, true), two_53 + 2.0);
248        // 2^53 + 3 is a tie whose even neighbor is the upper one.
249        assert_eq!(compose((1 << 53) + 3, 0, false), two_53 + 4.0);
250        // Under and over half round down and up: 2^53 + 0.5, 2^53 + 1.5.
251        assert_eq!(compose((1 << 54) + 1, -1, false), two_53);
252        assert_eq!(compose((1 << 54) + 3, -1, false), two_53 + 2.0);
253    }
254
255    #[test]
256    fn carry_out_of_the_significand() {
257        // 54 ones round up to 2^54: the significand overflows into the
258        // exponent.
259        assert_eq!(compose((1 << 54) - 1, 0, false), 18_014_398_509_481_984.0);
260        // 0b1111111101110111011101101101011011010110110100010000011, the
261        // 55-bit literal `number_coercion`'s tests pin: its two dropped
262        // bits are over half.
263        assert_eq!(
264            compose(0x7f_bbbb_6b6b_6883, 0, false),
265            35_953_735_677_732_996.0
266        );
267    }
268
269    #[test]
270    fn overflow() {
271        // Just under 2^1024 rounds up to it: Infinity.
272        assert_eq!(compose((1 << 54) - 1, 970, false), f64::INFINITY);
273        // 53 ones at the top exponent is the largest finite value.
274        assert_eq!(compose((1 << 53) - 1, 971, false), f64::MAX);
275        assert_eq!(compose(1, 1_024, false), f64::INFINITY);
276        assert_eq!(compose(1, i32::MAX, false), f64::INFINITY);
277        assert_eq!(compose(u128::MAX, i32::MAX, true), f64::INFINITY);
278        let negative = f64::from(Number::compose(Sign::Negative, 1, 1_024, false).unwrap());
279        assert_eq!(negative, f64::NEG_INFINITY);
280    }
281
282    /// A sticky bit under 53 significant bits or fewer has no rounded-away
283    /// bit to place the half, so the call is refused.
284    #[test]
285    fn narrow_sticky_refused() {
286        assert_eq!(Number::compose(Sign::Positive, 1, -1_074, true), None);
287        assert_eq!(Number::compose(Sign::Positive, 0, 0, true), None);
288        assert_eq!(
289            Number::compose(Sign::Positive, (1 << 53) - 1, 0, true),
290            None
291        );
292        assert!(Number::compose(Sign::Positive, 1 << 53, 0, true).is_some());
293    }
294
295    #[test]
296    fn subnormal() {
297        let min = f64::from_bits(1);
298        assert_eq!(compose(1, -1_074, false), min);
299        assert_eq!(compose(1 << 60, -1_134, false), min);
300        // Half the smallest subnormal is a tie with zero: zero is even.
301        assert_eq!(compose(1, -1_075, false), 0.0);
302        // A sticky bit under it rounds up.
303        assert_eq!(compose(1 << 60, -1_135, true), min);
304        // Three halves is a tie between one and two: two is even.
305        assert_eq!(compose(3, -1_075, false), f64::from_bits(2));
306        assert_eq!(compose(1, i32::MIN, false), 0.0);
307        // The largest subnormal plus a half rounds up into the normals.
308        assert_eq!(
309            compose(0x1f_ffff_ffff_ffff, -1_075, false),
310            f64::MIN_POSITIVE
311        );
312        // A subnormal significand at a normal exponent is a normal value.
313        assert_eq!(compose(1, -1_022, false), f64::MIN_POSITIVE);
314    }
315
316    /// Rust's own `u128 as f64` rounds to nearest even too: over widths up
317    /// to 128 bits and exponents scaling it by an exact power of two,
318    /// `compose` agrees with it.
319    #[test]
320    fn agrees_with_u128_as_f64() {
321        let mut state = 0x9e37_79b9_7f4a_7c15u64;
322        let mut next = || {
323            state ^= state << 13;
324            state ^= state >> 7;
325            state ^= state << 17;
326            state
327        };
328        for _ in 0..10_000 {
329            let wide = (u128::from(next()) << 64) | u128::from(next());
330            let significand = wide >> (next() % 128);
331            let exponent = (next() % 1_601) as i32 - 800;
332            let expected = (significand as f64) * 2f64.powi(exponent);
333            assert_eq!(
334                compose(significand, exponent, false),
335                expected,
336                "{significand:#x}"
337            );
338        }
339    }
340}
nanvm_lib/vm/number/binary64.rs

nanvm_lib/vm/number/
binary64.rs
