miden-crypto/src/dsa/rpo_falcon512/math/samplerz.rs
src/dsa/rpo_falcon512/math/samplerz.rsRaw
1use core::f64::consts::LN_2;
2
3#[cfg(not(feature = "std"))]
4use num::Float;
5use rand::Rng;
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 std::{thread::sleep, time::Duration};
120
121 use rand::RngCore;
122
123 use super::{approx_exp, ber_exp, sampler_z};
124
125 /// RNG used only for testing purposes, whereby the produced
126 /// string of random bytes is equal to the one it is initialized
127 /// with. Whatever you do, do not use this RNG in production.
128 struct UnsafeBufferRng {
129 buffer: Vec<u8>,
130 index: usize,
131 }
132
133 impl UnsafeBufferRng {
134 fn new(buffer: &[u8]) -> Self {
135 Self { buffer: buffer.to_vec(), index: 0 }
136 }
137
138 fn next(&mut self) -> u8 {
139 if self.buffer.len() <= self.index {
140 // panic!("Ran out of buffer.");
141 sleep(Duration::from_millis(10));
142 0u8
143 } else {
144 let return_value = self.buffer[self.index];
145 self.index += 1;
146 return_value
147 }
148 }
149 }
150
151 impl RngCore for UnsafeBufferRng {
152 fn next_u32(&mut self) -> u32 {
153 // let bytes: [u8; 4] = (0..4)
154 // .map(|_| self.next())
155 // .collect_vec()
156 // .try_into()
157 // .unwrap();
158 // u32::from_be_bytes(bytes)
159 u32::from_le_bytes([self.next(), 0, 0, 0])
160 }
161
162 fn next_u64(&mut self) -> u64 {
163 // let bytes: [u8; 8] = (0..8)
164 // .map(|_| self.next())
165 // .collect_vec()
166 // .try_into()
167 // .unwrap();
168 // u64::from_be_bytes(bytes)
169 u64::from_le_bytes([self.next(), 0, 0, 0, 0, 0, 0, 0])
170 }
171
172 fn fill_bytes(&mut self, dest: &mut [u8]) {
173 for d in dest.iter_mut() {
174 *d = self.next();
175 }
176 }
177
178 fn try_fill_bytes(&mut self, dest: &mut [u8]) -> Result<(), rand::Error> {
179 for d in dest.iter_mut() {
180 *d = self.next();
181 }
182 Ok(())
183 }
184 }
185
186 #[test]
187 fn test_unsafe_buffer_rng() {
188 let seed_bytes = hex::decode("7FFECD162AE2").unwrap();
189 let mut rng = UnsafeBufferRng::new(&seed_bytes);
190 let generated_bytes: Vec<u8> = (0..seed_bytes.len()).map(|_| rng.next()).collect();
191 assert_eq!(seed_bytes, generated_bytes);
192 }
193
194 #[test]
195 fn test_approx_exp() {
196 let precision = 1u64 << 14;
197 // known answers were generated with the following sage script:
198 //```sage
199 // num_samples = 10
200 // precision = 200
201 // R = Reals(precision)
202 //
203 // print(f"let kats : [(f64, f64, u64);{num_samples}] = [")
204 // for i in range(num_samples):
205 // x = RDF.random_element(0.0, 0.693147180559945)
206 // ccs = RDF.random_element(0.0, 1.0)
207 // res = round(2^63 * R(ccs) * exp(R(-x)))
208 // print(f"({x}, {ccs}, {res}),")
209 // print("];")
210 // ```
211 let kats: [(f64, f64, u64); 10] = [
212 (0.2314993926072656, 0.8148006314615972, 5962140072160879737),
213 (0.2648875572812225, 0.12769669655309035, 903712282351034505),
214 (0.11251957513682391, 0.9264611470305881, 7635725498677341553),
215 (0.04353439307256617, 0.5306497137523327, 4685877322232397936),
216 (0.41834495299784347, 0.879438856118578, 5338392138535350986),
217 (0.32579398973228557, 0.16513412873289002, 1099603299296456803),
218 (0.5939508073919817, 0.029776019144967303, 151637565622779016),
219 (0.2932367999399056, 0.37123847662857923, 2553827649386670452),
220 (0.5005699297417507, 0.31447208863888976, 1758235618083658825),
221 (0.4876437338498085, 0.6159515298936868, 3488632981903743976),
222 ];
223 for (x, ccs, answer) in kats {
224 let difference = (answer as i128) - (approx_exp(x, ccs) as i128);
225 assert!(
226 (difference * difference) as u64 <= precision * precision,
227 "answer: {answer} versus approximation: {}\ndifference: {} whereas precision: {}",
228 approx_exp(x, ccs),
229 difference,
230 precision
231 );
232 }
233 }
234
235 #[test]
236 fn test_ber_exp() {
237 let kats = [
238 (
239 1.268_314_048_020_498_4,
240 0.749_990_853_267_664_9,
241 hex::decode("ea000000000000").unwrap(),
242 false,
243 ),
244 (
245 0.001_563_917_959_143_409_6,
246 0.749_990_853_267_664_9,
247 hex::decode("6c000000000000").unwrap(),
248 true,
249 ),
250 (
251 0.017_921_215_753_999_235,
252 0.749_990_853_267_664_9,
253 hex::decode("c2000000000000").unwrap(),
254 false,
255 ),
256 (
257 0.776_117_648_844_980_6,
258 0.751_181_554_542_520_8,
259 hex::decode("58000000000000").unwrap(),
260 true,
261 ),
262 ];
263 for (x, ccs, bytes, answer) in kats {
264 assert_eq!(answer, ber_exp(x, ccs, bytes.try_into().unwrap()));
265 }
266 }
267
268 #[test]
269 fn test_sampler_z() {
270 let sigma_min = 1.277833697;
271 // known answers from the doc, table 3.2, page 44
272 // https://falcon-sign.info/falcon.pdf
273 // The zeros were added to account for dropped bytes.
274 let kats = [
275 (-91.90471153063714,1.7037990414754918,hex::decode("0fc5442ff043d66e91d1ea000000000000cac64ea5450a22941edc6c").unwrap(),-92),
276 (-8.322564895434937,1.7037990414754918,hex::decode("f4da0f8d8444d1a77265c2000000000000ef6f98bbbb4bee7db8d9b3").unwrap(),-8),
277 (-19.096516109216804,1.7035823083824078,hex::decode("db47f6d7fb9b19f25c36d6000000000000b9334d477a8bc0be68145d").unwrap(),-20),
278 (-11.335543982423326, 1.7035823083824078, hex::decode("ae41b4f5209665c74d00dc000000000000c1a8168a7bb516b3190cb42c1ded26cd52000000000000aed770eca7dd334e0547bcc3c163ce0b").unwrap(), -12),
279 (7.9386734193997555, 1.6984647769450156, hex::decode("31054166c1012780c603ae0000000000009b833cec73f2f41ca5807c000000000000c89c92158834632f9b1555").unwrap(), 8),
280 (-28.990850086867255, 1.6984647769450156, hex::decode("737e9d68a50a06dbbc6477").unwrap(), -30),
281 (-9.071257914091655, 1.6980782114808988, hex::decode("a98ddd14bf0bf22061d632").unwrap(), -10),
282 (-43.88754568839566, 1.6980782114808988, hex::decode("3cbf6818a68f7ab9991514").unwrap(), -41),
283 (-58.17435547946095,1.7010983419195522,hex::decode("6f8633f5bfa5d26848668e0000000000003d5ddd46958e97630410587c").unwrap(),-61),
284 (-43.58664906684732, 1.7010983419195522, hex::decode("272bc6c25f5c5ee53f83c40000000000003a361fbc7cc91dc783e20a").unwrap(), -46),
285 (-34.70565203313315, 1.7009387219711465, hex::decode("45443c59574c2c3b07e2e1000000000000d9071e6d133dbe32754b0a").unwrap(), -34),
286 (-44.36009577368896, 1.7009387219711465, hex::decode("6ac116ed60c258e2cbaeab000000000000728c4823e6da36e18d08da0000000000005d0cc104e21cc7fd1f5ca8000000000000d9dbb675266c928448059e").unwrap(), -44),
287 (-21.783037079346236, 1.6958406126012802, hex::decode("68163bc1e2cbf3e18e7426").unwrap(), -23),
288 (-39.68827784633828, 1.6958406126012802, hex::decode("d6a1b51d76222a705a0259").unwrap(), -40),
289 (-18.488607061056847, 1.6955259305261838, hex::decode("f0523bfaa8a394bf4ea5c10000000000000f842366fde286d6a30803").unwrap(), -22),
290 (-48.39610939101591, 1.6955259305261838, hex::decode("87bd87e63374cee62127fc0000000000006931104aab64f136a0485b").unwrap(), -50),
291 ];
292 for (mu, sigma, random_bytes, answer) in kats {
293 assert_eq!(
294 sampler_z(mu, sigma, sigma_min, &mut UnsafeBufferRng::new(&random_bytes)),
295 answer
296 );
297 }
298 }
299}