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