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::montgomery::Montgomery;
  3use crate::modulus::ONCE;
  4use crate::poly::Poly;
  5use crate::ring::Ring;
  6use std::cmp::min;
  7use std::rc::Rc;
  8
  9impl Ring<u64> {
 10
 11    pub fn pack<const ZEROGARBAGE: bool, const NTT: bool>(
 12        &self,
 13        polys: &mut Vec<Option<Poly<u64>>>,
 14        log_gap: usize,
 15    ) {
 16        let log_n: usize = self.log_n();
 17        let log_nth_root: usize = log_n + 1;
 18        let nth_root: usize = 1 << log_nth_root;
 19        let log_start: usize = log_n - log_gap;
 20        let mut log_end: usize = log_n;
 21
 22        let mut indices: Vec<usize> = Vec::<usize>::new();
 23
 24        // Retrives non-empty indexes
 25        polys.iter().enumerate().for_each(|(i, poly)| {
 26            if !poly.is_none() {
 27                indices.push(i);
 28            }
 29        });
 30
 31        let gap: usize = max_gap(&indices);
 32
 33        if !ZEROGARBAGE {
 34            if gap > 0 {
 35                log_end -= gap.trailing_zeros() as usize;
 36            }
 37        }
 38
 39        assert!(
 40            log_start < log_end,
 41            "invalid input polys: gap between non None value is smaller than 2^log_gap"
 42        );
 43
 44        let n_inv: Barrett<u64> = self
 45            .modulus
 46            .barrett
 47            .prepare(self.modulus.inv(1 << (log_end - log_start)));
 48
 49        indices.iter().for_each(|i| {
 50            if let Some(poly) = polys[*i].as_mut() {
 51                if !NTT {
 52                    self.ntt_inplace::<true>(poly);
 53                }
 54                self.a_mul_b_scalar_barrett_into_a::<ONCE>(&n_inv, poly);
 55            }
 56        });
 57
 58        let x_pow2: Vec<Poly<Montgomery<u64>>> = self.gen_x_pow_2::<true, false>(log_n);
 59        let mut tmpa: Poly<u64> = self.new_poly();
 60        let mut tmpb: Poly<u64> = self.new_poly();
 61
 62        for i in log_start..log_end {
 63            let t: usize = 1 << (log_n - 1 - i);
 64
 65            let (polys_lo, polys_hi) = polys.split_at_mut(t);
 66
 67            for j in 0..t {
 68                if let Some(poly_hi) = polys_hi[j].as_mut() {
 69                    self.a_mul_b_montgomery_into_a::<ONCE>(&x_pow2[log_n - i - 1], poly_hi);
 70
 71                    if let Some(poly_lo) = polys_lo[j].as_mut() {
 72                        self.a_sub_b_into_c::<1, ONCE>(poly_lo, poly_hi, &mut tmpa);
 73                        self.a_add_b_into_b::<ONCE>(poly_hi, poly_lo);
 74                    }
 75                }
 76
 77                if let Some(poly_lo) = polys_lo[j].as_mut() {
 78                    let gal_el: usize = self.galois_element((1 << i) >> 1, i == 0, log_nth_root);
 79
 80                    if !polys_hi[j].is_none() {
 81                        self.a_apply_automorphism_add_b_into_b::<ONCE, true>(&tmpa, gal_el, 2 << self.log_n(), poly_lo);
 82                    } else {
 83                        self.a_apply_automorphism_into_b::<true>(poly_lo, gal_el, nth_root, &mut tmpa);
 84                        self.a_add_b_into_b::<ONCE>(&tmpa, poly_lo);
 85                    }
 86                } else if let Some(poly_hi) = polys_hi[j].as_mut() {
 87                    let gal_el: usize = self.galois_element((1 << i) >> 1, i == 0, log_nth_root);
 88                    self.a_apply_automorphism_into_b::<true>(poly_hi, gal_el, nth_root, &mut tmpa);
 89                    self.a_sub_b_into_a::<1, ONCE>(&tmpa, poly_hi);
 90                    std::mem::swap(&mut polys_lo[j], &mut polys_hi[j]);
 91                }
 92            }
 93
 94            polys.truncate(t);
 95        }
 96
 97        if !NTT {
 98            if let Some(poly) = polys[0].as_mut() {
 99                self.intt_inplace::<false>(poly);
100            }
101        }
102    }
103}
104
105
106// Returns the largest gap between two values in an ordered array of distinct values.
107// Panics if the array is not ordered or values are not distincts.
108fn max_gap(vec: &[usize]) -> usize {
109    let mut gap: usize = usize::MAX;
110    for i in 1..vec.len() {
111        let (l, r) = (vec[i - 1], vec[i]);
112        assert!(r > l, "invalid input vec: not sorted or collision between indices");
113        gap = min(gap, r - l);
114        if gap == 1 {
115            break;
116        }
117    }
118    gap
119}
120
121
122pub struct StreamRepacker{
123    accumulators: Vec<Accumulator>,
124    buf0: Poly<u64>,
125    buf1: Poly<u64>,
126    buf_auto: Poly<u64>,
127    x_pow_2: Vec<Poly<Montgomery<u64>>>,
128    n_inv: Barrett<u64>,
129    counter: usize,
130}
131
132pub struct Accumulator{
133    buf: [Option<Rc<Poly<u64>>>; 2],
134    control: bool,
135}
136
137impl Accumulator{
138    pub fn new(r: &Ring<u64>) -> Self{
139        Self { buf: [Some(Rc::new(r.new_poly())), None], control: false }
140    }
141}
142
143impl StreamRepacker{
144    pub fn new(r: &Ring<u64>) -> Self{
145
146        let mut accumulators: Vec<Accumulator> = Vec::<Accumulator>::new();
147
148        (0..r.log_n()).for_each(|_|
149            accumulators.push(Accumulator::new(r))
150        );
151
152        Self{
153            accumulators: accumulators,
154            buf0: r.new_poly(),
155            buf1: r.new_poly(),
156            buf_auto: r.new_poly(),
157            x_pow_2: r.gen_x_pow_2::<true, false>(r.log_n()),
158            n_inv: r.modulus.barrett.prepare(r.modulus.inv(r.n() as u64)),
159            counter:0,
160        }
161    }
162
163    fn merge_ab(&mut self, r: &Ring<u64>, a: &Poly<u64>, b: &Poly<u64>, i: usize) -> &Poly<u64>{
164
165        let tmp_a: &mut Poly<u64> = &mut self.buf0;
166        let tmp_b: &mut Poly<u64> = &mut self.buf1;
167
168        r.a_mul_b_montgomery_into_c::<ONCE>(a, &self.x_pow_2[r.log_n()-i-1], tmp_a);
169        r.a_sub_b_into_c::<1, ONCE>(a, tmp_a, tmp_b);
170        r.a_add_b_into_b::<ONCE>(a, tmp_a);
171
172        if i == 0{
173            r.a_mul_b_scalar_barrett_into_a::<ONCE>(&self.n_inv, tmp_a);
174            r.a_mul_b_scalar_barrett_into_a::<ONCE>(&self.n_inv, tmp_b);
175        }
176
177        let log_nth_root = r.log_n()+1;
178        let nth_root = 1<<log_nth_root;
179            
180        let gal_el: usize = r.galois_element((1 << i) >> 1, i == 0, log_nth_root);
181
182        r.a_apply_automorphism_add_b_into_b::<ONCE, true>(tmp_b, gal_el, nth_root, tmp_a);
183
184        tmp_a
185    }
186
187    fn merge_a(&mut self, r: &Ring<u64>, a: &Poly<u64>, i: usize) -> &Poly<u64>{
188
189        let tmp_a: &mut Poly<u64> = &mut self.buf0;
190 
191        let log_nth_root = r.log_n()+1;
192        let nth_root = 1<<log_nth_root;
193        let gal_el: usize = r.galois_element((1 << i) >> 1, i == 0, log_nth_root);
194
195        if i == 0{
196            r.a_mul_b_scalar_barrett_into_a::<ONCE>(&self.n_inv, tmp_a);
197            r.a_apply_automorphism_into_b::<true>(tmp_a, gal_el, nth_root, &mut self.buf_auto)
198            r.a_add_b_into_c::<ONCE>(&self.buf_auto, a, tmp_a);
199        }else{
200            r.a_apply_automorphism_into_b::<true>(a, gal_el, nth_root, tmp_a);
201            r.a_add_b_into_b::<ONCE>(a, tmp_a);
202        }
203
204        tmp_a
205    }
206
207    fn merge_b(&mut self, r: &Ring<u64>, b: &Poly<u64>, i: usize) -> &Poly<u64>{
208
209        let tmp_a: &mut Poly<u64> = &mut self.buf0;
210        let tmp_b: &mut Poly<u64> = &mut self.buf1;
211
212        let log_nth_root = r.log_n()+1;
213        let nth_root = 1<<log_nth_root;
214        let gal_el: usize = r.galois_element((1 << i) >> 1, i == 0, log_nth_root);
215
216        if i == 0{
217            r.a_mul_b_scalar_barrett_into_c::<ONCE>(&self.n_inv, b, tmp_b);
218            r.a_mul_b_montgomery_into_a::<ONCE>(&self.x_pow_2[r.log_n()-i-1], tmp_b);
219        }else{
220            r.a_mul_b_montgomery_into_c::<ONCE>(b, &self.x_pow_2[r.log_n()-i-1], tmp_b);
221        }
222
223        r.a_apply_automorphism_into_b::<true>(tmp_b, gal_el, nth_root, &mut self.buf_auto);
224        r.a_sub_b_into_a::<1, ONCE>(&self.buf_auto, tmp_b);
225
226        tmp_b
227    }
228}