he_ring/number_ring/hypercube/
interpolate.rs

1
2use feanor_math::algorithms::extension_invert::invert_over_local_zn;
3use feanor_math::algorithms::linsolve::LinSolveRing;
4use feanor_math::homomorphism::CanIsoFromTo;
5use feanor_math::homomorphism::Homomorphism;
6use feanor_math::homomorphism::SelfIso;
7use feanor_math::local::PrincipalLocalRing;
8use feanor_math::ring::*;
9use feanor_math::rings::extension::extension_impl::FreeAlgebraImpl;
10use feanor_math::rings::extension::*;
11use feanor_math::rings::field::AsFieldBase;
12use feanor_math::rings::poly::*;
13use feanor_math::rings::zn::FromModulusCreateableZnRing;
14use feanor_math::rings::zn::ZnReductionMap;
15use feanor_math::rings::zn::ZnRing;
16
17///
18/// Interpolation data for a list of moduli `f1, ..., fn` that can be used
19/// to derive from remainders `r1, ..., rn` an "interpolation polynomial" `h`
20/// such that `h = ri mod fi`.
21/// 
22/// Clearly this requires that the moduli `fi` are pairwise coprime. Additionally,
23/// we currently require that all interpolation unit vectors `ei` (defined
24/// by `ei = 1 mod fi`, `ei = 0 mod fj` for `j != i`) exist over the base
25/// ring (e.g. they might not be integral, even if the base ring is `Z`).
26/// 
27pub struct FastPolyInterpolation<P>
28    where P: RingStore,
29        P::Type: PolyRing,
30        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing
31{
32    poly_ring: P,
33    input_degree: usize,
34    unit_vectors: Vec<Vec<(El<P>, El<P>)>>,
35    final_modulus: El<P>,
36    n: usize
37}
38
39///
40/// Computes a polynomial `h` of degree `< deg(fg)` such that `h = 1 mod f` and `h = 0 mod g`.
41/// 
42fn crt_unit_vectors<P>(poly_ring: P, f: &El<P>, g: &El<P>) -> El<P>
43    where P: RingStore,
44        P::Type: PolyRing,
45        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing + FromModulusCreateableZnRing + ZnRing + PrincipalLocalRing,
46        AsFieldBase<RingValue<<<P::Type as RingExtension>::BaseRing as RingStore>::Type>>: CanIsoFromTo<<<P::Type as RingExtension>::BaseRing as RingStore>::Type> + SelfIso
47{
48    assert!(poly_ring.base_ring().is_one(poly_ring.lc(f).unwrap()));
49    assert!(poly_ring.base_ring().is_one(poly_ring.lc(g).unwrap()));
50    let deg_f = poly_ring.degree(&f).unwrap();
51
52    let mod_f_ring = FreeAlgebraImpl::new(poly_ring.base_ring(), deg_f, (0..deg_f).map(|i| poly_ring.base_ring().negate(poly_ring.base_ring().clone_el(poly_ring.coefficient_at(&f, i)))).collect::<Vec<_>>());
53    let g_mod_f = poly_ring.div_rem_monic(poly_ring.clone_el(&g), &f).1;
54    let g_mod_f = mod_f_ring.from_canonical_basis((0..deg_f).map(|i| poly_ring.base_ring().clone_el(poly_ring.coefficient_at(&g_mod_f, i))));
55
56    let normalization_factor = invert_over_local_zn(RingRef::new(mod_f_ring.get_ring()), &g_mod_f);
57    
58    assert!(normalization_factor.is_some(), "crt unit vector modulo {} and {} does not exist", poly_ring.format(f), poly_ring.format(g));
59    debug_assert!(mod_f_ring.is_one(&mod_f_ring.mul_ref(normalization_factor.as_ref().unwrap(), &g_mod_f)));
60    let g_mod_f_inv = mod_f_ring.poly_repr(&poly_ring, &normalization_factor.unwrap(), poly_ring.base_ring().identity());
61    let first_unit_vector = poly_ring.mul_ref_snd(g_mod_f_inv, &g);
62
63    debug_assert!(poly_ring.is_one(&poly_ring.div_rem_monic(poly_ring.clone_el(&first_unit_vector), &f).1));
64    debug_assert!(poly_ring.is_zero(&poly_ring.div_rem_monic(poly_ring.clone_el(&first_unit_vector), &g).1));
65
66    return first_unit_vector;
67}
68
69impl<P> FastPolyInterpolation<P>
70    where P: RingStore,
71        P::Type: PolyRing,
72        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing + FromModulusCreateableZnRing + ZnRing + PrincipalLocalRing,
73        AsFieldBase<RingValue<<<P::Type as RingExtension>::BaseRing as RingStore>::Type>>: CanIsoFromTo<<<P::Type as RingExtension>::BaseRing as RingStore>::Type> + SelfIso
74{
75    pub fn new(poly_ring: P, moduli: Vec<El<P>>) -> Self {
76        let n = moduli.len();
77        let input_degree = moduli.iter().map(|f| poly_ring.degree(f).unwrap_or(0)).max().unwrap();
78        let mut current = moduli.iter().map(|f| poly_ring.clone_el(f)).collect::<Vec<_>>();
79        let mut result = Vec::new();
80        while current.len() > 1 {
81            let mut result_part = Vec::new();
82            let mut new = Vec::new();
83            let mut moduli_it = current.array_chunks::<2>();
84            for [f, g] in moduli_it.by_ref() {
85                let new_modulus = poly_ring.mul_ref(f, g);
86                let unit_vectors = (crt_unit_vectors(&poly_ring, f, g), crt_unit_vectors(&poly_ring, g, f));
87                result_part.push(unit_vectors);
88                new.push(new_modulus);
89            }
90            if let Some(last) = moduli_it.remainder().get(0) {
91                new.push(poly_ring.clone_el(last));
92            }
93            current = new;
94            result.push(result_part);
95        }
96        return Self {
97            final_modulus: current.pop().unwrap(),
98            n: n,
99            input_degree: input_degree,
100            poly_ring: poly_ring,
101            unit_vectors: result
102        };
103    }
104
105    pub fn change_modulus<PNew>(&self, new_poly_ring: PNew) -> FastPolyInterpolation<PNew>
106        where PNew: RingStore,
107            PNew::Type: PolyRing,
108            <<PNew::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing + FromModulusCreateableZnRing + ZnRing + PrincipalLocalRing,
109            AsFieldBase<RingValue<<<PNew::Type as RingExtension>::BaseRing as RingStore>::Type>>: CanIsoFromTo<<<PNew::Type as RingExtension>::BaseRing as RingStore>::Type> + SelfIso
110    {
111        let red_map = ZnReductionMap::new(self.poly_ring().base_ring(), new_poly_ring.base_ring()).unwrap();
112        let lifted_red_map = new_poly_ring.lifted_hom(self.poly_ring(), &red_map);
113        let unit_vectors = self.unit_vectors.iter().map(|list| list.iter().map(|(e0, e1)| (lifted_red_map.map_ref(e0), lifted_red_map.map_ref(e1))).collect()).collect();
114        let final_modulus = lifted_red_map.map_ref(&self.final_modulus);
115        return FastPolyInterpolation {
116            final_modulus: final_modulus,
117            input_degree: self.input_degree,
118            poly_ring: new_poly_ring,
119            n: self.n,
120            unit_vectors: unit_vectors
121        };
122    }
123
124    pub fn poly_ring(&self) -> &P {
125        &self.poly_ring
126    }
127
128    pub fn final_modulus(&self) -> &El<P> {
129        &self.final_modulus
130    }
131
132    ///
133    /// Computes a polynomial of degree `< 2 * deg(prod(moduli))` that is congruent
134    /// to `remainders[i]` modulo `moduli[i]`.
135    /// 
136    /// It is unreduced, since we can reduce its degree to `< deg(prod(moduli))` by
137    /// taking the remainder modulo `prod(moduli)`.
138    /// 
139    /// However, this can be computed really fast, in time `n log(n)^2` if FFT-based
140    /// polynomial multiplication is used by the underlying polynomial ring. It is also
141    /// very fast in practice, since we don't perform any polynomial division.
142    /// 
143    pub fn interpolate_unreduced(&self, remainders: Vec<El<P>>) -> El<P> {
144        assert_eq!(self.n, remainders.len());
145        for i in 0..self.n {
146            assert!(self.poly_ring.degree(&remainders[i]).unwrap_or(0) < self.input_degree);
147        }
148        let mut current = remainders;
149        for i in 0..self.unit_vectors.len() {
150            let mut new = Vec::new();
151            let mut current_it = current.array_chunks::<2>();
152            for ([f0, f1], (e0, e1)) in current_it.by_ref().zip(self.unit_vectors[i].iter()) {
153                new.push(self.poly_ring.add(
154                    self.poly_ring.mul_ref(f0, e0),
155                    self.poly_ring.mul_ref(f1, e1)
156                ));
157            }
158            if let Some(last) = current_it.remainder().get(0) {
159                new.push(self.poly_ring.clone_el(last));
160            }
161            current = new;
162        }
163        debug_assert_eq!(1, current.len());
164        let result = current.pop().unwrap();
165        debug_assert!(self.poly_ring.degree(&result).unwrap_or(0) < (self.input_degree << (self.unit_vectors.len() +  1)));
166        return result;
167    }
168}
169
170#[cfg(test)]
171use feanor_math::assert_el_eq;
172#[cfg(test)]
173use feanor_math::rings::local::AsLocalPIR;
174#[cfg(test)]
175use feanor_math::rings::poly::dense_poly::DensePolyRing;
176#[cfg(test)]
177use feanor_math::rings::zn::zn_64::Zn;
178
179#[test]
180fn test_interpolate() {
181    let base_ring = AsLocalPIR::from_zn(Zn::new(257)).unwrap();
182    let poly_ring = DensePolyRing::new(base_ring, "X");
183
184    let moduli = poly_ring.with_wrapped_indeterminate(|X| [
185        X.pow_ref(2) - 1,
186        X.pow_ref(3) + X - 1,
187        X.pow_ref(2) - X + 2,
188        X.pow_ref(2) - 2 * X
189    ]);
190    let remainders = poly_ring.with_wrapped_indeterminate(|X| [
191        -1 * X + 3,
192        -5 * X.pow_ref(2) + 21 * X - 12,
193        -728 * X + 16,
194        720896 * X
195    ]);
196    let interpolation = FastPolyInterpolation::new(&poly_ring, moduli.iter().map(|f| poly_ring.clone_el(f)).collect());
197    let result = interpolation.interpolate_unreduced(remainders.iter().map(|f| poly_ring.clone_el(f)).collect());
198    for i in 0..4 {
199        assert_el_eq!(&poly_ring, &remainders[i], poly_ring.div_rem_monic(poly_ring.clone_el(&result), &moduli[i]).1);
200    }
201
202    let moduli = poly_ring.with_wrapped_indeterminate(|X| [
203        X.pow_ref(2) - 1,
204        X.pow_ref(3) + X - 1,
205        X.pow_ref(2) - X + 2
206    ]);
207    let remainders = poly_ring.with_wrapped_indeterminate(|X| [
208        -1 * X + 3,
209        -5 * X.pow_ref(2) + 21 * X - 12,
210        -728 * X + 16
211    ]);
212    let interpolation = FastPolyInterpolation::new(&poly_ring, moduli.iter().map(|f| poly_ring.clone_el(f)).collect());
213    let result = interpolation.interpolate_unreduced(remainders.iter().map(|f| poly_ring.clone_el(f)).collect());
214    for i in 0..3 {
215        assert_el_eq!(&poly_ring, &remainders[i], poly_ring.div_rem_monic(poly_ring.clone_el(&result), &moduli[i]).1);
216    }
217    
218    let moduli = poly_ring.with_wrapped_indeterminate(|X| [
219        X.pow_ref(2) - 1,
220        X.pow_ref(2) - 2,
221        X.pow_ref(2) - 3,
222        X.pow_ref(2) - 4,
223        X.pow_ref(2) - 5,
224        X.pow_ref(2) - 6,
225        X.pow_ref(2) - 7,
226        X.pow_ref(2) - 8,
227        X.pow_ref(2) - 9,
228        X.pow_ref(2) - 10,
229    ]);
230    let remainders = poly_ring.with_wrapped_indeterminate(|X| [
231        5 * X + 10,
232        5 * X + 9,
233        5 * X + 8,
234        5 * X + 7,
235        5 * X + 6,
236        5 * X + 5,
237        5 * X + 4,
238        5 * X + 3,
239        5 * X + 2,
240        5 * X + 1,
241    ]);
242    let interpolation = FastPolyInterpolation::new(&poly_ring, moduli.iter().map(|f| poly_ring.clone_el(f)).collect());
243    let result = interpolation.interpolate_unreduced(remainders.iter().map(|f| poly_ring.clone_el(f)).collect());
244    for i in 0..10 {
245        assert_el_eq!(&poly_ring, &remainders[i], poly_ring.div_rem_monic(poly_ring.clone_el(&result), &moduli[i]).1);
246    }
247}
he_ring/number_ring/hypercube/interpolate.rs

he_ring/number_ring/hypercube/
interpolate.rs
