poulpy/math/src/ring/impl_u64/packing.rs
math/src/ring/impl_u64/packing.rsRaw
1use std::collections::HashMap;
2use crate::poly::Poly;
3use crate::ring::Ring;
4use crate::modulus::{ONCE, WordOps};
5use crate::modulus::barrett::Barrett;
6use std::cmp::min;
7use std::mem::transmute;
8
9
10impl Ring<u64>{
11
12 // Generates a vector storing {X^{2^0}, X^{2^1}, .., X^{2^log_n}}.
13 pub fn gen_x_pow_2<const NTT: bool, const INV: bool>(&self, log_n: usize) -> Vec<Poly<u64>>{
14 let mut x_pow: Vec<Poly<u64>> = Vec::<Poly<u64>>::with_capacity(log_n);
15
16 (0..log_n).for_each(|i|{
17 let mut idx: usize = 1<<i;
18
19 if INV{
20 idx = self.n() - idx;
21 }
22
23 x_pow.push(self.new_poly());
24
25 if i == 0{
26 x_pow[i].0[idx] = self.modulus.montgomery.one();
27 self.ntt_inplace::<false>(&mut x_pow[i]);
28 }else{
29 let (left, right) = x_pow.split_at_mut(i);
30 self.a_mul_b_montgomery_into_c::<ONCE>(&left[i-1], &left[i-1], &mut right[0]);
31 }
32 });
33
34 if INV{
35 self.a_neg_into_a::<1, ONCE>(&mut x_pow[0]);
36 }
37
38 if !NTT{
39 x_pow.iter_mut().for_each(|x| self.intt_inplace::<false>(x));
40 }
41
42 x_pow
43 }
44
45 pub fn pack<'a, const ZEROGARBAGE: bool, const NTT: bool>(&self, polys: &'a mut HashMap<usize, &'a mut Poly<u64>>, log_gap: usize) -> &'a Poly<u64>{
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 keys: Vec<usize> = polys.keys().copied().collect();
54 keys.sort();
55
56 let mut gap = 0usize;
57
58 if keys.len() > 1{
59 gap = max_pow2_gap(&keys);
60 }else{
61 gap = 1<<log_n;
62 }
63
64 let log_gap: usize = gap.log2();
65
66 if !ZEROGARBAGE{
67 if gap > 0 {
68 log_end -= log_gap;
69 }
70 }
71
72 let n_inv: Barrett<u64> = self.modulus.barrett.prepare(self.modulus.inv(1<<(log_end - log_start)));
73
74 for (_, poly) in polys.iter_mut() {
75 if !NTT{
76 self.ntt_inplace::<true>(poly);
77 }
78
79 self.a_mul_b_scalar_barrett_into_a::<ONCE>(&n_inv, poly);
80 }
81
82 let x_pow2: Vec<Poly<u64>> = self.gen_x_pow_2::<true, false>(log_n);
83 let mut tmpa: Poly<u64> = self.new_poly();
84 let mut tmpb: Poly<u64> = self.new_poly();
85
86 for i in log_start..log_end{
87
88 let t: usize = 1<<(log_n-1-i);
89
90 for j in 0..t{
91
92 let option_lo: Option<&&mut Poly<u64>> = polys.get(&i);
93 let option_hi: Option<&&mut Poly<u64>> = polys.get(&(i+t));
94 let mut hi_exists: bool = false;
95
96 match option_hi{
97 Some(hi) =>{
98
99 // Unsafe code is necessary because two mutables references are
100 // accessed from the map.
101 unsafe{
102 self.a_mul_b_montgomery_into_a::<ONCE>(&x_pow2[log_n-i-1], transmute(*hi as *const Poly<u64> as *mut Poly<u64>));
103 }
104
105 hi_exists = true;
106
107 match option_lo{
108 Some(lo) =>{
109
110 self.a_sub_b_into_c::<1, ONCE>(lo, hi, &mut tmpa);
111
112 // Ensures unsafe blocks are "safe".
113 let ptr_hi: *mut Poly<u64> = *hi as *const Poly<u64> as *mut Poly<u64>;
114 let ptr_lo: *mut Poly<u64> = *lo as *const Poly<u64> as *mut Poly<u64>;
115 assert!(ptr_hi != ptr_lo, "something went seriously wrong");
116
117 unsafe{
118 self.a_add_b_into_b::<ONCE>(hi, transmute(ptr_lo));
119 }
120 }
121 None =>{
122 unsafe{
123 polys.insert(j, transmute(*hi as *const Poly<u64> as *mut Poly<u64>));
124 }
125 },
126 }
127
128 polys.remove(&(j+t));
129 }
130
131 None =>{},
132 }
133
134 let option_lo: Option<&&mut Poly<u64>> = polys.get(&i);
135 let option_hi: Option<&&mut Poly<u64>> = polys.get(&(i+t));
136
137 match option_lo{
138 Some(lo) =>{
139
140 let gal_el: usize = self.galois_element(1<<(i-1), i == 0, log_nth_root);
141
142 if hi_exists{
143 self.automorphism::<true>(&tmpa, gal_el, 2<<self.log_n(), &mut tmpb);
144 }else{
145 self.automorphism::<true>(*lo, gal_el, nth_root, &mut tmpa);
146 }
147 unsafe{
148 self.a_add_b_into_b::<ONCE>(&tmpa, transmute(*lo as *const Poly<u64> as *mut Poly<u64>));
149 }
150 }
151
152 None =>{
153 match option_hi{
154 Some(hi) =>{
155 let gal_el: usize = self.galois_element(1<<(i-1), i == 0, log_nth_root);
156
157 self.automorphism::<true>(*hi, gal_el, nth_root, &mut tmpa);
158
159 unsafe{
160 self.a_sub_b_into_a::<1, ONCE>(&tmpa, transmute(*hi as *const Poly<u64> as *mut Poly<u64>))
161 }
162 }
163
164 None =>{}
165 }
166 }
167 }
168 }
169 }
170
171 *polys.get(&0).unwrap()
172 }
173}
174
175
176// Returns the largest
177fn max_pow2_gap(vec: &[usize]) -> usize{
178 let mut gap: usize = usize::MAX;
179 for i in 1..vec.len(){
180 let (l, r) = (vec[i-1], vec[i]);
181 assert!(l > r, "invalid input vec: not sorted");
182 gap = min(gap, r-l);
183 if gap == 1{
184 break;
185 }
186 };
187 1 << gap.trailing_zeros()
188}