qfall_math/integer/mat_poly_over_z/sample/
discrete_gauss.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 contains algorithms for sampling according to the discrete Gaussian distribution.
10
11use crate::{
12    error::MathError,
13    integer::{MatPolyOverZ, MatZ, PolyOverZ},
14    rational::{PolyOverQ, Q},
15    traits::{
16        Concatenate, FromCoefficientEmbedding, IntoCoefficientEmbedding, MatrixDimensions,
17        MatrixSetEntry, SetCoefficient,
18    },
19    utils::{
20        index::evaluate_index,
21        sample::discrete_gauss::{DiscreteGaussianIntegerSampler, LookupTableSetting, TAILCUT},
22    },
23};
24use std::fmt::Display;
25
26impl MatPolyOverZ {
27    /// Initializes a new matrix with dimensions `num_rows` x `num_columns` and with each entry
28    /// sampled independently according to the discrete Gaussian distribution,
29    /// using [`PolyOverZ::sample_discrete_gauss`].
30    ///
31    /// Parameters:
32    /// - `num_rows`: specifies the number of rows the new matrix should have
33    /// - `num_cols`: specifies the number of columns the new matrix should have
34    /// - `max_degree`: specifies the included maximal degree the created [`PolyOverZ`] should have
35    /// - `center`: specifies the positions of the center with peak probability
36    /// - `s`: specifies the Gaussian parameter, which is proportional
37    ///   to the standard deviation `sigma * sqrt(2 * pi) = s`
38    ///
39    /// Returns a [`MatPolyOverZ`] with each entry sampled independently from the
40    /// specified discrete Gaussian distribution or an error if `s < 0`.
41    ///
42    /// # Examples
43    /// ```
44    /// use qfall_math::integer::MatPolyOverZ;
45    ///
46    /// let matrix = MatPolyOverZ::sample_discrete_gauss(3, 1, 5, 0, 1.25f32).unwrap();
47    /// ```
48    ///
49    /// # Errors and Failures
50    /// - Returns a [`MathError`] of type [`InvalidIntegerInput`](MathError::InvalidIntegerInput)
51    ///   if `s < 0`.
52    ///
53    /// # Panics ...
54    /// - if the provided number of rows and columns are not suited to create a matrix.
55    ///   For further information see [`MatPolyOverZ::new`].
56    /// - if `max_degree` is negative, or does not fit into an [`i64`].
57    pub fn sample_discrete_gauss(
58        num_rows: impl TryInto<i64> + Display,
59        num_cols: impl TryInto<i64> + Display,
60        max_degree: impl TryInto<i64> + Display,
61        center: impl Into<Q>,
62        s: impl Into<Q>,
63    ) -> Result<MatPolyOverZ, MathError> {
64        let max_degree = evaluate_index(max_degree).unwrap();
65        let mut matrix = MatPolyOverZ::new(num_rows, num_cols);
66
67        let mut dgis = DiscreteGaussianIntegerSampler::init(
68            center,
69            s,
70            unsafe { TAILCUT },
71            LookupTableSetting::FillOnTheFly,
72        )?;
73
74        for row in 0..matrix.get_num_rows() {
75            for col in 0..matrix.get_num_columns() {
76                let mut entry = PolyOverZ::default();
77                for index in 0..=max_degree {
78                    let sample = dgis.sample_z();
79                    unsafe { entry.set_coeff_unchecked(index, &sample) };
80                }
81
82                unsafe { matrix.set_entry_unchecked(row, col, entry) };
83            }
84        }
85
86        Ok(matrix)
87    }
88
89    /// SampleD samples a discrete Gaussian from the lattice with a provided `basis`.
90    ///
91    /// We do not check whether `basis` is actually a basis. Hence, the callee is
92    /// responsible for making sure that `basis` provides a suitable basis.
93    ///
94    /// Parameters:
95    /// - `basis`: specifies a basis for the lattice from which is sampled
96    /// - `k`: the maximal length the polynomial can have
97    /// - `center`: specifies the positions of the center with peak probability
98    /// - `s`: specifies the Gaussian parameter, which is proportional
99    ///   to the standard deviation `sigma * sqrt(2 * pi) = s`
100    ///
101    /// Returns a vector of polynomials sampled according to the
102    /// discrete Gaussian distribution or an error if the basis is not a row vector,
103    /// `s < 0`, or the number of rows of the `basis` and `center` differ.
104    ///
105    /// # Example
106    /// ```
107    /// use qfall_math::{
108    ///     integer::MatPolyOverZ,
109    ///     rational::PolyOverQ,
110    /// };
111    /// use std::str::FromStr;
112    ///
113    /// let basis = MatPolyOverZ::from_str("[[1  1, 3  0 1 -1, 2  2 2]]").unwrap();
114    /// let center = vec![PolyOverQ::default()];
115    ///
116    /// let sample = MatPolyOverZ::sample_d(&basis, 3, &center, 10.5_f64).unwrap();
117    /// ```
118    ///
119    /// # Errors and Failures
120    /// - Returns a [`MathError`] of type [`VectorFunctionCalledOnNonVector`](MathError::VectorFunctionCalledOnNonVector),
121    ///   if the basis is not a row vector.
122    /// - Returns a [`MathError`] of type [`InvalidIntegerInput`](MathError::InvalidIntegerInput)
123    ///   if `s < 0`.
124    /// - Returns a [`MathError`] of type [`MismatchingMatrixDimension`](MathError::MismatchingMatrixDimension)
125    ///   if the number of rows of the `basis` and `center` differ.
126    ///
127    /// This function implements SampleD according to:
128    /// - \[1\] Gentry, Craig and Peikert, Chris and Vaikuntanathan, Vinod (2008).
129    ///   Trapdoors for hard lattices and new cryptographic constructions.
130    ///   In: Proceedings of the fortieth annual ACM symposium on Theory of computing.
131    ///   <https://dl.acm.org/doi/pdf/10.1145/1374376.1374407>
132    ///
133    /// # Panics ...
134    /// - if the polynomials have higher length than the provided upper bound `k`
135    pub fn sample_d(
136        basis: &Self,
137        k: impl Into<i64>,
138        center: &[PolyOverQ],
139        s: impl Into<Q>,
140    ) -> Result<MatPolyOverZ, MathError> {
141        let k = k.into();
142        // use coefficient embedding and then call sampleD for the matrix representation
143        let basis_embedded = basis.into_coefficient_embedding(k);
144
145        // use coefficient embedding to get center
146        let mut center_embedded = center[0].into_coefficient_embedding(k);
147        for row in center.iter().skip(1) {
148            let c_row = row.into_coefficient_embedding(k);
149            center_embedded = center_embedded.concat_vertical(&c_row)?;
150        }
151
152        let sample = MatZ::sample_d(&basis_embedded, &center_embedded, s)?;
153
154        Ok(MatPolyOverZ::from_coefficient_embedding((&sample, k - 1)))
155    }
156}
157
158#[cfg(test)]
159mod test_sample_discrete_gauss {
160    use crate::{integer::MatPolyOverZ, traits::MatrixGetEntry};
161
162    /// Checks whether `sample_discrete_gauss` is available for all types
163    /// implementing [`Into<Z>`], i.e. u8, u16, u32, u64, i8, ...
164    /// or [`Into<Q>`], i.e. u8, i16, f32, Z, Q, ...
165    #[test]
166    fn availability() {
167        let _ = MatPolyOverZ::sample_discrete_gauss(2, 3, 128, 1u16, 2);
168        let _ = MatPolyOverZ::sample_discrete_gauss(1, 3, 128, 1u8, 2.0);
169        let _ = MatPolyOverZ::sample_discrete_gauss(2_i64, 3, 128, 1u32, 2.0_f64);
170        let _ = MatPolyOverZ::sample_discrete_gauss(2_i32, 3, 128, 1u64, 1);
171        let _ = MatPolyOverZ::sample_discrete_gauss(2_i16, 3, 128, 1i64, 3_u64);
172        let _ = MatPolyOverZ::sample_discrete_gauss(2_i8, 3, 128, 1i32, 1);
173        let _ = MatPolyOverZ::sample_discrete_gauss(2_u64, 3, 128, 1i16, 1);
174        let _ = MatPolyOverZ::sample_discrete_gauss(2_u32, 3, 128, 1i8, 1);
175        let _ = MatPolyOverZ::sample_discrete_gauss(2_u16, 3, 128, 1i64, 2);
176        let _ = MatPolyOverZ::sample_discrete_gauss(2_u8, 3, 128, -2, 3);
177        let _ = MatPolyOverZ::sample_discrete_gauss(1, 3, 128, 4, 3);
178        let _ = MatPolyOverZ::sample_discrete_gauss(3, 3, 128, 1.25f64, 3);
179        let _ = MatPolyOverZ::sample_discrete_gauss(4, 3, 128, 15.75f32, 3);
180    }
181
182    /// Ensures that the resulting entries have correct degree.
183    #[test]
184    fn correct_degree_entries() {
185        let degrees = [1, 3, 7, 15, 32, 120];
186        for degree in degrees {
187            let res = MatPolyOverZ::sample_discrete_gauss(1, 1, degree, i64::MAX, 1).unwrap();
188
189            assert_eq!(
190                res.get_entry(0, 0).unwrap().get_degree(),
191                degree,
192                "Could fail with negligible probability."
193            );
194        }
195    }
196
197    /// Checks whether the maximum degree needs to be at least 0.
198    #[test]
199    #[should_panic]
200    fn invalid_max_degree() {
201        let _ = MatPolyOverZ::sample_discrete_gauss(2, 2, -1, 0, 1).unwrap();
202    }
203}
204
205#[cfg(test)]
206mod test_sample_d {
207    use crate::{
208        integer::{MatPolyOverZ, MatZ, Z},
209        rational::{PolyOverQ, Q},
210        traits::IntoCoefficientEmbedding,
211    };
212    use std::str::FromStr;
213
214    /// Ensure that the sample is from the base.
215    #[test]
216    fn ensure_sampled_from_base() {
217        let base = MatPolyOverZ::from_str("[[1  1, 3  0 1 -1]]").unwrap();
218        let center = vec![PolyOverQ::default()];
219
220        for _ in 0..10 {
221            let sample = MatPolyOverZ::sample_d(&base, 3, &center, 10.5_f64).unwrap();
222            let sample_vec = sample.into_coefficient_embedding(3);
223            let orthogonal = MatZ::from_str("[[0],[1],[1]]").unwrap();
224
225            assert_eq!(Z::ZERO, sample_vec.dot_product(&orthogonal).unwrap());
226        }
227    }
228
229    /// Ensure that the sample is from the base for higher dimensional bases.
230    #[test]
231    fn ensure_sampled_from_base_higher_dimension() {
232        let base = MatPolyOverZ::from_str("[[1  1, 3  0 1 -1],[3  0 1 -1, 1  1],[0, 0]]").unwrap();
233        let center = vec![
234            PolyOverQ::default(),
235            PolyOverQ::default(),
236            PolyOverQ::default(),
237        ];
238
239        let orthogonal = MatZ::from_str("[[0, 1, 1, 0, 1, 1, 0 , 0, 0]]").unwrap();
240        for _ in 0..10 {
241            let sample = MatPolyOverZ::sample_d(&base, 3, &center, 10.5_f64).unwrap();
242            let sample_embedded = sample.into_coefficient_embedding(3);
243
244            assert_eq!(MatZ::new(1, 1), &orthogonal * &sample_embedded);
245        }
246    }
247
248    /// Checks whether `sample_d` is available for all types
249    /// implementing [`Into<Z>`], i.e. u8, u16, u32, u64, i8, ...
250    /// or [`Into<Q>`], i.e. u8, i16, f32, Z, Q, ...
251    #[test]
252    fn availability() {
253        let basis = MatPolyOverZ::from_str("[[1  1, 2  0 1, 3  0 0 1]]").unwrap();
254        let center = vec![PolyOverQ::default()];
255        let n = Z::from(1024);
256        let s = Q::ONE;
257
258        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1u16);
259        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1u8);
260        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1u32);
261        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1u64);
262        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1i64);
263        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1i32);
264        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1i16);
265        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1i8);
266        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1i64);
267        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, &n);
268        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, &s);
269        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 1.25f64);
270        let _ = MatPolyOverZ::sample_d(&basis, 3, &center, 15.75f32);
271    }
272}