arnaucube
arnaucube/poulpyPublic Self-hosted repository
poulpy/math/src/ring/impl_u64/packing.rs
math/src/ring/impl_u64/packing.rsRaw
  1use crate::modulus::barrett::Barrett;
  2use crate::modulus::ONCE;
  3use crate::poly::Poly;
  4use crate::ring::Ring;
  5use std::cmp::min;
  6use std::collections::HashSet;
  7
  8impl Ring<u64> {
  9    // Generates a vector storing {X^{2^0}, X^{2^1}, .., X^{2^log_n}}.
 10    pub fn gen_x_pow_2<const NTT: bool, const INV: bool>(&self, log_n: usize) -> Vec<Poly<u64>> {
 11        let mut x_pow: Vec<Poly<u64>> = Vec::<Poly<u64>>::with_capacity(log_n);
 12
 13        (0..log_n).for_each(|i| {
 14            let mut idx: usize = 1 << i;
 15
 16            if INV {
 17                idx = self.n() - idx;
 18            }
 19
 20            x_pow.push(self.new_poly());
 21
 22            if i == 0 {
 23                x_pow[i].0[idx] = self.modulus.montgomery.one();
 24                self.ntt_inplace::<false>(&mut x_pow[i]);
 25            } else {
 26                let (left, right) = x_pow.split_at_mut(i);
 27                self.a_mul_b_montgomery_into_c::<ONCE>(&left[i - 1], &left[i - 1], &mut right[0]);
 28            }
 29        });
 30
 31        if INV {
 32            self.a_neg_into_a::<1, ONCE>(&mut x_pow[0]);
 33        }
 34
 35        if !NTT {
 36            x_pow.iter_mut().for_each(|x| self.intt_inplace::<false>(x));
 37        }
 38
 39        x_pow
 40    }
 41
 42    pub fn pack<const ZEROGARBAGE: bool, const NTT: bool>(
 43        &self,
 44        polys: &mut Vec<Option<Poly<u64>>>,
 45        log_gap: usize,
 46    ) {
 47        let log_n: usize = self.log_n();
 48        let log_nth_root: usize = log_n + 1;
 49        let nth_root: usize = 1 << log_nth_root;
 50        let log_start: usize = log_n - log_gap;
 51        let mut log_end: usize = log_n;
 52
 53        let mut indices: Vec<usize> = Vec::<usize>::new();
 54
 55        // Retrives non-empty indexes
 56        polys.iter().enumerate().for_each(|(i, poly)| {
 57            if !poly.is_none() {
 58                indices.push(i);
 59            }
 60        });
 61
 62        let gap: usize = max_gap(&indices);
 63
 64        let set: HashSet<_> = indices.into_iter().collect();
 65
 66        if !ZEROGARBAGE {
 67            if gap > 0 {
 68                log_end -= gap.trailing_zeros() as usize;
 69            }
 70        }
 71
 72        assert!(
 73            log_start < log_end,
 74            "invalid input polys: gap between non None value is smaller than 2^log_gap"
 75        );
 76
 77        let n_inv: Barrett<u64> = self
 78            .modulus
 79            .barrett
 80            .prepare(self.modulus.inv(1 << (log_end - log_start)));
 81
 82        set.iter().for_each(|i| {
 83            if let Some(poly) = polys[*i].as_mut() {
 84                if !NTT {
 85                    self.ntt_inplace::<true>(poly);
 86                }
 87                self.a_mul_b_scalar_barrett_into_a::<ONCE>(&n_inv, poly);
 88            }
 89        });
 90
 91        let x_pow2: Vec<Poly<u64>> = self.gen_x_pow_2::<true, false>(log_n);
 92        let mut tmpa: Poly<u64> = self.new_poly();
 93        let mut tmpb: Poly<u64> = self.new_poly();
 94
 95        for i in log_start..log_end {
 96            let t: usize = 1 << (log_n - 1 - i);
 97
 98            let (polys_lo, polys_hi) = polys.split_at_mut(t);
 99
100            for j in 0..t {
101                if let Some(poly_hi) = polys_hi[j].as_mut() {
102                    self.a_mul_b_montgomery_into_a::<ONCE>(&x_pow2[log_n - i - 1], poly_hi);
103
104                    if let Some(poly_lo) = polys_lo[j].as_mut() {
105                        self.a_sub_b_into_c::<1, ONCE>(poly_lo, poly_hi, &mut tmpa);
106                        self.a_add_b_into_b::<ONCE>(poly_hi, poly_lo);
107                    }
108                }
109
110                if let Some(poly_lo) = polys_lo[j].as_mut() {
111                    let gal_el: usize = self.galois_element((1 << i) >> 1, i == 0, log_nth_root);
112
113                    if !polys_hi[j].is_none() {
114                        self.automorphism::<true>(&tmpa, gal_el, 2 << self.log_n(), &mut tmpb);
115                        self.a_add_b_into_b::<ONCE>(&tmpb, poly_lo);
116                    } else {
117                        self.automorphism::<true>(poly_lo, gal_el, nth_root, &mut tmpa);
118                        self.a_add_b_into_b::<ONCE>(&tmpa, poly_lo);
119                    }
120                } else if let Some(poly_hi) = polys_hi[j].as_mut() {
121                    let gal_el: usize = self.galois_element((1 << i) >> 1, i == 0, log_nth_root);
122                    self.automorphism::<true>(poly_hi, gal_el, nth_root, &mut tmpa);
123                    self.a_sub_b_into_a::<1, ONCE>(&tmpa, poly_hi);
124                    std::mem::swap(&mut polys_lo[j], &mut polys_hi[j]);
125                }
126            }
127
128            polys.truncate(t);
129        }
130
131        if !NTT {
132            if let Some(poly) = polys[0].as_mut() {
133                self.intt_inplace::<false>(poly);
134            }
135        }
136    }
137}
138
139// Returns the largest gap.
140fn max_gap(vec: &[usize]) -> usize {
141    let mut gap: usize = usize::MAX;
142    for i in 1..vec.len() {
143        let (l, r) = (vec[i - 1], vec[i]);
144        assert!(r > l, "invalid input vec: not sorted");
145        gap = min(gap, r - l);
146        if gap == 1 {
147            break;
148        }
149    }
150    gap
151}