miden-crypto/src/dsa/rpo_falcon512/math/samplerz.rs
src/dsa/rpo_falcon512/math/samplerz.rsRaw
1use core::f64::consts::LN_2;
2use rand::Rng;
3
4#[cfg(not(feature = "std"))]
5use num::Float;
6
7/// Samples an integer from {0, ..., 18} according to the distribution χ, which is close to
8/// the half-Gaussian distribution on the natural numbers with mean 0 and standard deviation
9/// equal to sigma_max.
10fn base_sampler(bytes: [u8; 9]) -> i16 {
11 const RCDT: [u128; 18] = [
12 3024686241123004913666,
13 1564742784480091954050,
14 636254429462080897535,
15 199560484645026482916,
16 47667343854657281903,
17 8595902006365044063,
18 1163297957344668388,
19 117656387352093658,
20 8867391802663976,
21 496969357462633,
22 20680885154299,
23 638331848991,
24 14602316184,
25 247426747,
26 3104126,
27 28824,
28 198,
29 1,
30 ];
31 let u = u128::from_be_bytes([vec![0u8; 7], bytes.to_vec()].concat().try_into().unwrap());
32 RCDT.into_iter().filter(|r| u < *r).count() as i16
33}
34
35/// Computes an integer approximation of 2^63 * ccs * exp(-x).
36fn approx_exp(x: f64, ccs: f64) -> u64 {
37 // The constants C are used to approximate exp(-x); these
38 // constants are taken from FACCT (up to a scaling factor
39 // of 2^63):
40 // https://eprint.iacr.org/2018/1234
41 // https://github.com/raykzhao/gaussian
42 const C: [u64; 13] = [
43 0x00000004741183A3u64,
44 0x00000036548CFC06u64,
45 0x0000024FDCBF140Au64,
46 0x0000171D939DE045u64,
47 0x0000D00CF58F6F84u64,
48 0x000680681CF796E3u64,
49 0x002D82D8305B0FEAu64,
50 0x011111110E066FD0u64,
51 0x0555555555070F00u64,
52 0x155555555581FF00u64,
53 0x400000000002B400u64,
54 0x7FFFFFFFFFFF4800u64,
55 0x8000000000000000u64,
56 ];
57
58 let mut z: u64;
59 let mut y: u64;
60 let twoe63 = 1u64 << 63;
61
62 y = C[0];
63 z = f64::floor(x * (twoe63 as f64)) as u64;
64 for cu in C.iter().skip(1) {
65 let zy = (z as u128) * (y as u128);
66 y = cu - ((zy >> 63) as u64);
67 }
68
69 z = f64::floor((twoe63 as f64) * ccs) as u64;
70
71 (((z as u128) * (y as u128)) >> 63) as u64
72}
73
74/// A random bool that is true with probability ≈ ccs · exp(-x).
75fn ber_exp(x: f64, ccs: f64, random_bytes: [u8; 7]) -> bool {
76 // 0.69314718055994530941 = ln(2)
77 let s = f64::floor(x / LN_2) as usize;
78 let r = x - LN_2 * (s as f64);
79 let shamt = usize::min(s, 63);
80 let z = ((((approx_exp(r, ccs) as u128) << 1) - 1) >> shamt) as u64;
81 let mut w = 0i16;
82 for (index, i) in (0..64).step_by(8).rev().enumerate() {
83 let byte = random_bytes[index];
84 w = (byte as i16) - (((z >> i) & 0xff) as i16);
85 if w != 0 {
86 break;
87 }
88 }
89 w < 0
90}
91
92/// Samples an integer from the Gaussian distribution with given mean (mu) and standard deviation
93/// (sigma).
94pub(crate) fn sampler_z<R: Rng>(mu: f64, sigma: f64, sigma_min: f64, rng: &mut R) -> i16 {
95 const SIGMA_MAX: f64 = 1.8205;
96 const INV_2SIGMA_MAX_SQ: f64 = 1f64 / (2f64 * SIGMA_MAX * SIGMA_MAX);
97 let isigma = 1f64 / sigma;
98 let dss = 0.5f64 * isigma * isigma;
99 let s = f64::floor(mu);
100 let r = mu - s;
101 let ccs = sigma_min * isigma;
102 loop {
103 let z0 = base_sampler(rng.gen());
104 let random_byte: u8 = rng.gen();
105 let b = (random_byte & 1) as i16;
106 let z = b + ((b << 1) - 1) * z0;
107 let zf_min_r = (z as f64) - r;
108 // x = ((z-r)^2)/(2*sigma^2) - ((z-b)^2)/(2*sigma0^2)
109 let x = zf_min_r * zf_min_r * dss - (z0 * z0) as f64 * INV_2SIGMA_MAX_SQ;
110 if ber_exp(x, ccs, rng.gen()) {
111 return z + (s as i16);
112 }
113 }
114}
115
116#[cfg(all(test, feature = "std"))]
117mod test {
118 use alloc::vec::Vec;
119 use rand::RngCore;
120 use std::{thread::sleep, time::Duration};
121
122 use super::{approx_exp, ber_exp, sampler_z};
123
124 /// RNG used only for testing purposes, whereby the produced
125 /// string of random bytes is equal to the one it is initialized
126 /// with. Whatever you do, do not use this RNG in production.
127 struct UnsafeBufferRng {
128 buffer: Vec<u8>,
129 index: usize,
130 }
131
132 impl UnsafeBufferRng {
133 fn new(buffer: &[u8]) -> Self {
134 Self { buffer: buffer.to_vec(), index: 0 }
135 }
136
137 fn next(&mut self) -> u8 {
138 if self.buffer.len() <= self.index {
139 // panic!("Ran out of buffer.");
140 sleep(Duration::from_millis(10));
141 0u8
142 } else {
143 let return_value = self.buffer[self.index];
144 self.index += 1;
145 return_value
146 }
147 }
148 }
149
150 impl RngCore for UnsafeBufferRng {
151 fn next_u32(&mut self) -> u32 {
152 // let bytes: [u8; 4] = (0..4)
153 // .map(|_| self.next())
154 // .collect_vec()
155 // .try_into()
156 // .unwrap();
157 // u32::from_be_bytes(bytes)
158 u32::from_le_bytes([self.next(), 0, 0, 0])
159 }
160
161 fn next_u64(&mut self) -> u64 {
162 // let bytes: [u8; 8] = (0..8)
163 // .map(|_| self.next())
164 // .collect_vec()
165 // .try_into()
166 // .unwrap();
167 // u64::from_be_bytes(bytes)
168 u64::from_le_bytes([self.next(), 0, 0, 0, 0, 0, 0, 0])
169 }
170
171 fn fill_bytes(&mut self, dest: &mut [u8]) {
172 for d in dest.iter_mut() {
173 *d = self.next();
174 }
175 }
176
177 fn try_fill_bytes(&mut self, dest: &mut [u8]) -> Result<(), rand::Error> {
178 for d in dest.iter_mut() {
179 *d = self.next();
180 }
181 Ok(())
182 }
183 }
184
185 #[test]
186 fn test_unsafe_buffer_rng() {
187 let seed_bytes = hex::decode("7FFECD162AE2").unwrap();
188 let mut rng = UnsafeBufferRng::new(&seed_bytes);
189 let generated_bytes: Vec<u8> = (0..seed_bytes.len()).map(|_| rng.next()).collect();
190 assert_eq!(seed_bytes, generated_bytes);
191 }
192
193 #[test]
194 fn test_approx_exp() {
195 let precision = 1u64 << 14;
196 // known answers were generated with the following sage script:
197 //```sage
198 // num_samples = 10
199 // precision = 200
200 // R = Reals(precision)
201 //
202 // print(f"let kats : [(f64, f64, u64);{num_samples}] = [")
203 // for i in range(num_samples):
204 // x = RDF.random_element(0.0, 0.693147180559945)
205 // ccs = RDF.random_element(0.0, 1.0)
206 // res = round(2^63 * R(ccs) * exp(R(-x)))
207 // print(f"({x}, {ccs}, {res}),")
208 // print("];")
209 // ```
210 let kats: [(f64, f64, u64); 10] = [
211 (0.2314993926072656, 0.8148006314615972, 5962140072160879737),
212 (0.2648875572812225, 0.12769669655309035, 903712282351034505),
213 (0.11251957513682391, 0.9264611470305881, 7635725498677341553),
214 (0.04353439307256617, 0.5306497137523327, 4685877322232397936),
215 (0.41834495299784347, 0.879438856118578, 5338392138535350986),
216 (0.32579398973228557, 0.16513412873289002, 1099603299296456803),
217 (0.5939508073919817, 0.029776019144967303, 151637565622779016),
218 (0.2932367999399056, 0.37123847662857923, 2553827649386670452),
219 (0.5005699297417507, 0.31447208863888976, 1758235618083658825),
220 (0.4876437338498085, 0.6159515298936868, 3488632981903743976),
221 ];
222 for (x, ccs, answer) in kats {
223 let difference = (answer as i128) - (approx_exp(x, ccs) as i128);
224 assert!(
225 (difference * difference) as u64 <= precision * precision,
226 "answer: {answer} versus approximation: {}\ndifference: {} whereas precision: {}",
227 approx_exp(x, ccs),
228 difference,
229 precision
230 );
231 }
232 }
233
234 #[test]
235 fn test_ber_exp() {
236 let kats = [
237 (
238 1.268_314_048_020_498_4,
239 0.749_990_853_267_664_9,
240 hex::decode("ea000000000000").unwrap(),
241 false,
242 ),
243 (
244 0.001_563_917_959_143_409_6,
245 0.749_990_853_267_664_9,
246 hex::decode("6c000000000000").unwrap(),
247 true,
248 ),
249 (
250 0.017_921_215_753_999_235,
251 0.749_990_853_267_664_9,
252 hex::decode("c2000000000000").unwrap(),
253 false,
254 ),
255 (
256 0.776_117_648_844_980_6,
257 0.751_181_554_542_520_8,
258 hex::decode("58000000000000").unwrap(),
259 true,
260 ),
261 ];
262 for (x, ccs, bytes, answer) in kats {
263 assert_eq!(answer, ber_exp(x, ccs, bytes.try_into().unwrap()));
264 }
265 }
266
267 #[test]
268 fn test_sampler_z() {
269 let sigma_min = 1.277833697;
270 // known answers from the doc, table 3.2, page 44
271 // https://falcon-sign.info/falcon.pdf
272 // The zeros were added to account for dropped bytes.
273 let kats = [
274 (-91.90471153063714,1.7037990414754918,hex::decode("0fc5442ff043d66e91d1ea000000000000cac64ea5450a22941edc6c").unwrap(),-92),
275 (-8.322564895434937,1.7037990414754918,hex::decode("f4da0f8d8444d1a77265c2000000000000ef6f98bbbb4bee7db8d9b3").unwrap(),-8),
276 (-19.096516109216804,1.7035823083824078,hex::decode("db47f6d7fb9b19f25c36d6000000000000b9334d477a8bc0be68145d").unwrap(),-20),
277 (-11.335543982423326, 1.7035823083824078, hex::decode("ae41b4f5209665c74d00dc000000000000c1a8168a7bb516b3190cb42c1ded26cd52000000000000aed770eca7dd334e0547bcc3c163ce0b").unwrap(), -12),
278 (7.9386734193997555, 1.6984647769450156, hex::decode("31054166c1012780c603ae0000000000009b833cec73f2f41ca5807c000000000000c89c92158834632f9b1555").unwrap(), 8),
279 (-28.990850086867255, 1.6984647769450156, hex::decode("737e9d68a50a06dbbc6477").unwrap(), -30),
280 (-9.071257914091655, 1.6980782114808988, hex::decode("a98ddd14bf0bf22061d632").unwrap(), -10),
281 (-43.88754568839566, 1.6980782114808988, hex::decode("3cbf6818a68f7ab9991514").unwrap(), -41),
282 (-58.17435547946095,1.7010983419195522,hex::decode("6f8633f5bfa5d26848668e0000000000003d5ddd46958e97630410587c").unwrap(),-61),
283 (-43.58664906684732, 1.7010983419195522, hex::decode("272bc6c25f5c5ee53f83c40000000000003a361fbc7cc91dc783e20a").unwrap(), -46),
284 (-34.70565203313315, 1.7009387219711465, hex::decode("45443c59574c2c3b07e2e1000000000000d9071e6d133dbe32754b0a").unwrap(), -34),
285 (-44.36009577368896, 1.7009387219711465, hex::decode("6ac116ed60c258e2cbaeab000000000000728c4823e6da36e18d08da0000000000005d0cc104e21cc7fd1f5ca8000000000000d9dbb675266c928448059e").unwrap(), -44),
286 (-21.783037079346236, 1.6958406126012802, hex::decode("68163bc1e2cbf3e18e7426").unwrap(), -23),
287 (-39.68827784633828, 1.6958406126012802, hex::decode("d6a1b51d76222a705a0259").unwrap(), -40),
288 (-18.488607061056847, 1.6955259305261838, hex::decode("f0523bfaa8a394bf4ea5c10000000000000f842366fde286d6a30803").unwrap(), -22),
289 (-48.39610939101591, 1.6955259305261838, hex::decode("87bd87e63374cee62127fc0000000000006931104aab64f136a0485b").unwrap(), -50),
290 ];
291 for (mu, sigma, random_bytes, answer) in kats {
292 assert_eq!(
293 sampler_z(mu, sigma, sigma_min, &mut UnsafeBufferRng::new(&random_bytes)),
294 answer
295 );
296 }
297 }
298}