poulpy/base2k/src/stats.rs
base2k/src/stats.rsRaw
1use crate::{Encoding, VecZnx, ZnxInfos};
2use rug::Float;
3use rug::float::Round;
4use rug::ops::{AddAssignRound, DivAssignRound, SubAssignRound};
5
6pub trait Stats {
7 /// Returns the standard devaition of the i-th polynomial.
8 fn std(&self, col_i: usize, log_base2k: usize) -> f64;
9}
10
11impl Stats for VecZnx {
12 fn std(&self, col_i: usize, log_base2k: usize) -> f64 {
13 let prec: u32 = (self.limbs() * log_base2k) as u32;
14 let mut data: Vec<Float> = (0..self.n()).map(|_| Float::with_val(prec, 0)).collect();
15 self.decode_vec_float(col_i, log_base2k, &mut data);
16 // std = sqrt(sum((xi - avg)^2) / n)
17 let mut avg: Float = Float::with_val(prec, 0);
18 data.iter().for_each(|x| {
19 avg.add_assign_round(x, Round::Nearest);
20 });
21 avg.div_assign_round(Float::with_val(prec, data.len()), Round::Nearest);
22 data.iter_mut().for_each(|x| {
23 x.sub_assign_round(&avg, Round::Nearest);
24 });
25 let mut std: Float = Float::with_val(prec, 0);
26 data.iter().for_each(|x| std += x * x);
27 std.div_assign_round(Float::with_val(prec, data.len()), Round::Nearest);
28 std = std.sqrt();
29 std.to_f64()
30 }
31}