competitive/num/
dyadic.rs

1use super::{One, Signed, Zero};
2use std::{
3    cmp::Ordering,
4    fmt,
5    ops::{Add, Neg, Sub},
6};
7
8/// A reduced `numerator / 2^exponent`.
9#[derive(Clone, Copy, PartialEq, Eq, Hash)]
10pub struct Dyadic<T> {
11    numerator: T,
12    exponent: u32,
13}
14
15impl<T> Dyadic<T> {
16    pub fn numerator(self) -> T {
17        self.numerator
18    }
19
20    pub fn exponent(self) -> u32 {
21        self.exponent
22    }
23}
24
25impl<T> Dyadic<T>
26where
27    T: Signed,
28{
29    pub fn new(mut numerator: T, mut exponent: u32) -> Self {
30        let two = T::one() + T::one();
31        if numerator.is_zero() {
32            exponent = 0;
33        }
34        while exponent > 0 && (numerator % two).is_zero() {
35            numerator /= two;
36            exponent -= 1;
37        }
38        Self {
39            numerator,
40            exponent,
41        }
42    }
43
44    /// The simplest number strictly between the bounds; absent bounds are unbounded.
45    /// Returns `None` for an empty interval.
46    pub fn simplest_between(left: Option<Self>, right: Option<Self>) -> Option<Self> {
47        if left.zip(right).is_some_and(|(l, r)| l >= r) {
48            return None;
49        }
50        if left.is_none_or(|l| l.numerator < T::zero())
51            && right.is_none_or(|r| r.numerator > T::zero())
52        {
53            return Some(Self::zero());
54        }
55        let negative = right.is_some_and(|r| r.numerator <= T::zero());
56        let (lower, upper) = if negative {
57            (-right.unwrap(), left.map(Neg::neg))
58        } else {
59            (left.unwrap(), right)
60        };
61        let two = T::one() + T::one();
62        for exponent in 0.. {
63            let floor = if exponent < lower.exponent {
64                lower.numerator / two.pow(lower.exponent - exponent)
65            } else {
66                lower.numerator * two.pow(exponent - lower.exponent)
67            };
68            let value = Self::new(floor + T::one(), exponent);
69            if upper.is_none_or(|u| value < u) {
70                return Some(if negative { -value } else { value });
71            }
72        }
73        unreachable!()
74    }
75}
76
77impl<T> From<T> for Dyadic<T>
78where
79    T: Signed,
80{
81    fn from(numerator: T) -> Self {
82        Self::new(numerator, 0)
83    }
84}
85impl<T> Zero for Dyadic<T>
86where
87    T: Signed,
88{
89    fn zero() -> Self {
90        T::zero().into()
91    }
92}
93impl<T> One for Dyadic<T>
94where
95    T: Signed,
96{
97    fn one() -> Self {
98        T::one().into()
99    }
100}
101impl<T> Default for Dyadic<T>
102where
103    T: Signed,
104{
105    fn default() -> Self {
106        Self::zero()
107    }
108}
109impl<T> Neg for Dyadic<T>
110where
111    T: Signed,
112{
113    type Output = Self;
114    fn neg(self) -> Self {
115        Self {
116            numerator: -self.numerator,
117            exponent: self.exponent,
118        }
119    }
120}
121impl<T> Add for Dyadic<T>
122where
123    T: Signed,
124{
125    type Output = Self;
126    fn add(self, rhs: Self) -> Self {
127        let exponent = self.exponent.max(rhs.exponent);
128        let two = T::one() + T::one();
129        let left = self.numerator * two.pow(exponent - self.exponent);
130        let right = rhs.numerator * two.pow(exponent - rhs.exponent);
131        Self::new(left + right, exponent)
132    }
133}
134impl<T> Sub for Dyadic<T>
135where
136    T: Signed,
137{
138    type Output = Self;
139    fn sub(self, rhs: Self) -> Self {
140        self + -rhs
141    }
142}
143impl<T> Ord for Dyadic<T>
144where
145    T: Signed,
146{
147    fn cmp(&self, rhs: &Self) -> Ordering {
148        let exponent = self.exponent.max(rhs.exponent);
149        let two = T::one() + T::one();
150        let left = self.numerator * two.pow(exponent - self.exponent);
151        let right = rhs.numerator * two.pow(exponent - rhs.exponent);
152        left.cmp(&right)
153    }
154}
155impl<T> PartialOrd for Dyadic<T>
156where
157    T: Signed,
158{
159    fn partial_cmp(&self, rhs: &Self) -> Option<Ordering> {
160        Some(self.cmp(rhs))
161    }
162}
163impl<T> fmt::Display for Dyadic<T>
164where
165    T: Signed,
166{
167    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
168        if self.exponent == 0 {
169            write!(f, "{}", self.numerator)
170        } else {
171            let two = T::one() + T::one();
172            write!(f, "{}/{}", self.numerator, two.pow(self.exponent))
173        }
174    }
175}
176impl<T> fmt::Debug for Dyadic<T>
177where
178    T: Signed,
179{
180    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
181        fmt::Display::fmt(self, f)
182    }
183}
184
185#[cfg(kani)]
186mod verification {
187    use super::*;
188
189    macro_rules! proof_normalization {
190        ($name:ident, $ty:ty, $unwind:expr) => {
191            #[kani::proof]
192            #[kani::unwind($unwind)]
193            fn $name() {
194                let numerator: $ty = kani::any();
195                let exponent: u32 = kani::any();
196                let actual = Dyadic::new(numerator, exponent);
197                let expected = if numerator == 0 {
198                    (0, 0)
199                } else {
200                    let removed = exponent.min(numerator.trailing_zeros());
201                    (numerator >> removed, exponent - removed)
202                };
203                kani::assert(
204                    (actual.numerator(), actual.exponent()) == expected,
205                    concat!(
206                        "For every ",
207                        stringify!($ty),
208                        " numerator and u32 exponent, dyadic normalization divides out exactly the largest permitted power of two, preserving the value and reducing zero to 0/1",
209                    ),
210                );
211                kani::cover!(numerator == <$ty>::MIN && exponent == u32::MAX);
212                kani::cover!(numerator == <$ty>::MAX && exponent == u32::MAX);
213                kani::cover!(numerator == 0 && exponent == u32::MAX);
214                kani::cover!(numerator == <$ty>::MIN && exponent == 0);
215            }
216        };
217    }
218
219    proof_normalization!(proof_i8_normalization, i8, 9);
220    proof_normalization!(proof_i16_normalization, i16, 17);
221    proof_normalization!(proof_i32_normalization, i32, 33);
222}
223
224#[cfg(test)]
225mod tests {
226    use super::*;
227    use crate::tools::{Xorshift, testutil::integer_pairs};
228
229    #[test]
230    fn arithmetic_and_order() {
231        let mut rng = Xorshift::default();
232        const MAX_EXPONENT: u32 = i32::BITS / 2;
233        let bound = i32::MAX >> (MAX_EXPONENT + 1);
234        for exponent in 0..=MAX_EXPONENT {
235            for (left, right) in integer_pairs(-bound..=bound).chain(
236                rng.random_iter((-bound..=bound, -bound..=bound))
237                    .take(60_000),
238            ) {
239                for other_exponent in [0, exponent, MAX_EXPONENT] {
240                    let x = Dyadic::new(left, exponent);
241                    let y = Dyadic::new(right, other_exponent);
242                    let scaled = |v: Dyadic<i32>| v.numerator() << (MAX_EXPONENT - v.exponent());
243                    let a = left << (MAX_EXPONENT - exponent);
244                    let b = right << (MAX_EXPONENT - other_exponent);
245                    assert_eq!((scaled(x), scaled(y)), (a, b));
246                    assert_eq!(x.cmp(&y), a.cmp(&b));
247                    assert_eq!(scaled(x + y), a + b);
248                    assert_eq!(scaled(x - y), a - b);
249                    assert_eq!(scaled(-x), -a);
250                    for v in [x, y, x + y, x - y, -x] {
251                        assert!(v.exponent() == 0 || v.numerator() % 2 != 0);
252                    }
253                    let wide = Dyadic::new(i64::from(left), exponent)
254                        + Dyadic::new(i64::from(right), other_exponent);
255                    assert_eq!(
256                        (i64::from((x + y).numerator()), (x + y).exponent()),
257                        (wide.numerator(), wide.exponent())
258                    );
259                }
260            }
261        }
262    }
263
264    #[test]
265    fn simplest_between() {
266        // Bounds on a 1/8 grid. Enumerate by denominator, then distance from zero.
267        let bounds: Vec<_> = std::iter::once(None).chain((-32..=32).map(Some)).collect();
268        for &left in &bounds {
269            for &right in &bounds {
270                let actual = Dyadic::simplest_between(
271                    left.map(|n| Dyadic::new(n, 3)),
272                    right.map(|n| Dyadic::new(n, 3)),
273                );
274                let expected = (0..=4).find_map(|e| {
275                    (-80i32..=80)
276                        .filter(|&n| {
277                            left.is_none_or(|l| l * (1 << e) < n * 8)
278                                && right.is_none_or(|r| n * 8 < r * (1 << e))
279                        })
280                        .min_by_key(|n| n.abs())
281                        .map(|n| Dyadic::new(n, e))
282                });
283                assert_eq!(actual, expected, "{left:?} .. {right:?}");
284            }
285        }
286    }
287}
competitive/num/dyadic.rs

competitive/num/
dyadic.rs
