discrete_gaussian/
vtime.rs

1//! Variable time discrete sampling over algebraic structures.
2//!
3use rand::distributions::Bernoulli;
4use rand::distributions::Distribution;
5
6use super::ctime::pow2_1024;
7use crate::THETA_0;
8
9/// Sample a u32 integer from a gaussian distribution. Automatically
10/// approximates a `k` parameter based on the requested theta.
11pub fn sample_vartime<R: rand::Rng>(theta: f64, rng: &mut R) -> u32 {
12    // Approximate a k value
13    let k: u32 = super::k_from_theta(theta);
14    sample_vartime_k(k, rng)
15}
16
17/// Vartime u32 gaussian sampling based on [DDLL13](https://eprint.iacr.org/2013/383.pdf).
18pub fn sample_vartime_k<R: rand::Rng>(k: u32, rng: &mut R) -> u32 {
19    let bits = 32;
20    // check that k is <= 2^52? somehow without doing it vartime?
21    let theta: f64 = f64::from(k) * THETA_0;
22    let x = sample_theta_0_vartime(rng);
23    let y = rng.gen_range(0..k);
24    let z = u64::from(k) * u64::from(x) + u64::from(y);
25    let t = u64::from(y) * (u64::from(y) + 2 * u64::from(k) * u64::from(x));
26    let mut b = 1;
27    for i in (0..bits).rev() {
28        let p_i = euler_50_approx(-1.0 * pow2_1024(f64::from(i) / (2.0 * theta * theta)));
29        let t_i = (t >> i) & 1;
30        let d = Bernoulli::new(p_i).unwrap();
31        let v = d.sample(&mut rand::thread_rng());
32        b = b * (1 - t_i + u64::from(v) * t_i);
33    }
34    if b == 0 {
35        return sample_vartime_k(k, rng);
36    }
37    u32::try_from(z).unwrap()
38}
39
40/// Vartime probabilistic sampling from the theta_0 discrete
41/// gaussian distribution evaluated over the positive integers.
42///
43/// `theta_0 = sqrt(1/(2*ln(2)))` as defined in [FACCT](https://eprint.iacr.org/2018/1234.pdf)
44/// page 6
45pub fn sample_theta_0_vartime<R: rand::Rng>(rng: &mut R) -> u32 {
46    let b = rng.gen_range(0..=1);
47    if b == 0 {
48        return 0;
49    }
50    let mut i = 1;
51    loop {
52        for _ in 0..(2 * i - 2) {
53            if rng.gen_range(0..=1) != 0 {
54                return sample_theta_0_vartime(rng);
55            }
56        }
57        if rng.gen_range(0..=1) == 0 {
58            return i;
59        }
60        i += 1;
61    }
62}
63
64/// 50 bit euler approximation in f64
65/// returns `e^x`
66///
67/// variable time implementation due to division operator
68pub fn euler_50_approx(x: f64) -> f64 {
69    let p1: f64 = 1.66666666666666019037 * 10_f64.powf(-1_f64);
70    let p2: f64 = -2.77777777770155933842 * 10_f64.powf(-3_f64);
71    let p3: f64 = 6.61375632143793436117 * 10_f64.powf(-5_f64);
72    let p4: f64 = -1.65339022054652515390_f64 * 10_f64.powf(-6_f64);
73    let p5 = 4.13813679705723846039 * 10_f64.powf(-8_f64);
74    let s = x / 2_f64;
75    let t = s * s;
76    // Let c = s−t·(p1 +t·(p2 +t·(p3 +t·(p4 +t·p5))))
77    // explicitly write horner evaluation i guess
78    let c = s - t * (p1 + t * (p2 + t * (p3 + t * (p4 + t * p5))));
79    // Let r = 1 − ((s · c) / (c − 2) − s).
80    let r = 1_f64 - ((s * c) / (c - 2_f64) - s);
81    return r * r;
82}
83
84#[test]
85fn generate_samples_with_k() {
86    for _ in 0..1000 {
87        sample_vartime_k(40, &mut rand::thread_rng());
88    }
89}
90
91#[test]
92fn generate_samples() {
93    for _ in 0..1000 {
94        sample_vartime(10.0, &mut rand::thread_rng());
95    }
96}
discrete_gaussian/vtime.rs

discrete_gaussian/
vtime.rs
