qfall_math/utils/sample/discrete_gauss.rs
1// Copyright 2025 Jan Niklas Siemer
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 core functionality to sample according to the
10//! discrete gaussian distribution.
11//!
12//! The main references are listed in the following
13//! and will be further referenced in submodules by these numbers:
14//! - \[1\] Gentry, Craig and Peikert, Chris and Vaikuntanathan, Vinod (2008).
15//! Trapdoors for hard lattices and new cryptographic constructions.
16//! In: Proceedings of the fortieth annual ACM symposium on Theory of computing.
17//! <https://citeseerx.ist.psu.edu/document?doi=d9f54077d568784c786f7b1d030b00493eb3ae35>
18
19use super::uniform::{SamplerRng, UniformIntegerSampler};
20use crate::{
21 error::{MathError, StringConversionError},
22 integer::{MatZ, Z},
23 rational::{MatQ, Q},
24 traits::{MatrixDimensions, MatrixGetSubmatrix, Pow},
25};
26use rand::RngExt;
27use serde::{Deserialize, Serialize};
28use std::collections::HashMap;
29
30/// Defines whether a lookup-table should be precomputed, filled on-the-fly,
31/// or not used at all for a discrete Gaussian sampler.
32#[derive(PartialEq, Clone, Copy, Serialize, Deserialize, Debug)]
33pub enum LookupTableSetting {
34 Precompute,
35 FillOnTheFly,
36 NoLookup,
37}
38
39/// This is the global variable used in all `sample_discrete_gauss` and `sample_d`
40/// functions. Its value should be in `ω(log(sqrt(n)))`. We set it (as most other libraries)
41/// statically to `6.0`.
42///
43/// You can use and change in an `unsafe` environment.
44/// ```compile_fail
45/// unsafe { TAILCUT = 4.0 };
46/// ```
47/// Make sure that the tailcut stays positive and large enough for your purposes.
48/// If you use multi-threading, read up on the behaviour of a `static mut` variable.
49/// Our tests only cover cases where `TAILCUT = 6.0`.
50pub static mut TAILCUT: f64 = 6.0;
51
52/// Enables for discrete Gaussian sampling out of
53/// `[⌈center - s * tailcut⌉ , ⌊center + s * tailcut⌋ ]`.
54///
55/// **WARNING:** If the attributes are not set using [`DiscreteGaussianIntegerSampler::init`],
56/// we can't guarantee sampling from the correct discrete Gaussian distribution.
57/// Altering any value will invalidate the [`HashMap`] in `table` and might invalidate
58/// other attributes, too.
59///
60/// Attributes:
61/// - `center`: specifies the position of the center with peak probability
62/// - `s`: specifies the Gaussian parameter, which is proportional
63/// to the standard deviation `sigma * sqrt(2 * pi) = s`
64/// - `lower_bound`: specifies the lower bound to sample uniformly from
65/// - `interval_size`: specifies the interval size to sample uniformly from
66/// - `lookup_table_setting`: Specifies whether a lookup-table should be used and
67/// how it should be filled, i.e. lazily on-the-fly (impacting sampling time slightly) or precomputed
68/// - `table`: the lookup-table if one is used
69/// - `rng`: defines the [`ThreadRng`](rand::rngs::ThreadRng) or seeded
70/// [`StdRng`](rand::rngs::StdRng) that's used for sampling.
71///
72/// # Examples
73/// ```
74/// use qfall_math::utils::sample::discrete_gauss::{DiscreteGaussianIntegerSampler, LookupTableSetting};
75/// let center = 0.0;
76/// let gaussian_parameter = 1.0;
77/// let tailcut = 6.0;
78///
79/// let mut dgis = DiscreteGaussianIntegerSampler::init(center, gaussian_parameter, tailcut, LookupTableSetting::NoLookup, None).unwrap();
80///
81/// let sample = dgis.sample_z();
82/// ```
83#[derive(Debug, Serialize, Deserialize, Clone)]
84pub struct DiscreteGaussianIntegerSampler {
85 pub center: Q,
86 pub s: Q,
87 pub lower_bound: Z,
88 pub interval_size: Z,
89 pub lookup_table_setting: LookupTableSetting,
90 #[serde(with = "table_as_pairs")]
91 pub table: HashMap<Z, f64>,
92 #[serde(skip)]
93 pub rng: SamplerRng,
94}
95
96impl DiscreteGaussianIntegerSampler {
97 /// Initializes the [`DiscreteGaussianIntegerSampler`] with
98 /// - `center` as the center of the discrete Gaussian to sample from,
99 /// - `s` defining the Gaussian parameter, which is proportional
100 /// to the standard deviation `sigma * sqrt(2 * pi) = s`,
101 /// - `lower_bound` as `⌈center - 6 * s⌉`,
102 /// - `interval_size` as `⌊center + 6 * s⌋ - ⌈center - 6 * s⌉ + 1`, and
103 /// - `table` as an empty [`HashMap`] to store evaluations of the Gaussian function, and
104 /// - `rng` as a [`StdRng`](rand::rngs::StdRng) seeded with `seed` if `seed` is provided,
105 /// and as a fresh [`ThreadRng`](rand::rngs::ThreadRng) otherwise.
106 ///
107 /// Parameters:
108 /// - `center`: as the center of the discrete Gaussian to sample from
109 /// - `s`: specifies the Gaussian parameter, which is proportional
110 /// to the standard deviation `sigma * sqrt(2 * pi) = s`
111 /// - `tailcut`: specifies the number of Gaussian parameters `s` the interval
112 /// to sample from reaches from `center` to both sides
113 /// - `lookup_table_setting`: specifies whether a lookup-table should be used and
114 /// how it should be filled
115 /// - `seed`: specifies an optional 256-bit seed for the internal [`StdRng`](rand::rngs::StdRng).
116 /// If `None` is provided, a fresh [`ThreadRng`](rand::rngs::ThreadRng) is used instead.
117 ///
118 /// Returns a [`DiscreteGaussianIntegerSampler`] or a [`MathError`]
119 /// if the specified parameters were not chosen appropriately,
120 /// i.e. `tailcut < 0` or `s < 0`.
121 ///
122 /// # Examples
123 /// ```
124 /// use qfall_math::utils::sample::discrete_gauss::{DiscreteGaussianIntegerSampler, LookupTableSetting};
125 /// let center = 0.0;
126 /// let gaussian_parameter = 1.0;
127 /// let tailcut = 6.0;
128 ///
129 /// let dgis = DiscreteGaussianIntegerSampler::init(center, gaussian_parameter, tailcut, LookupTableSetting::Precompute, None).unwrap();
130 ///
131 /// let mut dgis_seeded_0 = DiscreteGaussianIntegerSampler::init(center, gaussian_parameter, tailcut, LookupTableSetting::Precompute, Some([42; 32])).unwrap();
132 /// let mut dgis_seeded_1 = DiscreteGaussianIntegerSampler::init(center, gaussian_parameter, tailcut, LookupTableSetting::Precompute, Some([42; 32])).unwrap();
133 /// assert_eq!(dgis_seeded_0.sample_z(), dgis_seeded_1.sample_z());
134 /// ```
135 ///
136 /// # Errors and Failures
137 /// - Returns a [`MathError`] of type [`InvalidIntegerInput`](MathError::InvalidIntegerInput)
138 /// if `tailcut < 0` or `s < 0`.
139 pub fn init(
140 center: impl Into<Q>,
141 s: impl Into<Q>,
142 tailcut: impl Into<Q>,
143 lookup_table_setting: LookupTableSetting,
144 seed: Option<[u8; 32]>,
145 ) -> Result<Self, MathError> {
146 let center = center.into();
147 let mut s = s.into();
148 let tailcut = tailcut.into();
149 if tailcut < Q::ZERO {
150 return Err(MathError::InvalidIntegerInput(format!(
151 "The value {tailcut} was provided for parameter tailcut of the function sample_z.
152 This function expects this input no smaller than 0."
153 )));
154 }
155 if s < Q::ZERO {
156 return Err(MathError::InvalidIntegerInput(format!(
157 "The value {s} was provided for parameter s of the function sample_z.
158 This function expects this input to be no smaller than 0."
159 )));
160 }
161 if s == Q::ZERO {
162 // Ensure that s != 0 s.t. we can divide by s^2 in the Gaussian function
163 s = Q::from(0.00001);
164 }
165
166 let lower_bound = (¢er - &s * &tailcut).ceil();
167 let upper_bound = (¢er + &s * tailcut).floor();
168 // interval [lower_bound, upper_bound] has upper_bound - lower_bound + 1 elements in it
169 let interval_size = &upper_bound - &lower_bound + Z::ONE;
170
171 let mut table = HashMap::new();
172
173 if lookup_table_setting != LookupTableSetting::NoLookup && interval_size > u16::MAX {
174 println!(
175 "WARNING: A completely filled lookup table will exceed 2^16 entries. You should reconsider your sampling method for discrete Gaussians."
176 )
177 }
178
179 if lookup_table_setting == LookupTableSetting::Precompute {
180 let mut i = lower_bound.clone();
181 while i <= upper_bound {
182 let evaluated_gauss_function = gaussian_function(&i, ¢er, &s);
183 table.insert(i.clone(), evaluated_gauss_function);
184 i += Z::ONE;
185 }
186 }
187
188 Ok(Self {
189 center,
190 s,
191 lower_bound,
192 interval_size,
193 lookup_table_setting,
194 table,
195 rng: SamplerRng::new(seed),
196 })
197 }
198
199 /// Chooses a sample according to the discrete Gaussian distribution out of
200 /// `[lower_bound , lower_bound + interval_size ]`.
201 ///
202 /// This function implements discrete Gaussian sampling according to the definition of
203 /// SampleZ as in [\[1\]](<index.html#:~:text=[1]>).
204 ///
205 /// # Examples
206 /// ```
207 /// use qfall_math::utils::sample::discrete_gauss::{DiscreteGaussianIntegerSampler, LookupTableSetting};
208 /// let center = 0.0;
209 /// let gaussian_parameter = 1.0;
210 /// let tailcut = 6.0;
211 ///
212 /// let mut dgis = DiscreteGaussianIntegerSampler::init(center, gaussian_parameter, tailcut, LookupTableSetting::Precompute, None).unwrap();
213 ///
214 /// let sample = dgis.sample_z();
215 /// ```
216 pub fn sample_z(&mut self) -> Z {
217 // the clone of `self.rng` shares its state with `self.rng`
218 let mut uis =
219 UniformIntegerSampler::init_with_rng(&self.interval_size, self.rng.clone()).unwrap();
220 loop {
221 // sample x in [c - s * tailcut, c + s * tailcut]
222 let sample = &self.lower_bound + uis.sample();
223
224 let evaluated_gauss_function: &f64 = match self.lookup_table_setting {
225 LookupTableSetting::NoLookup => &gaussian_function(&sample, &self.center, &self.s),
226 LookupTableSetting::FillOnTheFly => {
227 let pot_evaluated_gauss_function = self.table.get(&sample);
228 match pot_evaluated_gauss_function {
229 Some(x) => x,
230 None => &{
231 // if the entry doesn't exist yet, compute and insert it
232 let evaluated_function =
233 gaussian_function(&sample, &self.center, &self.s);
234 self.table.insert(sample.clone(), evaluated_function);
235 evaluated_function
236 },
237 }
238 }
239 LookupTableSetting::Precompute => self.table.get(&sample).unwrap(),
240 };
241
242 let random_f64: f64 = self.rng.random();
243 if evaluated_gauss_function >= &random_f64 {
244 return sample;
245 }
246 }
247 }
248}
249
250/// Serializes and deserializes the lookup table of a [`DiscreteGaussianIntegerSampler`]
251/// as a sequence of pairs, as [`Z`] can not be used as a key of a map in every format, e.g. JSON.
252mod table_as_pairs {
253 use crate::integer::Z;
254 use serde::{Deserialize, Deserializer, Serializer};
255 use std::collections::HashMap;
256
257 /// Serializes `table` as a sequence of `(key, value)` pairs.
258 pub(super) fn serialize<S: Serializer>(
259 table: &HashMap<Z, f64>,
260 serializer: S,
261 ) -> Result<S::Ok, S::Error> {
262 serializer.collect_seq(table.iter())
263 }
264
265 /// Deserializes a sequence of `(key, value)` pairs into a [`HashMap`].
266 pub(super) fn deserialize<'de, D: Deserializer<'de>>(
267 deserializer: D,
268 ) -> Result<HashMap<Z, f64>, D::Error> {
269 let pairs = Vec::<(Z, f64)>::deserialize(deserializer)?;
270 Ok(pairs.into_iter().collect())
271 }
272}
273
274/// Computes the value of the Gaussian function for `x`.
275///
276/// **Warning**: This functions assumes `s != 0`.
277///
278/// Parameters:
279/// - `x`: specifies the value/ sample for which the Gaussian function's value is computed
280/// - `c`: specifies the position of the center with peak probability
281/// - `s`: specifies the Gaussian parameter, which is proportional
282/// to the standard deviation `sigma * sqrt(2 * pi) = s`
283///
284/// Returns the computed value of the Gaussian function for `x`.
285///
286/// # Examples
287/// ```
288/// use qfall_math::{integer::Z, rational::Q};
289/// use qfall_math::utils::sample::discrete_gauss::gaussian_function;
290/// let sample = Z::ONE;
291/// let center = Q::ZERO;
292/// let gaussian_parameter = Q::ONE;
293///
294/// let probability = gaussian_function(&sample, ¢er, &gaussian_parameter);
295/// ```
296///
297/// # Panics ...
298/// - if `s = 0`.
299/// - if `-π * (x - c)^2 / s^2` is larger than [`f64::MAX`]
300pub fn gaussian_function(x: &Z, c: &Q, s: &Q) -> f64 {
301 let num = Q::MINUS_ONE * Q::PI * (x - c).pow(2).unwrap();
302 let den = s.pow(2).unwrap();
303 let res = f64::from(&(num / den));
304 res.exp()
305}
306
307/// SampleD samples a discrete Gaussian from the lattice with `basis` using [`sample_z`] as a subroutine.
308///
309/// We do not check whether `basis` is actually a basis. Hence, the callee is
310/// responsible for making sure that `basis` provides a suitable basis.
311///
312/// Parameters:
313/// - `basis`: specifies a basis for the lattice from which is sampled
314/// - `n`: specifies the range from which [`sample_z`] samples
315/// - `center`: specifies the positions of the center with peak probability
316/// - `s`: specifies the Gaussian parameter, which is proportional
317/// to the standard deviation `sigma * sqrt(2 * pi) = s`
318/// - `seed`: specifies an optional 256-bit seed for the PRNG used during sampling.
319/// If `None` is provided, a fresh [`ThreadRng`](rand::rngs::ThreadRng) is used instead.
320///
321/// Returns a vector with discrete gaussian error based on a lattice point
322/// as in [\[1\]](<index.html#:~:text=[1]>): SampleD or a [`MathError`], if the
323/// `n <= 1` or `s <= 0`, the number of rows of the `basis` and `center` differ,
324/// or `center` is not a column vector.
325///
326/// # Examples
327/// ```compile_fail
328/// use qfall_math::{integer::{MatZ, Z}, rational::{MatQ, Q}};
329/// use qfall_math::utils::sample::discrete_gauss::sample_d;
330/// let basis = MatZ::identity(5, 5);
331/// let n = Z::from(1024);
332/// let center = MatQ::new(5, 1);
333/// let gaussian_parameter = Q::ONE;
334///
335/// let sample = sample_d(basis, &n, ¢er, &gaussian_parameter, None).unwrap();
336/// ```
337///
338/// # Errors and Failures
339/// - Returns a [`MathError`] of type [`InvalidIntegerInput`](MathError::InvalidIntegerInput)
340/// if `n <= 1` or `s <= 0`.
341/// - Returns a [`MathError`] of type [`MismatchingMatrixDimension`](MathError::MismatchingMatrixDimension)
342/// if the number of rows of the `basis` and `center` differ.
343/// - Returns a [`MathError`] of type [`StringConversionError`](MathError::StringConversionError)
344/// if `center` is not a column vector.
345pub(crate) fn sample_d(
346 basis: &MatZ,
347 center: &MatQ,
348 s: &Q,
349 seed: Option<[u8; 32]>,
350) -> Result<MatZ, MathError> {
351 let basis_gso = MatQ::from(basis).gso();
352 sample_d_precomputed_gso(basis, &basis_gso, center, s, seed)
353}
354
355/// SampleD samples a discrete Gaussian from the lattice with `basis` using [`sample_z`] as a subroutine.
356///
357/// We do not check whether `basis` is actually a basis or whether `basis_gso` is
358/// actually the gso of `basis`. Hence, the callee is responsible for making sure that
359/// `basis` provides a suitable basis and `basis_gso` is a corresponding GSO.
360///
361/// Parameters:
362/// - `basis`: specifies a basis for the lattice from which is sampled
363/// - `basis_gso`: specifies the precomputed gso for `basis`
364/// - `n`: specifies the range from which [`sample_z`] samples
365/// - `center`: specifies the positions of the center with peak probability
366/// - `s`: specifies the Gaussian parameter, which is proportional
367/// to the standard deviation `sigma * sqrt(2 * pi) = s`
368/// - `seed`: specifies an optional 256-bit seed for the PRNG used during sampling.
369/// If `None` is provided, a fresh [`ThreadRng`](rand::rngs::ThreadRng) is used instead.
370///
371/// Returns a vector with discrete gaussian error based on a lattice point
372/// as in [\[1\]](<index.html#:~:text=[1]>): SampleD or a [`MathError`], if the
373/// `n <= 1` or `s <= 0`, the number of rows of the `basis` and `center` differ,
374/// or `center` is not a column vector.
375///
376/// # Examples
377/// ```compile_fail
378/// use qfall_math::{integer::{MatZ, Z}, rational::{MatQ, Q}};
379/// use qfall_math::utils::sample::discrete_gauss::sample_d;
380/// let basis = MatZ::identity(5, 5);
381/// let n = Z::from(1024);
382/// let center = MatQ::new(5, 1);
383/// let gaussian_parameter = Q::ONE;
384///
385/// let basis_gso = basis.gso();
386///
387/// let sample = sample_d(basis, &basis_gso, &n, ¢er, &gaussian_parameter, None).unwrap();
388/// ```
389///
390/// # Errors and Failures
391/// - Returns a [`MathError`] of type [`InvalidIntegerInput`](MathError::InvalidIntegerInput)
392/// if `n <= 1` or `s <= 0`.
393/// - Returns a [`MathError`] of type [`MismatchingMatrixDimension`](MathError::MismatchingMatrixDimension)
394/// if the number of rows of the `basis` and `center` differ.
395/// - Returns a [`MathError`] of type [`StringConversionError`](MathError::StringConversionError)
396/// if `center` is not a column vector.
397///
398/// # Panics ...
399/// - if the number of rows/columns of `basis_gso` and `basis` mismatch.
400pub(crate) fn sample_d_precomputed_gso(
401 basis: &MatZ,
402 basis_gso: &MatQ,
403 center: &MatQ,
404 s: &Q,
405 seed: Option<[u8; 32]>,
406) -> Result<MatZ, MathError> {
407 let mut center = center.clone();
408 assert_eq!(
409 basis.get_num_rows(),
410 basis_gso.get_num_rows(),
411 "The provided gso can not be based on the provided base, \
412 as they do not have the same number of rows."
413 );
414 assert_eq!(
415 basis.get_num_columns(),
416 basis_gso.get_num_columns(),
417 "The provided gso can not be based on the provided base, \
418 as they do not have the same number of columns."
419 );
420 if center.get_num_rows() != basis.get_num_rows() {
421 return Err(MathError::MismatchingMatrixDimension(format!(
422 "sample_d requires center and basis to have the same number of columns, but they were {} and {}.",
423 center.get_num_rows(),
424 basis.get_num_rows()
425 )));
426 }
427 if !center.is_column_vector() {
428 Err(StringConversionError::InvalidMatrix(format!(
429 "sample_d expects center to be a column vector, but it has dimensions {}x{}.",
430 center.get_num_rows(),
431 center.get_num_columns()
432 )))?;
433 }
434 if s < &Q::ZERO {
435 return Err(MathError::InvalidIntegerInput(format!(
436 "The value {s} was provided for parameter s of the function sample_z.
437 This function expects this input to be larger than 0."
438 )));
439 }
440
441 let mut out = MatZ::new(basis_gso.get_num_rows(), 1);
442
443 // derives an independent seed for the sampler of each dimension if a seed is provided,
444 // see the documentation of `seedable_function` in `src/macros/seeded.rs`
445 let mut rng = SamplerRng::new(seed);
446
447 for i in (0..basis_gso.get_num_columns()).rev() {
448 // basisvector_i = b_tilde[i]
449 let basisvector_orth_i = unsafe { basis_gso.get_column_unchecked(i) };
450
451 // define the center for sample_z as c_2 = <c, b_tilde[i]> / <b_tilde[i], b_tilde[i]>;
452 let c_2 = center.dot_product(&basisvector_orth_i).unwrap()
453 / basisvector_orth_i.dot_product(&basisvector_orth_i).unwrap();
454
455 // Defines the gaussian parameter to be normalized along the basis vector: s2 = s / ||b_tilde[i]||
456 let s_2 = s / (basisvector_orth_i.norm_eucl_sqrd().unwrap().sqrt());
457
458 // sample z ~ D_{Z, s2, c2}
459 let mut dgis = DiscreteGaussianIntegerSampler::init(
460 &c_2,
461 &s_2,
462 unsafe { TAILCUT },
463 LookupTableSetting::FillOnTheFly,
464 rng.derive_seed(),
465 )?;
466 let z = dgis.sample_z();
467
468 // update the center c = c - z * b[i]
469 let basisvector_i = unsafe { basis.get_column_unchecked(i) };
470 center -= MatQ::from(&(&z * &basisvector_i));
471
472 // out = out + z * b[i]
473 out = &out + &z * &basisvector_i;
474 }
475
476 Ok(out)
477}
478
479#[cfg(test)]
480mod test_discrete_gaussian_integer_sampler {
481 use super::DiscreteGaussianIntegerSampler;
482 use crate::utils::sample::test_seed;
483 use crate::{
484 rational::Q,
485 utils::sample::discrete_gauss::{LookupTableSetting, TAILCUT},
486 };
487
488 /// Checks whether samples are kept in correct interval for a small interval.
489 #[test]
490 fn small_interval() {
491 for seeded in [false, true] {
492 let center = Q::from(15);
493 let gaussian_parameter = Q::from((1, 2));
494
495 let mut dgis = DiscreteGaussianIntegerSampler::init(
496 ¢er,
497 &gaussian_parameter,
498 8.0,
499 LookupTableSetting::FillOnTheFly,
500 seeded.then(test_seed),
501 )
502 .unwrap();
503
504 for _ in 0..64 {
505 let sample = dgis.sample_z();
506
507 assert!(10 <= sample);
508 assert!(sample <= 20);
509 }
510 }
511 }
512
513 /// Checks whether samples are kept in correct interval for a large interval.
514 #[test]
515 fn large_interval() {
516 for seeded in [false, true] {
517 let center = Q::MINUS_ONE;
518 let gaussian_parameter = Q::ONE;
519
520 let mut dgis = DiscreteGaussianIntegerSampler::init(
521 ¢er,
522 &gaussian_parameter,
523 unsafe { TAILCUT },
524 LookupTableSetting::FillOnTheFly,
525 seeded.then(test_seed),
526 )
527 .unwrap();
528
529 for _ in 0..256 {
530 let sample = dgis.sample_z();
531
532 assert!(-64 <= sample);
533 assert!(sample <= 62);
534 }
535 }
536 }
537
538 /// Checks whether `sample_z` returns an error if the gaussian parameter `s < 0`.
539 #[test]
540 fn invalid_gaussian_parameter() {
541 for seeded in [false, true] {
542 let center = Q::ZERO;
543
544 assert!(
545 DiscreteGaussianIntegerSampler::init(
546 ¢er,
547 &Q::MINUS_ONE,
548 6.0,
549 LookupTableSetting::FillOnTheFly,
550 seeded.then(test_seed)
551 )
552 .is_err()
553 );
554 assert!(
555 DiscreteGaussianIntegerSampler::init(
556 ¢er,
557 Q::from(i64::MIN),
558 6.0,
559 LookupTableSetting::FillOnTheFly,
560 seeded.then(test_seed)
561 )
562 .is_err()
563 );
564 }
565 }
566
567 /// Checks whether `sample_z` returns an error if `n < 0`.
568 #[test]
569 fn invalid_tailcut() {
570 for seeded in [false, true] {
571 let center = Q::MINUS_ONE;
572 let gaussian_parameter = Q::ONE;
573
574 assert!(
575 DiscreteGaussianIntegerSampler::init(
576 ¢er,
577 &gaussian_parameter,
578 -0.1,
579 LookupTableSetting::FillOnTheFly,
580 seeded.then(test_seed)
581 )
582 .is_err()
583 );
584 assert!(
585 DiscreteGaussianIntegerSampler::init(
586 ¢er,
587 &gaussian_parameter,
588 i64::MIN,
589 LookupTableSetting::FillOnTheFly,
590 seeded.then(test_seed)
591 )
592 .is_err()
593 );
594 }
595 }
596
597 /// Checks whether a serialized and deserialized sampler keeps its lookup table
598 /// and still samples from the same interval.
599 #[test]
600 fn serialize_deserialize() {
601 let dgis = DiscreteGaussianIntegerSampler::init(
602 Q::from(15),
603 Q::from((1, 2)),
604 8.0,
605 LookupTableSetting::Precompute,
606 Some([42; 32]),
607 )
608 .unwrap();
609
610 let string = serde_json::to_string(&dgis).unwrap();
611 let dgis_deserialized: DiscreteGaussianIntegerSampler =
612 serde_json::from_str(&string).unwrap();
613
614 assert!(!dgis.table.is_empty());
615 assert_eq!(dgis.table, dgis_deserialized.table);
616 let mut dgis = dgis_deserialized;
617
618 for _ in 0..64 {
619 let sample = dgis.sample_z();
620
621 assert!(10 <= sample);
622 assert!(sample <= 20);
623 }
624 }
625
626 /// Checks whether two samplers with the same seed output the same samples
627 /// for every [`LookupTableSetting`].
628 #[test]
629 fn same_seed_same_samples() {
630 let settings = [
631 LookupTableSetting::NoLookup,
632 LookupTableSetting::FillOnTheFly,
633 LookupTableSetting::Precompute,
634 ];
635
636 for setting in settings {
637 let mut dgis_0 =
638 DiscreteGaussianIntegerSampler::init(Q::MINUS_ONE, 4, 6, setting, Some([42; 32]))
639 .unwrap();
640 let mut dgis_1 =
641 DiscreteGaussianIntegerSampler::init(Q::MINUS_ONE, 4, 6, setting, Some([42; 32]))
642 .unwrap();
643
644 for _ in 0..u8::MAX {
645 assert_eq!(dgis_0.sample_z(), dgis_1.sample_z());
646 }
647 }
648 }
649
650 /// Checks whether a clone of a seeded sampler continues the stream of its original,
651 /// i.e. alternating between both yields the same samples as a single seeded sampler.
652 #[test]
653 fn clone_shares_state() {
654 let mut dgis = DiscreteGaussianIntegerSampler::init(
655 0,
656 16,
657 6,
658 LookupTableSetting::NoLookup,
659 Some([42; 32]),
660 )
661 .unwrap();
662 let mut dgis_clone = dgis.clone();
663 let mut dgis_cmp = DiscreteGaussianIntegerSampler::init(
664 0,
665 16,
666 6,
667 LookupTableSetting::NoLookup,
668 Some([42; 32]),
669 )
670 .unwrap();
671
672 for _ in 0..16 {
673 assert_eq!(dgis_cmp.sample_z(), dgis.sample_z());
674 assert_eq!(dgis_cmp.sample_z(), dgis_clone.sample_z());
675 }
676 }
677
678 /// Checks whether the seeded sampler outputs the known samples for a fixed seed.
679 #[test]
680 fn known_answer() {
681 let mut dgis = DiscreteGaussianIntegerSampler::init(
682 Q::MINUS_ONE,
683 16,
684 6,
685 LookupTableSetting::FillOnTheFly,
686 Some([42; 32]),
687 )
688 .unwrap();
689
690 let samples: Vec<String> = (0..8).map(|_| dgis.sample_z().to_string()).collect();
691
692 assert_eq!("3, -3, 5, 3, -12, 13, -1, 0", samples.join(", "));
693 }
694}
695
696#[cfg(test)]
697mod test_gaussian_function {
698 use super::{Q, Z, gaussian_function};
699 use crate::traits::Distance;
700
701 /// Ensures that the doc test would run properly.
702 #[test]
703 fn doc_test() {
704 let sample = Z::ONE;
705 let center = Q::ZERO;
706 let gaussian_parameter = Q::ONE;
707 // result roughly 0.0432139 computed via WolframAlpha
708 let cmp = Q::from((43214, 1_000_000));
709
710 let value = gaussian_function(&sample, ¢er, &gaussian_parameter);
711
712 assert!(cmp.distance(&Q::from(value)) < Q::from((1, 1_000_000)));
713 }
714
715 /// Checks whether the values for small values are computed appropriately
716 /// and with appropriate precision.
717 #[test]
718 fn small_values() {
719 let sample_0 = Z::ZERO;
720 let sample_1 = Z::MINUS_ONE;
721 let center = Q::MINUS_ONE;
722 let gaussian_parameter_0 = Q::from((1, 2));
723 let gaussian_parameter_1 = Q::from((3, 2));
724 // result roughly 0.00000348734 computed via WolframAlpha
725 let cmp_0 = Q::from((349, 100_000_000));
726 // result 0.247520 computed via WolframAlpha
727 let cmp_1 = Q::from((24752, 100_000));
728
729 let res_0 = gaussian_function(&sample_0, ¢er, &gaussian_parameter_0);
730 let res_1 = gaussian_function(&sample_0, ¢er, &gaussian_parameter_1);
731 let res_2 = gaussian_function(&sample_1, ¢er, &gaussian_parameter_0);
732 let res_3 = gaussian_function(&sample_1, ¢er, &gaussian_parameter_1);
733
734 assert!(cmp_0.distance(&Q::from(res_0)) < Q::from((3, 1_000_000_000)));
735 assert!(cmp_1.distance(&Q::from(res_1)) < Q::from((1, 1_000_000)));
736 assert_eq!(1.0, res_2);
737 assert_eq!(1.0, res_3);
738 }
739
740 /// Checks whether the values for large values are computed appropriately
741 /// and with appropriate precision.
742 #[test]
743 fn large_values() {
744 let sample = Z::from(i64::MAX);
745 let center = Q::from(i64::MAX as u64 + 1);
746 let gaussian_parameter = Q::from((1, 2));
747 // result roughly 0.00000348734 computed via WolframAlpha
748 let cmp = Q::from((349, 100_000_000));
749
750 let res = gaussian_function(&sample, ¢er, &gaussian_parameter);
751
752 assert!(cmp.distance(&Q::from(res)) < Q::from((3, 1_000_000_000)));
753 }
754
755 /// Checks whether `s = 0` results in a panic.
756 #[test]
757 #[should_panic]
758 fn invalid_s() {
759 let sample = Z::from(i64::MAX);
760 let center = Q::from(i64::MAX as u64 + 1);
761 let gaussian_parameter = Q::ZERO;
762
763 let _ = gaussian_function(&sample, ¢er, &gaussian_parameter);
764 }
765}
766
767#[cfg(test)]
768mod test_sample_d {
769 use super::sample_d_precomputed_gso;
770 use crate::traits::{Concatenate, MatrixDimensions, MatrixGetSubmatrix, Pow};
771 use crate::utils::sample::discrete_gauss::sample_d;
772 use crate::utils::sample::test_seed;
773 use crate::{
774 integer::{MatZ, Z},
775 rational::{MatQ, Q},
776 };
777 use flint3_sys::fmpz_mat_hnf;
778 use std::str::FromStr;
779
780 /// Ensures that the doc-test compiles and runs properly.
781 #[test]
782 fn doc_test() {
783 for seeded in [false, true] {
784 let basis = MatZ::identity(5, 5);
785 let center = MatQ::new(5, 1);
786 let gaussian_parameter = Q::ONE;
787 let basis_gso = MatQ::from(&basis).gso();
788
789 let _ = sample_d(&basis, ¢er, &gaussian_parameter, seeded.then(test_seed)).unwrap();
790 let _ = sample_d_precomputed_gso(
791 &basis,
792 &basis_gso,
793 ¢er,
794 &gaussian_parameter,
795 seeded.then(test_seed),
796 )
797 .unwrap();
798 }
799 }
800
801 /// Ensures that `sample_d` works properly for a non-zero center.
802 #[test]
803 fn non_zero_center() {
804 for seeded in [false, true] {
805 let basis = MatZ::identity(5, 5);
806 let center = MatQ::identity(5, 1);
807 let gaussian_parameter = Q::ONE;
808 let basis_gso = MatQ::from(&basis).gso();
809
810 let _ = sample_d(&basis, ¢er, &gaussian_parameter, seeded.then(test_seed)).unwrap();
811 let _ = sample_d_precomputed_gso(
812 &basis,
813 &basis_gso,
814 ¢er,
815 &gaussian_parameter,
816 seeded.then(test_seed),
817 )
818 .unwrap();
819 }
820 }
821
822 /// Ensures that `sample_d` works properly for a different basis.
823 #[test]
824 fn non_identity_basis() {
825 for seeded in [false, true] {
826 let basis = MatZ::from_str("[[2, 1],[1, 2]]").unwrap();
827 let center = MatQ::new(2, 1);
828 let gaussian_parameter = Q::ONE;
829 let basis_gso = MatQ::from(&basis).gso();
830
831 let _ = sample_d(&basis, ¢er, &gaussian_parameter, seeded.then(test_seed)).unwrap();
832 let _ = sample_d_precomputed_gso(
833 &basis,
834 &basis_gso,
835 ¢er,
836 &gaussian_parameter,
837 seeded.then(test_seed),
838 )
839 .unwrap();
840 }
841 }
842
843 /// Ensures that `sample_d` outputs a vector that's part of the specified lattice.
844 ///
845 /// Checks whether the Hermite Normal Form HNF of the basis is equal to the HNF of
846 /// the basis concatenated with the sampled vector. If it is part of the lattice, it
847 /// should become a zero vector at the end of the matrix.
848 #[test]
849 fn point_of_lattice() {
850 for seeded in [false, true] {
851 let basis = MatZ::from_str("[[7, 0],[7, 3]]").unwrap();
852 let center = MatQ::new(2, 1);
853 let gaussian_parameter = Q::ONE;
854 let basis_gso = MatQ::from(&basis).gso();
855
856 let sample =
857 sample_d(&basis, ¢er, &gaussian_parameter, seeded.then(test_seed)).unwrap();
858 let sample_prec = sample_d_precomputed_gso(
859 &basis,
860 &basis_gso,
861 ¢er,
862 &gaussian_parameter,
863 seeded.then(test_seed),
864 )
865 .unwrap();
866
867 // check whether hermite normal form of HNF(b) = HNF([b|sample_vector])
868 let basis_concat_sample = basis.concat_horizontal(&sample).unwrap();
869 let basis_concat_sample_prec = basis.concat_horizontal(&sample_prec).unwrap();
870 let mut hnf_basis = MatZ::new(2, 2);
871 unsafe { fmpz_mat_hnf(&mut hnf_basis.matrix, &basis.matrix) };
872 let mut hnf_basis_concat_sample = MatZ::new(2, 3);
873 let mut hnf_basis_concat_sample_prec = MatZ::new(2, 3);
874 unsafe {
875 fmpz_mat_hnf(
876 &mut hnf_basis_concat_sample.matrix,
877 &basis_concat_sample.matrix,
878 )
879 };
880 unsafe {
881 fmpz_mat_hnf(
882 &mut hnf_basis_concat_sample_prec.matrix,
883 &basis_concat_sample_prec.matrix,
884 )
885 };
886 assert_eq!(
887 hnf_basis.get_column(0).unwrap(),
888 hnf_basis_concat_sample.get_column(0).unwrap()
889 );
890 assert_eq!(
891 hnf_basis.get_column(0).unwrap(),
892 hnf_basis_concat_sample_prec.get_column(0).unwrap()
893 );
894 assert_eq!(
895 hnf_basis.get_column(1).unwrap(),
896 hnf_basis_concat_sample.get_column(1).unwrap()
897 );
898 assert_eq!(
899 hnf_basis.get_column(1).unwrap(),
900 hnf_basis_concat_sample_prec.get_column(1).unwrap()
901 );
902 // check whether last vector is zero, i.e. was linearly dependent and part of lattice
903 assert!(hnf_basis_concat_sample.get_column(2).unwrap().is_zero());
904 assert!(
905 hnf_basis_concat_sample_prec
906 .get_column(2)
907 .unwrap()
908 .is_zero()
909 );
910 }
911 }
912
913 /// Checks whether `sample_d` returns an error if the gaussian parameter `s < 0`.
914 #[test]
915 fn invalid_gaussian_parameter() {
916 for seeded in [false, true] {
917 let basis = MatZ::identity(5, 5);
918 let center = MatQ::new(5, 1);
919 let basis_gso = MatQ::from(&basis).gso();
920
921 assert!(sample_d(&basis, ¢er, &Q::MINUS_ONE, seeded.then(test_seed)).is_err());
922 assert!(sample_d(&basis, ¢er, &Q::from(i64::MIN), seeded.then(test_seed)).is_err());
923
924 assert!(
925 sample_d_precomputed_gso(
926 &basis,
927 &basis_gso,
928 ¢er,
929 &Q::MINUS_ONE,
930 seeded.then(test_seed)
931 )
932 .is_err()
933 );
934 assert!(
935 sample_d_precomputed_gso(
936 &basis,
937 &basis_gso,
938 ¢er,
939 &Q::from(i64::MIN),
940 seeded.then(test_seed)
941 )
942 .is_err()
943 );
944 }
945 }
946
947 /// Checks whether `sample_d` returns an error if the basis and center number of rows differs.
948 #[test]
949 fn mismatching_matrix_dimensions() {
950 for seeded in [false, true] {
951 let basis = MatZ::identity(3, 5);
952 let center = MatQ::new(4, 1);
953 let gaussian_parameter = Q::ONE;
954 let basis_gso = MatQ::from(&basis).gso();
955
956 let res = sample_d(&basis, ¢er, &gaussian_parameter, seeded.then(test_seed));
957 let res_prec = sample_d_precomputed_gso(
958 &basis,
959 &basis_gso,
960 ¢er,
961 &gaussian_parameter,
962 seeded.then(test_seed),
963 );
964
965 assert!(res.is_err());
966 assert!(res_prec.is_err());
967 }
968 }
969
970 /// Checks whether `sample_d` returns an error if center isn't a column vector.
971 #[test]
972 fn center_not_column_vector() {
973 for seeded in [false, true] {
974 let basis = MatZ::identity(2, 2);
975 let center = MatQ::new(2, 2);
976 let gaussian_parameter = Q::ONE;
977 let basis_gso = MatQ::from(&basis).gso();
978
979 let res = sample_d(&basis, ¢er, &gaussian_parameter, seeded.then(test_seed));
980 let res_prec = sample_d_precomputed_gso(
981 &basis,
982 &basis_gso,
983 ¢er,
984 &gaussian_parameter,
985 seeded.then(test_seed),
986 );
987
988 assert!(res.is_err());
989 assert!(res_prec.is_err());
990 }
991 }
992
993 /// Ensures that the concentration bound holds.
994 #[test]
995 fn concentration_bound() {
996 for seeded in [false, true] {
997 let n = Z::from(20);
998 let basis = MatZ::sample_uniform(&n, &n, 0, 5000).unwrap();
999 let orth = MatQ::from(&basis).gso();
1000 let mut len = Q::ZERO;
1001 for i in 0..orth.get_num_columns() {
1002 let column = orth.get_column(i).unwrap();
1003 let column_len = column.norm_eucl_sqrd().unwrap().sqrt();
1004 if column_len > len {
1005 len = column_len
1006 }
1007 }
1008
1009 let expl_text = String::from("This test can fail with probability close to 0.
1010 It fails if the length of the sampled is longer than expected.
1011 If this happens, rerun the tests several times and check whether this issue comes up again.");
1012
1013 let center = MatQ::new(&n, 1);
1014 let gaussian_parameter =
1015 len * n.log(2).unwrap().sqrt() * (n.log(2).unwrap().log(2).unwrap());
1016
1017 for _ in 0..20 {
1018 let res =
1019 sample_d(&basis, ¢er, &gaussian_parameter, seeded.then(test_seed)).unwrap();
1020 let res_prec = sample_d_precomputed_gso(
1021 &basis,
1022 &orth,
1023 ¢er,
1024 &gaussian_parameter,
1025 seeded.then(test_seed),
1026 )
1027 .unwrap();
1028
1029 assert!(
1030 res.norm_eucl_sqrd().unwrap()
1031 <= gaussian_parameter.pow(2).unwrap().round() * &n,
1032 "{expl_text}"
1033 );
1034 assert!(
1035 res_prec.norm_eucl_sqrd().unwrap()
1036 <= gaussian_parameter.pow(2).unwrap().round() * &n,
1037 "{expl_text}"
1038 );
1039 }
1040 }
1041 }
1042
1043 /// Ensures that `sample_d` and `sample_d_precomputed_gso` output the same vector
1044 /// for the same seed.
1045 #[test]
1046 fn same_seed_same_sample() {
1047 let basis = MatZ::from_str("[[7, 0, 1],[7, 3, 2],[0, 1, 5]]").unwrap();
1048 let center = MatQ::from_str("[[1/2],[-3],[5]]").unwrap();
1049 let gaussian_parameter = Q::from(100);
1050 let basis_gso = MatQ::from(&basis).gso();
1051
1052 let sample_0 = sample_d(&basis, ¢er, &gaussian_parameter, Some([42; 32])).unwrap();
1053 let sample_1 = sample_d(&basis, ¢er, &gaussian_parameter, Some([42; 32])).unwrap();
1054 let sample_prec = sample_d_precomputed_gso(
1055 &basis,
1056 &basis_gso,
1057 ¢er,
1058 &gaussian_parameter,
1059 Some([42; 32]),
1060 )
1061 .unwrap();
1062
1063 assert_eq!(sample_0, sample_1);
1064 assert_eq!(sample_0, sample_prec);
1065 }
1066
1067 /// Ensure that an orthogonalized base with not matching rows panics.
1068 #[test]
1069 fn precomputed_gso_mismatching_rows() {
1070 for seeded in [false, true] {
1071 let result = std::panic::catch_unwind(|| {
1072 let n = Z::from(20);
1073 let basis = MatZ::sample_uniform(&n, &n, 0, 5000).unwrap();
1074 let center = MatQ::new(&n, 1);
1075 let false_gso = MatQ::new(basis.get_num_rows() + 1, basis.get_num_columns());
1076
1077 let _ = sample_d_precomputed_gso(
1078 &basis,
1079 &false_gso,
1080 ¢er,
1081 &Q::from(5),
1082 seeded.then(test_seed),
1083 )
1084 .unwrap();
1085 });
1086 assert!(result.is_err());
1087 }
1088 }
1089 /// Ensure that an orthogonalized base with not matching columns panics.
1090 #[test]
1091 fn precomputed_gso_mismatching_columns() {
1092 for seeded in [false, true] {
1093 let result = std::panic::catch_unwind(|| {
1094 let n = Z::from(20);
1095 let basis = MatZ::sample_uniform(&n, &n, 0, 5000).unwrap();
1096 let center = MatQ::new(&n, 1);
1097 let false_gso = MatQ::new(basis.get_num_rows(), basis.get_num_columns() + 1);
1098
1099 let _ = sample_d_precomputed_gso(
1100 &basis,
1101 &false_gso,
1102 ¢er,
1103 &Q::from(5),
1104 seeded.then(test_seed),
1105 )
1106 .unwrap();
1107 });
1108 assert!(result.is_err());
1109 }
1110 }
1111
1112 /// Checks whether `sample_d` outputs the known sample for a fixed seed.
1113 #[test]
1114 fn known_answer() {
1115 let basis = MatZ::from_str("[[7, 0, 1],[7, 3, 2],[0, 1, 5]]").unwrap();
1116 let center = MatQ::from_str("[[1/2],[-3],[5]]").unwrap();
1117
1118 let sample = sample_d(&basis, ¢er, &Q::from(100), Some([42; 32])).unwrap();
1119
1120 assert_eq!("[[28],[-44],[-24]]", sample.to_string());
1121 }
1122}