qfall_math/rational/poly_over_q/sample/
gauss.rs

1// Copyright © 2024 Marcel Luca Schmidt
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 contains algorithms for sampling according
10//! to the Gaussian distribution.
11
12use crate::{
13    error::MathError,
14    rational::{PolyOverQ, Q},
15    traits::SetCoefficient,
16    utils::index::evaluate_index,
17};
18use probability::{
19    prelude::{Gaussian, Sample},
20    source,
21};
22use rand::RngCore;
23use std::fmt::Display;
24
25impl PolyOverQ {
26    /// Initializes a new [`PolyOverQ`] with maximum degree `max_degree`
27    /// and with each coefficient sampled independently according to the
28    /// Gaussian distribution, using [`Q::sample_gauss`].
29    ///
30    /// Parameters:
31    /// - `max_degree`: specifies the included maximal degree the created [`PolyOverQ`] should have
32    /// - `center`: specifies the center for each coefficient of the polynomial
33    /// - `sigma`: specifies the standard deviation
34    ///
35    /// Returns a fresh [`PolyOverQ`] instance of maximum degree `max_degree`
36    /// with coefficients chosen independently according the Gaussian
37    /// distribution or a [`MathError`] if the specified parameters were not chosen
38    /// appropriately (`sigma > 0`).
39    ///
40    /// # Examples
41    /// ```
42    /// use qfall_math::rational::PolyOverQ;
43    ///
44    /// let sample = PolyOverQ::sample_gauss(2, 0, 1).unwrap();
45    /// ```
46    ///
47    /// # Errors and Failures
48    /// - Returns a [`MathError`] of type [`NonPositive`](MathError::NonPositive)
49    ///   if `sigma <= 0`.
50    ///
51    /// # Panics ...
52    /// - if `max_degree` is negative, or does not fit into an [`i64`].
53    pub fn sample_gauss(
54        max_degree: impl TryInto<i64> + Display,
55        center: impl Into<Q>,
56        sigma: impl Into<f64>,
57    ) -> Result<Self, MathError> {
58        let max_degree = evaluate_index(max_degree).unwrap();
59
60        let center = center.into();
61        let sigma: f64 = sigma.into();
62        let mut poly = PolyOverQ::default();
63
64        if sigma <= 0.0 {
65            return Err(MathError::NonPositive(format!(
66                "The sigma has to be positive and not zero, but the provided value is {sigma}."
67            )));
68        }
69        let mut rng = rand::rng();
70        let mut source = source::default(rng.next_u64());
71
72        // Instead of sampling with a center of c, we sample with center 0 and add the
73        // center later. These are equivalent and this way we can sample in larger ranges
74        let sampler = Gaussian::new(0.0, sigma);
75
76        for index in 0..=max_degree {
77            let mut sample = Q::from(sampler.sample(&mut source));
78            sample += &center;
79            unsafe { poly.set_coeff_unchecked(index, &sample) };
80        }
81        Ok(poly)
82    }
83}
84
85#[cfg(test)]
86mod test_sample_gauss {
87    use crate::{
88        integer::Z,
89        rational::{PolyOverQ, Q},
90    };
91
92    /// Checks whether `sample_gauss` is available for all types
93    /// implementing [`Into<Z>`], i.e. u8, u16, u32, u64, i8, ...
94    /// or [`Into<Q>`], i.e. u8, i16, f32, Z, Q, ...
95    #[test]
96    fn availability() {
97        let center_q = Q::ZERO;
98        let center_z = Z::ZERO;
99
100        let _ = PolyOverQ::sample_gauss(1u8, 16u8, 0f32);
101        let _ = PolyOverQ::sample_gauss(1u16, 16u16, 0f64);
102        let _ = PolyOverQ::sample_gauss(1u32, 16u32, 0f32);
103        let _ = PolyOverQ::sample_gauss(1u64, 16u64, 0f64);
104        let _ = PolyOverQ::sample_gauss(1i8, 16u8, 0f32);
105        let _ = PolyOverQ::sample_gauss(1i16, 16i16, 0f32);
106        let _ = PolyOverQ::sample_gauss(1i32, 16i32, 0f32);
107        let _ = PolyOverQ::sample_gauss(1i64, 16i64, 0f64);
108        let _ = PolyOverQ::sample_gauss(1u8, center_q, 0f32);
109        let _ = PolyOverQ::sample_gauss(1u8, center_z, 0f64);
110    }
111
112    /// Checks whether the number of coefficients is correct.
113    #[test]
114    fn nr_coeffs() {
115        let degrees = [1, 3, 7, 15, 32, 120];
116        for degree in degrees {
117            let res = PolyOverQ::sample_gauss(degree, i64::MAX, 1).unwrap();
118
119            assert_eq!(
120                res.get_degree(),
121                degree,
122                "Could fail with negligible probability."
123            );
124        }
125    }
126
127    /// Checks is as error is returned, if sigma is less than 1.
128    #[test]
129    fn invalid_sigma() {
130        assert!(PolyOverQ::sample_gauss(1, 0, 0).is_err());
131        assert!(PolyOverQ::sample_gauss(1, 0, -1).is_err());
132    }
133
134    /// Checks whether the maximum degree needs to be at least 0.
135    #[test]
136    #[should_panic]
137    fn invalid_max_degree() {
138        let _ = PolyOverQ::sample_gauss(-1, 0, 1).unwrap();
139    }
140}
qfall_math/rational/poly_over_q/sample/gauss.rs

qfall_math/rational/poly_over_q/sample/
gauss.rs
