qfall_math/rational/q/
rounding.rs

1// Copyright 2023 Marvin Beckmann
2//
3// This file is part of qFALL-math.
4//
5// qFALL-math is free software: you can redistribute it and/or modify it under
6// the terms of the Mozilla Public License Version 2.0 as published by the
7// Mozilla Foundation. See <https://mozilla.org/en-US/MPL/2.0/>.
8
9//! This module includes functionality for rounding instances of [`Q`].
10
11use super::Q;
12use crate::{error::MathError, integer::Z, macros::seeded::seedable_function, traits::Distance};
13use flint3_sys::{
14    {fmpq_sgn, fmpq_simplest_between}, {fmpz_cdiv_q, fmpz_fdiv_q},
15};
16
17impl Q {
18    /// Rounds the given rational [`Q`] down to the next integer [`Z`].
19    ///
20    /// # Examples
21    /// ```
22    /// use qfall_math::rational::Q;
23    /// use qfall_math::integer::Z;
24    ///
25    /// let value = Q::from((5, 2));
26    /// assert_eq!(Z::from(2), value.floor());
27    ///
28    /// let value = Q::from((-5, 2));
29    /// assert_eq!(Z::from(-3), value.floor());
30    ///
31    /// let value = Q::from(2);
32    /// assert_eq!(Z::from(2), value.floor());
33    /// ```
34    pub fn floor(&self) -> Z {
35        let mut out = Z::default();
36        unsafe { fmpz_fdiv_q(&mut out.value, &self.value.num, &self.value.den) };
37        out
38    }
39
40    /// Rounds the given rational [`Q`] up to the next integer [`Z`].
41    ///
42    /// # Examples
43    /// ```
44    /// use qfall_math::rational::Q;
45    /// use qfall_math::integer::Z;
46    ///
47    /// let value = Q::from((5, 2));
48    /// assert_eq!(Z::from(3), value.ceil());
49    ///
50    /// let value = Q::from((-5, 2));
51    /// assert_eq!(Z::from(-2), value.ceil());
52    ///
53    /// let value = Q::from(2);
54    /// assert_eq!(Z::from(2), value.ceil());
55    /// ```
56    pub fn ceil(&self) -> Z {
57        let mut out = Z::default();
58        unsafe { fmpz_cdiv_q(&mut out.value, &self.value.num, &self.value.den) };
59        out
60    }
61
62    /// Rounds the given rational [`Q`] to the closest integer [`Z`].
63    /// If the distance is equal, it rounds up.
64    ///
65    /// # Examples
66    /// ```
67    /// use qfall_math::rational::Q;
68    /// use qfall_math::integer::Z;
69    ///
70    /// let value = Q::from((5, 2));
71    /// assert_eq!(Z::from(3), value.round());
72    ///
73    /// let value = Q::from((-5, 2));
74    /// assert_eq!(Z::from(-2), value.round());
75    ///
76    /// let value = Q::from(2);
77    /// assert_eq!(Z::from(2), value.round());
78    /// ```
79    pub fn round(&self) -> Z {
80        if Q::from(self.floor()).distance(self) < Q::from(0.5) {
81            self.floor()
82        } else {
83            self.ceil()
84        }
85    }
86
87    /// Returns the smallest rational with the smallest denominator in the range
88    /// `\[self - |precision|, self + |precision|\]`.
89    ///
90    /// This function allows to free memory in exchange for the specified loss of
91    /// precision (see Example 3). Be aware that this loss of precision is propagated by
92    /// arithmetic operations and can be significantly increased depending on the
93    /// performed operations.
94    ///
95    /// This function ensures that simplifying does not change the sign of `self`.
96    ///
97    /// Parameters:
98    /// - `precision`: the precision the new value can differ from `self`.
99    ///   Note that the absolute value is relevant, not the sign.
100    ///
101    /// Returns the simplest [`Q`] within the defined range.
102    ///
103    /// # Examples
104    /// ```
105    /// use qfall_math::rational::Q;
106    ///
107    /// let value = Q::from((17, 20));
108    /// let precision = Q::from((1, 20));
109    ///
110    /// let simplified = Q::from((4, 5));
111    /// assert_eq!(simplified, value.simplify(&precision));
112    /// ```
113    ///
114    /// ```
115    /// use qfall_math::rational::Q;
116    ///
117    /// let value = Q::from((3, 2));
118    ///
119    /// assert_eq!(Q::ONE, value.simplify(0.5));
120    /// ```
121    ///
122    /// ## Simplify with reasonable precision loss
123    /// This example uses [`Q::INV_MAX32`], i.e. a loss of precision of at most `1 / 2^31 - 2` behind the decimal point.
124    /// If you require higher precision, [`Q::INV_MAX62`] is available.
125    /// ```
126    /// use qfall_math::rational::Q;
127    /// let value = Q::PI;
128    ///
129    /// let simplified = value.simplify(Q::INV_MAX32);
130    ///
131    /// assert_ne!(&Q::PI, &simplified);
132    /// assert!(&simplified >= &(Q::PI - Q::INV_MAX32));
133    /// assert!(&simplified <= &(Q::PI + Q::INV_MAX32));
134    /// ```
135    pub fn simplify(&self, precision: impl Into<Q>) -> Self {
136        let precision = precision.into();
137
138        let lower = self - &precision;
139        let upper = self + &precision;
140        let mut out = Q::default();
141        unsafe { fmpq_simplest_between(&mut out.value, &lower.value, &upper.value) };
142
143        if unsafe { fmpq_sgn(&self.value) != fmpq_sgn(&out.value) } {
144            return Q::MINUS_ONE * out;
145        }
146
147        out
148    }
149
150    seedable_function!(
151        /// Performs the randomized rounding algorithm
152        /// by sampling from a discrete Gaussian over the integers shifted
153        /// by `self` with gaussian parameter `r`.
154        ///
155        /// Parameters:
156        /// - `r`: specifies the Gaussian parameter, which is proportional
157        ///   to the standard deviation `sigma * sqrt(2 * pi) = r`
158        #[seeded]
159        /// - `seed`: specifies the 256-bit seed of the PRNG used for sampling
160        #[optionally_seeded]
161        /// - `seed`: specifies an optional 256-bit seed of the PRNG used for sampling.
162        #[optionally_seeded]
163        ///   If `None` is provided, a fresh [`ThreadRng`](rand::rngs::ThreadRng) is used instead.
164        ///
165        /// Returns the rounded value as an [`Z`] or an error if `r < 0`.
166        ///
167        /// # Examples
168        /// ```
169        /// use qfall_math::rational::Q;
170        ///
171        /// let value = Q::from((5, 2));
172        #[unseeded]
173        /// let rounded = value.randomized_rounding(3).unwrap();
174        #[seeded]
175        /// let rounded = value.randomized_rounding_seeded(3, [42; 32]).unwrap();
176        /// ```
177        ///
178        /// # Errors and Failures
179        /// - Returns a [`MathError`] of type [`InvalidIntegerInput`](MathError::InvalidIntegerInput)
180        ///   if `r < 0`.
181        ///
182        /// This function implements randomized rounding according to:
183        /// - Peikert, C. (2010, August).
184        ///   An efficient and parallel Gaussian sampler for lattices.
185        ///   In: Annual Cryptology Conference (pp. 80-97).
186        ///   <https://link.springer.com/chapter/10.1007/978-3-642-14623-7_5>
187        pub(crate) fn randomized_rounding(
188            &self,
189            r: impl Into<Q>,
190            seed: Option<[u8; 32]>,
191        ) -> Result<Z, MathError> {
192            Z::sample_discrete_gauss_optionally_seeded(self, r, seed)
193        }
194    );
195}
196
197#[cfg(test)]
198mod test_floor {
199    use crate::{integer::Z, rational::Q};
200
201    // Ensure that positive rationals are rounded correctly
202    #[test]
203    fn positive() {
204        let val_1 = Q::from((i64::MAX, 2));
205        let val_2 = Q::from((1, u64::MAX));
206
207        assert_eq!(Z::from((i64::MAX - 1) / 2), val_1.floor());
208        assert_eq!(Z::ZERO, val_2.floor());
209    }
210
211    // Ensure that negative rationals are rounded correctly
212    #[test]
213    fn negative() {
214        let val_1 = Q::from((-i64::MAX, 2));
215        let val_2 = Q::from((-1, u64::MAX));
216
217        assert_eq!(Z::from((-i64::MAX - 1) / 2), val_1.floor());
218        assert_eq!(Z::MINUS_ONE, val_2.floor());
219    }
220}
221
222#[cfg(test)]
223mod test_ceil {
224    use crate::{integer::Z, rational::Q};
225
226    // Ensure that positive rationals are rounded correctly
227    #[test]
228    fn positive() {
229        let val_1 = Q::from((i64::MAX, 2));
230        let val_2 = Q::from((1, u64::MAX));
231
232        assert_eq!(Z::from((i64::MAX - 1) / 2 + 1), val_1.ceil());
233        assert_eq!(Z::ONE, val_2.ceil());
234    }
235
236    // Ensure that negative rationals are rounded correctly
237    #[test]
238    fn negative() {
239        let val_1 = Q::from((-i64::MAX, 2));
240        let val_2 = Q::from((-1, u64::MAX));
241
242        assert_eq!(Z::from((-i64::MAX - 1) / 2 + 1), val_1.ceil());
243        assert_eq!(Z::ZERO, val_2.ceil());
244    }
245}
246
247#[cfg(test)]
248mod test_round {
249    use crate::{integer::Z, rational::Q};
250
251    // Ensure that positive rationals are rounded correctly
252    #[test]
253    fn positive() {
254        let val_1 = Q::from((i64::MAX, 2));
255        let val_2 = Q::from((1, u64::MAX));
256
257        assert_eq!(Z::from((i64::MAX - 1) / 2 + 1), val_1.round());
258        assert_eq!(Z::ZERO, val_2.round());
259    }
260
261    // Ensure that negative rationals are rounded correctly
262    #[test]
263    fn negative() {
264        let val_1 = Q::from((-i64::MAX, 2));
265        let val_2 = Q::from((-1, u64::MAX));
266
267        assert_eq!(Z::from((-i64::MAX - 1) / 2 + 1), val_1.round());
268        assert_eq!(Z::ZERO, val_2.round());
269    }
270}
271
272#[cfg(test)]
273mod test_simplify {
274    use crate::{integer::Z, rational::Q, traits::Distance};
275
276    /// Ensure that negative precision works as expected
277    #[test]
278    fn precision_absolute_value() {
279        let value_1 = Q::from((17, 20));
280        let value_2 = Q::from((-17, 20));
281        let precision = Q::from((-1, 20));
282
283        let simplified_1 = Q::from((4, 5));
284        let simplified_2 = Q::from((-4, 5));
285        assert_eq!(simplified_1, value_1.simplify(&precision));
286        assert_eq!(simplified_2, value_2.simplify(precision));
287    }
288
289    /// Ensure that large values with pointer representations are reduced
290    #[test]
291    fn large_pointer_representation() {
292        let value = Q::from((u64::MAX - 1, u64::MAX));
293        let precision = Q::from((1, u64::MAX));
294
295        assert_eq!(Q::ONE, value.simplify(&precision));
296    }
297
298    /// Ensure that the simplified value stays in range
299    #[test]
300    fn stay_in_precision() {
301        let value = Q::from((i64::MAX - 1, i64::MAX));
302        let precision = Q::from((1, u64::MAX - 1));
303
304        let simplified = value.simplify(&precision);
305        assert!(&value - &precision <= simplified && simplified <= &value + &precision);
306        assert!(Q::from((i64::MAX - 2, i64::MAX)) <= simplified && simplified <= 1);
307    }
308
309    /// Ensure max_bits of denominator are not bigger than 1/2 * precision
310    #[test]
311    fn max_bits_denominator() {
312        let value = Q::PI;
313        let precisions = [Q::INV_MAX8, Q::INV_MAX16, Q::INV_MAX32];
314
315        for precision in precisions {
316            let inv_precision = precision.get_denominator();
317            let inv_precision = inv_precision.div_ceil(2);
318
319            let simplified = value.simplify(&precision);
320            let denominator = simplified.get_denominator();
321
322            assert!(denominator.distance(Z::ZERO) < inv_precision);
323        }
324    }
325
326    /// Ensure that a value which can not be simplified is not changed
327    #[test]
328    fn no_change() {
329        let precision = Q::from((1, u64::MAX - 1));
330
331        assert_eq!(Q::ONE, Q::ONE.simplify(&precision));
332        assert_eq!(Q::MINUS_ONE, Q::MINUS_ONE.simplify(&precision));
333        assert_eq!(Q::ZERO, Q::ZERO.simplify(precision));
334    }
335
336    /// Ensure that no sign change can occurr.
337    #[test]
338    fn no_change_of_sign() {
339        assert!(Q::ZERO < Q::ONE.simplify(2));
340        assert!(Q::ZERO > Q::MINUS_ONE.simplify(2));
341    }
342}
343
344#[cfg(test)]
345mod test_randomized_rounding {
346    use crate::rational::Q;
347    use crate::utils::sample::test_seed;
348
349    /// Ensure that a `r < 0` throws an error
350    #[test]
351    fn negative_r() {
352        for seeded in [false, true] {
353            let value = Q::from((2, 3));
354            assert!(
355                value
356                    .randomized_rounding_optionally_seeded(-1, seeded.then(test_seed))
357                    .is_err()
358            );
359        }
360    }
361
362    /// Checks whether the same seed results in the same sample.
363    #[test]
364    fn same_seed_same_sample() {
365        use crate::rational::Q;
366        let value = Q::from((5, 2));
367        let sample_0 = value.randomized_rounding_seeded(3, [42; 32]).unwrap();
368        let sample_1 = value.randomized_rounding_seeded(3, [42; 32]).unwrap();
369
370        assert_eq!(sample_0, sample_1);
371    }
372
373    /// Checks whether the seeded variant outputs the known sample for a fixed seed.
374    #[test]
375    fn known_answer() {
376        use crate::rational::Q;
377        let value = Q::from((5, 2));
378        let sample = value.randomized_rounding_seeded(3, [42; 32]).unwrap();
379
380        assert_eq!("1", sample.to_string());
381    }
382
383    /// Checks whether the unseeded and the seeded function are available.
384    #[test]
385    fn availability() {
386        use crate::rational::Q;
387        let value = Q::from((5, 2));
388
389        // explicitly checks the existence of the unseeded and the seeded function
390        let _ = value.randomized_rounding(3);
391        let _ = value.randomized_rounding_seeded(3, test_seed());
392    }
393}
qfall_math/rational/q/rounding.rs

qfall_math/rational/q/
rounding.rs
