fixed_point/
fixed_point.rs

1// Copyright 2016 The Fuchsia Authors
2// Copyright (c) 2013, Google Inc. All rights reserved.
3//
4// Use of this source code is governed by a MIT-style
5// license that can be found in the LICENSE file or at
6// https://opensource.org/licenses/MIT
7
8//! 32.64 unsigned fixed point arithmetic, mirroring `<lib/fixed_point.h>`.
9
10use core::ptr::NonNull;
11
12use crate::fixed_point_debug::{
13    debug_mul_u32_u32, debug_u32_mul_u64_fp32_64, debug_u64_mul_u32_fp32_64,
14    debug_u64_mul_u64_fp32_64,
15};
16
17/// A 32.64 unsigned fixed point number: 32 integer bits followed by 64
18/// fractional bits, stored as three little-endian-by-significance limbs.
19///
20/// Layout-compatible with the C `struct fp_32_64`.
21#[repr(C)]
22#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)]
23pub struct Fp3264 {
24    /// unshifted value
25    pub l0: u32,
26    /// value shifted left 32 bits (or bit -1 to -32)
27    pub l32: u32,
28    /// value shifted left 64 bits (or bit -33 to -64)
29    pub l64: u32,
30}
31
32zr::static_assert_size_and_align!(Fp3264, 12, 4);
33
34/// Computes `dividend / divisor` as a 32.64 fixed point value, truncating
35/// any fraction below the 64th fractional bit, and stores the result into `result`.
36///
37/// # Panics
38///
39/// Panics if `divisor` is zero (undefined behaviour in the C original).
40///
41/// # Safety
42///
43/// `result` must be valid for writes and properly aligned.
44#[inline]
45pub unsafe fn fp_32_64_div_32_32(result: NonNull<Fp3264>, dividend: u32, divisor: u32) {
46    let mut tmp = (u64::from(dividend) << 32) / u64::from(divisor);
47    let rem = ((u64::from(dividend) << 32) % u64::from(divisor)) as u32;
48    // SAFETY: Caller guarantees `result` is properly aligned and valid for writes.
49    unsafe {
50        let res = result.as_ptr();
51        (*res).l0 = (tmp >> 32) as u32;
52        (*res).l32 = tmp as u32;
53        tmp = (u64::from(rem) << 32) / u64::from(divisor);
54        (*res).l64 = tmp as u32;
55    }
56}
57
58/// Multiplies two 32-bit operands into a 64-bit product. `a_shift` and
59/// `b_shift` describe the bit position of each operand relative to the binary
60/// point and are used purely for debug tracing.
61#[inline]
62pub(crate) fn mul_u32_u32(a: u32, b: u32, a_shift: i32, b_shift: i32) -> u64 {
63    let ret = u64::from(a) * u64::from(b);
64    debug_mul_u32_u32(a, b, a_shift, b_shift, ret);
65    ret
66}
67
68/// Multiplies a 32-bit integer `a` by the 32.64 fixed point value `b`,
69/// returning the product rounded to the nearest 64-bit integer.
70#[inline]
71pub fn u64_mul_u32_fp32_64(a: u32, b: Fp3264) -> u64 {
72    let mut res_0 = mul_u32_u32(a, b.l0, 0, 0);
73    let tmp = mul_u32_u32(a, b.l32, 0, -32);
74    res_0 = res_0.wrapping_add(tmp >> 32);
75    let mut res_l32 = u64::from(tmp as u32);
76    res_l32 = res_l32.wrapping_add(mul_u32_u32(a, b.l64, 0, -64) >> 32); // Improve rounding accuracy
77    res_0 = res_0.wrapping_add(res_l32 >> 32);
78    let res_l32_32 = res_l32 as u32;
79    let ret = res_0.wrapping_add(u64::from(res_l32_32 >> 31)); // Round to nearest integer
80
81    debug_u64_mul_u32_fp32_64(a, b, res_0, res_l32_32, ret);
82
83    ret
84}
85
86/// Multiplies a 64-bit integer `a` by the 32.64 fixed point value `b`,
87/// returning the low 32 bits of the product rounded to the nearest integer.
88#[inline]
89pub fn u32_mul_u64_fp32_64(a: u64, b: Fp3264) -> u32 {
90    let a_r32 = (a >> 32) as u32;
91    let a_0 = a as u32;
92
93    // mul_u32_u32(a_r32, b.l0, 32, 0) does not affect result
94    let mut res_l32 = mul_u32_u32(a_0, b.l0, 0, 0) << 32;
95    res_l32 = res_l32.wrapping_add(mul_u32_u32(a_r32, b.l32, 32, -32) << 32);
96    res_l32 = res_l32.wrapping_add(mul_u32_u32(a_0, b.l32, 0, -32));
97    res_l32 = res_l32.wrapping_add(mul_u32_u32(a_r32, b.l64, 32, -64));
98    res_l32 = res_l32.wrapping_add(mul_u32_u32(a_0, b.l64, 0, -64) >> 32); // Improve rounding accuracy
99    let ret = (res_l32 >> 32).wrapping_add(u64::from((res_l32 as u32) >> 31)) as u32; // Round to nearest integer
100
101    debug_u32_mul_u64_fp32_64(a, b, res_l32, ret);
102
103    ret
104}
105
106/// Multiplies a 64-bit integer `a` by the 32.64 fixed point value `b`,
107/// returning the low 64 bits of the product rounded to the nearest integer.
108#[inline]
109pub fn u64_mul_u64_fp32_64(a: u64, b: Fp3264) -> u64 {
110    let a_r32 = (a >> 32) as u32;
111    let a_0 = a as u32;
112
113    let mut tmp = mul_u32_u32(a_r32, b.l0, 32, 0);
114    let mut res_0 = tmp << 32;
115    tmp = mul_u32_u32(a_0, b.l0, 0, 0);
116    res_0 = res_0.wrapping_add(tmp);
117    tmp = mul_u32_u32(a_r32, b.l32, 32, -32);
118    res_0 = res_0.wrapping_add(tmp);
119    tmp = mul_u32_u32(a_0, b.l32, 0, -32);
120    res_0 = res_0.wrapping_add(tmp >> 32);
121    let mut res_l32 = u64::from(tmp as u32);
122    tmp = mul_u32_u32(a_r32, b.l64, 32, -64);
123    res_0 = res_0.wrapping_add(tmp >> 32);
124    res_l32 = res_l32.wrapping_add(u64::from(tmp as u32));
125    tmp = mul_u32_u32(a_0, b.l64, 0, -64); // Improve rounding accuracy
126    res_l32 = res_l32.wrapping_add(tmp >> 32);
127    res_0 = res_0.wrapping_add(res_l32 >> 32);
128    let res_l32_32 = res_l32 as u32;
129    let ret = res_0.wrapping_add(u64::from(res_l32_32 >> 31)); // Round to nearest integer
130
131    debug_u64_mul_u64_fp32_64(a, b, res_0, res_l32_32, ret);
132
133    ret
134}
135
136/// FFI trampoline for C++ `fp_32_64_div_32_32`.
137///
138/// Called from C via `<lib/fixed_point.h>`; the C caller must pass a valid,
139/// aligned, non-null pointer to a `struct fp_32_64`.
140///
141/// # Safety
142///
143/// `result` must be valid for writes and properly aligned.
144#[unsafe(no_mangle)]
145pub unsafe extern "C" fn rust_fixed_point_fp_32_64_div_32_32(
146    result: NonNull<Fp3264>,
147    dividend: u32,
148    divisor: u32,
149) {
150    // SAFETY: Forwarding caller's safety contract to `fp_32_64_div_32_32`.
151    unsafe { fp_32_64_div_32_32(result, dividend, divisor) }
152}
153
154/// FFI trampoline for C++ `u64_mul_u32_fp32_64`.
155///
156/// Called from C via `<lib/fixed_point.h>`; `b` is passed by value matching C++.
157#[unsafe(no_mangle)]
158pub extern "C" fn rust_fixed_point_u64_mul_u32_fp32_64(a: u32, b: Fp3264) -> u64 {
159    u64_mul_u32_fp32_64(a, b)
160}
161
162/// FFI trampoline for C++ `u32_mul_u64_fp32_64`.
163///
164/// Called from C via `<lib/fixed_point.h>`; `b` is passed by value matching C++.
165#[unsafe(no_mangle)]
166pub extern "C" fn rust_fixed_point_u32_mul_u64_fp32_64(a: u64, b: Fp3264) -> u32 {
167    u32_mul_u64_fp32_64(a, b)
168}
169
170/// FFI trampoline for C++ `u64_mul_u64_fp32_64`.
171///
172/// Called from C via `<lib/fixed_point.h>`; `b` is passed by value matching C++.
173#[unsafe(no_mangle)]
174pub extern "C" fn rust_fixed_point_u64_mul_u64_fp32_64(a: u64, b: Fp3264) -> u64 {
175    u64_mul_u64_fp32_64(a, b)
176}
177
178#[cfg(test)]
179mod tests {
180    use super::*;
181
182    const ONE: Fp3264 = Fp3264 { l0: 1, l32: 0, l64: 0 };
183    const HALF: Fp3264 = Fp3264 { l0: 0, l32: 0x8000_0000, l64: 0 };
184    const MAX: Fp3264 = Fp3264 { l0: u32::MAX, l32: u32::MAX, l64: u32::MAX };
185
186    /// The 96-bit integer `l0:l32:l64` (i.e. the fixed point value scaled by 2^64).
187    fn raw(b: Fp3264) -> u128 {
188        (u128::from(b.l0) << 64) | (u128::from(b.l32) << 32) | u128::from(b.l64)
189    }
190
191    /// Reference for `fp_32_64_div_32_32`: floor(dividend * 2^64 / divisor) split into limbs.
192    fn div_reference(dividend: u32, divisor: u32) -> Fp3264 {
193        let q = (u128::from(dividend) << 64) / u128::from(divisor);
194        Fp3264 { l0: (q >> 64) as u32, l32: (q >> 32) as u32, l64: q as u32 }
195    }
196
197    /// Bit-exact reference model of the multiply routines.
198    ///
199    /// All three multiply helpers compute `a * b` exactly except that the low 32 bits of the
200    /// `a_0 * b.l64` partial product are discarded before rounding ("Improve rounding accuracy"
201    /// keeps only the high half). The intermediate `X` therefore equals the exact product scaled
202    /// by 2^-32 minus that discarded term, and the result is `X` rounded to nearest at bit 32.
203    /// Wrapping arithmetic mirrors the C unsigned semantics for inputs that overflow.
204    fn mul_reference(a: u64, b: Fp3264) -> u128 {
205        let a_0 = a as u32;
206        let discarded = (u64::from(a_0) * u64::from(b.l64)) & 0xffff_ffff;
207        let x = (u128::from(a).wrapping_mul(raw(b)).wrapping_sub(u128::from(discarded))) >> 32;
208        x.wrapping_add(1 << 31) >> 32
209    }
210
211    /// Mathematically exact `a * b` rounded to nearest integer.
212    ///
213    /// The full product is up to 160 bits wide, so this is computed modulo 2^128; that still
214    /// determines the rounded value modulo 2^64, which is all [`within_one_mod_u64`] needs.
215    fn exact_rounded(a: u64, b: Fp3264) -> u128 {
216        u128::from(a).wrapping_mul(raw(b)).wrapping_add(1u128 << 63) >> 64
217    }
218
219    /// Returns true if `got` is `exact` or `exact - 1`, compared modulo 2^64 so that the check
220    /// stays meaningful when the true product does not fit in the 64-bit result.
221    fn within_one_mod_u64(exact: u128, got: u64) -> bool {
222        (exact.wrapping_sub(u128::from(got)) & u128::from(u64::MAX)) <= 1
223    }
224
225    /// Deterministic pseudo-random sequence (64-bit LCG) for cross-check vectors.
226    fn next_random(state: &mut u64) -> u64 {
227        *state = state.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
228        *state ^ (*state >> 29)
229    }
230
231    #[test]
232    fn div_known_values() {
233        let mut result = Fp3264::default();
234
235        // SAFETY: `result` is a valid, aligned mutable reference.
236        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 1, 1) };
237        assert_eq!(result, ONE);
238
239        // SAFETY: `result` is a valid, aligned mutable reference.
240        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 1, 2) };
241        assert_eq!(result, HALF);
242
243        // SAFETY: `result` is a valid, aligned mutable reference.
244        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 3, 2) };
245        assert_eq!(result, Fp3264 { l0: 1, l32: 0x8000_0000, l64: 0 });
246
247        // SAFETY: `result` is a valid, aligned mutable reference.
248        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 1, 3) };
249        assert_eq!(result, Fp3264 { l0: 0, l32: 0x5555_5555, l64: 0x5555_5555 });
250
251        // SAFETY: `result` is a valid, aligned mutable reference.
252        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 2, 3) };
253        assert_eq!(result, Fp3264 { l0: 0, l32: 0xaaaa_aaaa, l64: 0xaaaa_aaaa });
254
255        // SAFETY: `result` is a valid, aligned mutable reference.
256        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 1, 4) };
257        assert_eq!(result, Fp3264 { l0: 0, l32: 0x4000_0000, l64: 0 });
258
259        // SAFETY: `result` is a valid, aligned mutable reference.
260        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 0, 7) };
261        assert_eq!(result, Fp3264::default());
262
263        // SAFETY: `result` is a valid, aligned mutable reference.
264        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), u32::MAX, 1) };
265        assert_eq!(result, Fp3264 { l0: u32::MAX, l32: 0, l64: 0 });
266
267        // SAFETY: `result` is a valid, aligned mutable reference.
268        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 1, u32::MAX) };
269        assert_eq!(result, Fp3264 { l0: 0, l32: 1, l64: 1 });
270
271        // SAFETY: `result` is a valid, aligned mutable reference.
272        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), u32::MAX, u32::MAX) };
273        assert_eq!(result, ONE);
274
275        // 1_000_000 / 3_000_000 == 1/3: the "ns per TSC tick" ratio for a 3GHz TSC.
276        let mut expected = Fp3264::default();
277        // SAFETY: `expected` is a valid, aligned mutable reference.
278        unsafe { fp_32_64_div_32_32(NonNull::from(&mut expected), 1, 3) };
279        // SAFETY: `result` is a valid, aligned mutable reference.
280        unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), 1_000_000, 3_000_000) };
281        assert_eq!(result, expected);
282    }
283
284    #[test]
285    fn div_matches_u128_reference() {
286        let mut state = 0x1234_5678_9abc_def0;
287        let mut result = Fp3264::default();
288        for _ in 0..10_000 {
289            let dividend = next_random(&mut state) as u32;
290            let divisor = (next_random(&mut state) as u32).max(1);
291            // SAFETY: `result` is a valid, aligned mutable reference.
292            unsafe { fp_32_64_div_32_32(NonNull::from(&mut result), dividend, divisor) };
293            assert_eq!(result, div_reference(dividend, divisor), "{dividend} / {divisor}");
294        }
295    }
296
297    #[test]
298    fn mul_u32_u32_products() {
299        assert_eq!(mul_u32_u32(0, u32::MAX, 0, 0), 0);
300        assert_eq!(mul_u32_u32(1, 1, 0, 0), 1);
301        assert_eq!(mul_u32_u32(0x1_0000, 0x1_0000, 0, -32), 1 << 32);
302        assert_eq!(mul_u32_u32(u32::MAX, u32::MAX, 32, -64), 0xffff_fffe_0000_0001);
303    }
304
305    #[test]
306    fn mul_by_zero() {
307        let zero = Fp3264::default();
308        assert_eq!(u64_mul_u32_fp32_64(u32::MAX, zero), 0);
309        assert_eq!(u32_mul_u64_fp32_64(u64::MAX, zero), 0);
310        assert_eq!(u64_mul_u64_fp32_64(u64::MAX, zero), 0);
311        assert_eq!(u64_mul_u32_fp32_64(0, MAX), 0);
312        assert_eq!(u32_mul_u64_fp32_64(0, MAX), 0);
313        assert_eq!(u64_mul_u64_fp32_64(0, MAX), 0);
314    }
315
316    #[test]
317    fn mul_by_one_is_identity() {
318        for a in [0u64, 1, 2, 0x7fff_ffff, 0x8000_0000, 0xffff_ffff, 0x1_0000_0000, u64::MAX] {
319            assert_eq!(u64_mul_u64_fp32_64(a, ONE), a);
320            // The 32-bit result keeps the low 32 bits of the product.
321            assert_eq!(u32_mul_u64_fp32_64(a, ONE), a as u32);
322            if let Ok(a32) = u32::try_from(a) {
323                assert_eq!(u64_mul_u32_fp32_64(a32, ONE), a);
324            }
325        }
326    }
327
328    #[test]
329    fn mul_rounds_to_nearest() {
330        // 0.5 * n rounds half up: 1*0.5 -> 1, 2*0.5 -> 1, 3*0.5 -> 2.
331        assert_eq!(u64_mul_u32_fp32_64(1, HALF), 1);
332        assert_eq!(u64_mul_u32_fp32_64(2, HALF), 1);
333        assert_eq!(u64_mul_u32_fp32_64(3, HALF), 2);
334        assert_eq!(u32_mul_u64_fp32_64(1, HALF), 1);
335        assert_eq!(u32_mul_u64_fp32_64(2, HALF), 1);
336        assert_eq!(u32_mul_u64_fp32_64(3, HALF), 2);
337        assert_eq!(u64_mul_u64_fp32_64(1, HALF), 1);
338        assert_eq!(u64_mul_u64_fp32_64(2, HALF), 1);
339        assert_eq!(u64_mul_u64_fp32_64(3, HALF), 2);
340
341        // Just under one half (0.5 - 2^-64) rounds down for a == 1...
342        let just_under_half = Fp3264 { l0: 0, l32: 0x7fff_ffff, l64: u32::MAX };
343        assert_eq!(u64_mul_u32_fp32_64(1, just_under_half), 0);
344        assert_eq!(u32_mul_u64_fp32_64(1, just_under_half), 0);
345        assert_eq!(u64_mul_u64_fp32_64(1, just_under_half), 0);
346        // ...and 3 * (0.5 - 2^-64) = 1.4999... rounds to 1.
347        assert_eq!(u64_mul_u32_fp32_64(3, just_under_half), 1);
348        assert_eq!(u32_mul_u64_fp32_64(3, just_under_half), 1);
349        assert_eq!(u64_mul_u64_fp32_64(3, just_under_half), 1);
350
351        // 1/3 scaled by 3_000_000_000 ticks (one second at 3GHz) is 1e9 ns.
352        let mut ns_per_tsc = Fp3264::default();
353        // SAFETY: `ns_per_tsc` is a valid, aligned mutable reference.
354        unsafe { fp_32_64_div_32_32(NonNull::from(&mut ns_per_tsc), 1_000_000, 3_000_000) };
355        assert_eq!(u64_mul_u64_fp32_64(3_000_000_000, ns_per_tsc), 1_000_000_000);
356        assert_eq!(u32_mul_u64_fp32_64(3_000_000_000, ns_per_tsc), 1_000_000_000);
357        assert_eq!(u64_mul_u32_fp32_64(3_000_000_000, ns_per_tsc), 1_000_000_000);
358        assert_eq!(u64_mul_u64_fp32_64(1, ns_per_tsc), 0);
359        assert_eq!(u64_mul_u64_fp32_64(2, ns_per_tsc), 1);
360    }
361
362    #[test]
363    fn mul_max_values() {
364        // u32::MAX * (2^32 - 2^-64) = 2^64 - 2^32 - (2^32 - 1) * 2^-64, which rounds to
365        // 2^64 - 2^32 and comfortably fits in a u64.
366        assert_eq!(u64_mul_u32_fp32_64(u32::MAX, MAX), 0xffff_ffff_0000_0000);
367        assert_eq!(
368            u64_mul_u32_fp32_64(u32::MAX, MAX),
369            mul_reference(u64::from(u32::MAX), MAX) as u64
370        );
371
372        // u64::MAX * (2^32 - 2^-64) overflows; the C routines wrap modulo 2^64 / 2^32.
373        assert_eq!(u64_mul_u64_fp32_64(u64::MAX, MAX), mul_reference(u64::MAX, MAX) as u64);
374        assert_eq!(u32_mul_u64_fp32_64(u64::MAX, MAX), mul_reference(u64::MAX, MAX) as u32);
375        assert!(within_one_mod_u64(
376            exact_rounded(u64::MAX, MAX),
377            u64_mul_u64_fp32_64(u64::MAX, MAX)
378        ));
379
380        // Largest non-overflowing 64-bit case: 2^32 * (2^32 - 1) = 2^64 - 2^32.
381        let int_max = Fp3264 { l0: u32::MAX, l32: 0, l64: 0 };
382        assert_eq!(u64_mul_u64_fp32_64(1 << 32, int_max), 0xffff_ffff_0000_0000);
383    }
384
385    #[test]
386    fn mul_matches_u128_reference() {
387        let mut state = 0xdead_beef_cafe_f00d;
388        for i in 0..20_000 {
389            let a = next_random(&mut state);
390            let b = Fp3264 {
391                l0: next_random(&mut state) as u32,
392                l32: next_random(&mut state) as u32,
393                l64: next_random(&mut state) as u32,
394            };
395            // Mix in some small operands so that the products do not always overflow.
396            let (a, b) = match i % 4 {
397                0 => (a, b),
398                1 => (a >> 40, b),
399                2 => (a, Fp3264 { l0: b.l0 >> 24, ..b }),
400                _ => (a >> 48, Fp3264 { l0: 0, ..b }),
401            };
402            let a32 = a as u32;
403
404            assert_eq!(u64_mul_u64_fp32_64(a, b), mul_reference(a, b) as u64, "{a:#x} * {b:x?}");
405            assert_eq!(u32_mul_u64_fp32_64(a, b), mul_reference(a, b) as u32, "{a:#x} * {b:x?}");
406            assert_eq!(
407                u64_mul_u32_fp32_64(a32, b),
408                mul_reference(u64::from(a32), b) as u64,
409                "{a32:#x} * {b:x?}"
410            );
411
412            // Independently of the bit-exact model, the truncated `a_0 * l64` term can perturb the
413            // result by at most one unit below the exactly rounded product.
414            assert!(
415                within_one_mod_u64(exact_rounded(a, b), u64_mul_u64_fp32_64(a, b)),
416                "{a:#x} * {b:x?}"
417            );
418            assert!(
419                within_one_mod_u64(exact_rounded(u64::from(a32), b), u64_mul_u32_fp32_64(a32, b)),
420                "{a32:#x} * {b:x?}"
421            );
422        }
423    }
424
425    #[test]
426    fn ffi_trampolines() {
427        let mut result = Fp3264::default();
428        // SAFETY: `result` is a valid, aligned mutable reference.
429        unsafe { rust_fixed_point_fp_32_64_div_32_32(NonNull::from(&mut result), 1, 3) };
430        let mut expected = Fp3264::default();
431        // SAFETY: `expected` is a valid, aligned mutable reference.
432        unsafe { fp_32_64_div_32_32(NonNull::from(&mut expected), 1, 3) };
433        assert_eq!(result, expected);
434        assert_eq!(rust_fixed_point_u64_mul_u32_fp32_64(3, result), u64_mul_u32_fp32_64(3, result));
435        assert_eq!(
436            rust_fixed_point_u32_mul_u64_fp32_64(3_000_000_000, result),
437            u32_mul_u64_fp32_64(3_000_000_000, result)
438        );
439        assert_eq!(
440            rust_fixed_point_u64_mul_u64_fp32_64(3_000_000_000, result),
441            u64_mul_u64_fp32_64(3_000_000_000, result)
442        );
443        assert_eq!(rust_fixed_point_u64_mul_u64_fp32_64(3_000_000_000, result), 1_000_000_000);
444    }
445}
fixed_point/fixed_point.rs

fixed_point/
fixed_point.rs
