feanor_math/algorithms/
interpolate.rs

1use crate::divisibility::{DivisibilityRing, DivisibilityRingStore, Domain};
2use crate::integer::IntegerRingStore;
3use crate::iters::multi_cartesian_product;
4use crate::primitive_int::StaticRing;
5use crate::ring::*;
6use crate::rings::multivariate::*;
7use crate::rings::poly::*;
8use crate::seq::*;
9use crate::homomorphism::Homomorphism;
10use crate::rings::poly::dense_poly::DensePolyRing;
11
12use std::alloc::Allocator;
13use std::cmp::min;
14use std::ops::Range;
15
16///
17/// Computes `out[i] = prod_(j != i) values[j]`.
18/// 
19/// This algorithm recursively halfes the input, and thus it implicitly pads the input to the
20/// next power-of-two. The time complexity is then `O(n sum_(0 <= i <= log n) T(2^i))`
21/// where `T(d)` is the complexity of multiplying two products of `d` input elements. If the cost
22/// of multiplication is constant, this becomes `O(n log n T)`. In the other common case, where
23/// the input elements are degree-1 polynomials, this becomes `O(n sum_(0 <= i <= log n) 2^(i c)) = O(n^(c + 1))`
24/// where `c` is the multiplication exponent, e.g. `c = 1.58...` for Karatsuba multiplication.
25/// 
26/// # Example
27/// ```rust
28/// # use feanor_math::ring::*;
29/// # use feanor_math::primitive_int::*;
30/// # use feanor_math::algorithms::interpolate::*;
31/// # use feanor_math::seq::*;
32/// let ring = StaticRing::<i64>::RING;
33/// let mut result = [0; 6];
34/// product_except_one(ring, (1..7).map_fn(|x| x as i64), &mut result);
35/// let factorial_6 = 6 * 5 * 4 * 3 * 2 * 1;
36/// // `product_except_one()` computes exactly these values, but without using any divisions
37/// let expected = [factorial_6 / 1, factorial_6 / 2, factorial_6 / 3, factorial_6 / 4, factorial_6 / 5, factorial_6 / 6];
38/// assert_eq!(expected, result);
39/// ```
40/// 
41#[stability::unstable(feature = "enable")]
42pub fn product_except_one<V, R>(ring: R, values: V, out: &mut [El<R>])
43    where R: RingStore,
44        V: VectorFn<El<R>>
45{
46    assert_eq!(values.len(), out.len());
47    let n = values.len();
48    let log2_n = StaticRing::<i64>::RING.abs_log2_ceil(&n.try_into().unwrap()).unwrap();
49    assert!(n <= (1 << log2_n));
50    if n % 2 == 0 {
51        for i in 0..n {
52            out[i] = values.at(i ^ 1);
53        }
54    } else {
55        for i in 0..(n - 1) {
56            out[i] = values.at(i ^ 1);
57        }
58        out[n - 1] = ring.one();
59    }
60    for s in 1..log2_n {
61        for j in 0..(1 << (log2_n - s - 1)) {
62            let block_index = j << (s + 1);
63            if block_index + (1 << s) < n {
64                let (fst, snd) = (&mut out[block_index..min(n, block_index + (1 << (s + 1)))]).split_at_mut(1 << s);
65                let snd_block_prod = ring.mul_ref_fst(&snd[0], values.at(block_index + (1 << s)));
66                let fst_block_prod = ring.mul_ref_fst(&fst[0], values.at(block_index));
67                for i in 0..(1 << s) {
68                    ring.mul_assign_ref(&mut fst[i], &snd_block_prod);
69                }
70                for i in 0..snd.len() {
71                    ring.mul_assign_ref(&mut snd[i], &fst_block_prod);
72                }
73            }
74        }
75    }
76}
77
78#[stability::unstable(feature = "enable")]
79#[derive(PartialEq, Eq, Hash, Debug, Clone, Copy)]
80pub enum InterpolationError {
81    NotInvertible
82}
83
84///
85/// Uses Lagrange interpolation to compute the interpolation polynomial of the given values.
86/// Concretely, this is the univariate polynomial `f` of degree `< x.len()` such that `f(x[i]) = y[i]`
87/// for all `i`.
88/// 
89/// If no such polynomial exists (this is only possible if the base ring is not a field), an 
90/// error is returned.
91/// 
92/// The complexity is `O(n T(n))` where `T(n)` is the cost of multiplying two degree-`n` polynomials.
93/// 
94/// # Example
95/// ```rust
96/// # #![feature(allocator_api)]
97/// # use std::alloc::Global;
98/// # use feanor_math::ring::*;
99/// # use feanor_math::seq::*;
100/// # use feanor_math::assert_el_eq;
101/// # use feanor_math::algorithms::interpolate::*;
102/// # use feanor_math::rings::poly::*;
103/// # use feanor_math::primitive_int::*;
104/// # use feanor_math::rings::poly::dense_poly::*;
105/// let ZZX = DensePolyRing::new(StaticRing::<i64>::RING, "X");
106/// let [expected] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) + 1]);
107/// let actual = interpolate(&ZZX, [1, 2, 6].copy_els(), [2, 5, 37].copy_els(), Global).unwrap();
108/// assert_el_eq!(&ZZX, expected, actual);
109/// ```
110/// In some cases the interpolation polynomial does not exist.
111/// ```rust
112/// # #![feature(allocator_api)]
113/// # use std::alloc::Global;
114/// # use feanor_math::ring::*;
115/// # use feanor_math::primitive_int::StaticRing;
116/// # use feanor_math::homomorphism::Homomorphism;
117/// # use feanor_math::seq::*;
118/// # use feanor_math::algorithms::interpolate::*;
119/// # use feanor_math::rings::poly::*;
120/// # use feanor_math::rings::poly::dense_poly::*;
121/// let ZnX = DensePolyRing::new(StaticRing::<i64>::RING, "X");
122/// let actual = interpolate(&ZnX, [-2, 0, 2].copy_els(), [1, 0, 1].copy_els(), Global);
123/// assert!(actual.is_err());
124/// ```
125/// 
126#[stability::unstable(feature = "enable")]
127pub fn interpolate<P, V1, V2, A: Allocator>(poly_ring: P, x: V1, y: V2, allocator: A) -> Result<El<P>, InterpolationError>
128    where P: RingStore,
129        P::Type: PolyRing,
130        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing + Domain,
131        V1: VectorFn<El<<P::Type as RingExtension>::BaseRing>>,
132        V2: VectorFn<El<<P::Type as RingExtension>::BaseRing>>
133{
134    assert_eq!(x.len(), y.len());
135    let R = poly_ring.base_ring();
136    let null_poly = poly_ring.prod(x.iter().map(|x| poly_ring.from_terms([(R.one(), 1), (R.negate(x), 0)])));
137    let mut nums = Vec::with_capacity_in(x.len(), &allocator);
138    let div_linear = |poly: &El<P>, a: &El<<P::Type as RingExtension>::BaseRing>| if let Some(d) = poly_ring.degree(poly) {
139        poly_ring.from_terms((0..d).rev().scan(R.zero(), |current, i| {
140            R.add_assign_ref(current, poly_ring.coefficient_at(poly, i + 1));
141            let result = R.clone_el(current);
142            R.mul_assign_ref(current, a);
143            return Some((result, i));
144        }))
145    } else { poly_ring.zero() };
146    nums.extend(x.iter().map(|x| div_linear(&null_poly, &x)));
147
148    let mut denoms = Vec::with_capacity_in(x.len(), &allocator);
149    denoms.extend((0..x.len()).map(|i| poly_ring.evaluate(&nums[i], &x.at(i), &R.identity())));
150    let mut factors = Vec::with_capacity_in(x.len(), &allocator);
151    factors.resize_with(x.len(), || R.zero());
152    product_except_one(R, (&denoms[..]).into_clone_ring_els(R), &mut factors);
153    let denominator = R.mul_ref(&factors[0], &denoms[0]);
154    for i in 0..x.len() {
155        R.mul_assign(&mut factors[i], y.at(i));
156    }
157
158    if let Some(inv) = R.invert(&denominator) {
159        Ok(poly_ring.inclusion().mul_map(<_ as RingStore>::sum(&poly_ring, nums.into_iter().zip(factors.into_iter()).map(|(num, c)| poly_ring.inclusion().mul_map(num, c))), inv))
160    } else {
161        let scaled_result = <_ as RingStore>::sum(&poly_ring, nums.into_iter().zip(factors.into_iter()).map(|(num, c)| poly_ring.inclusion().mul_map(num, c)));
162        poly_ring.try_from_terms(poly_ring.terms(&scaled_result).map(|(c, i)| R.checked_div(&c, &denominator).map(|c| (c, i)).ok_or(InterpolationError::NotInvertible)))
163    }
164}
165
166#[stability::unstable(feature = "enable")]
167pub fn interpolate_multivariate<P, V1, V2, A, A2>(poly_ring: P, interpolation_points: V1, mut values: Vec<El<<P::Type as RingExtension>::BaseRing>, A2>, allocator: A) -> Result<El<P>, InterpolationError>
168    where P: RingStore,
169        P::Type: MultivariatePolyRing,
170        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing + Domain,
171        V1: VectorFn<V2>,
172        V2: VectorFn<El<<P::Type as RingExtension>::BaseRing>>,
173        A: Allocator,
174        A2: Allocator
175{
176    let dim_prod = |range: Range<usize>| <_ as RingStore>::prod(&StaticRing::<i64>::RING, range.map(|i| interpolation_points.at(i).len().try_into().unwrap())) as usize;
177    assert_eq!(interpolation_points.len(), poly_ring.indeterminate_count());
178    let n = poly_ring.indeterminate_count();
179    assert_eq!(values.len(), dim_prod(0..n));
180
181    let uni_poly_ring = DensePolyRing::new_with_convolution(poly_ring.base_ring(), "X", &allocator, STANDARD_CONVOLUTION);
182
183    for i in (0..n).rev() {
184        let leading_dim = dim_prod((i + 1)..n);
185        let outer_block_count = dim_prod(0..i);
186        let len = interpolation_points.at(i).len();
187        let outer_block_size = leading_dim * len;
188        for outer_block_index in 0..outer_block_count {
189            for inner_block_index in 0..leading_dim {
190                let block_start = inner_block_index + outer_block_index * outer_block_size;
191                let poly = interpolate(&uni_poly_ring, interpolation_points.at(i), (&values[..]).into_clone_ring_els(poly_ring.base_ring()).restrict(block_start..(block_start + outer_block_size + 1 - leading_dim)).step_by_fn(leading_dim), &allocator)?;
192                for j in 0..len {
193                    values[block_start + leading_dim * j] = poly_ring.base_ring().clone_el(uni_poly_ring.coefficient_at(&poly, j));
194                }
195            }
196        }
197    }
198    return Ok(poly_ring.from_terms(
199        multi_cartesian_product((0..n).map(|i| 0..interpolation_points.at(i).len()), |idxs| poly_ring.get_ring().create_monomial(idxs.iter().map(|e| *e)), |_, x| *x)
200            .zip(values.into_iter())
201            .map(|(m, c)| (c, m))
202    ));
203}
204
205#[cfg(test)]
206use crate::rings::zn::zn_64::Zn;
207#[cfg(test)]
208use crate::rings::zn::ZnRingStore;
209#[cfg(test)]
210use std::alloc::Global;
211#[cfg(test)]
212use multivariate_impl::MultivariatePolyRingImpl;
213
214use super::convolution::STANDARD_CONVOLUTION;
215
216#[test]
217fn test_product_except_one() {
218    let ring = StaticRing::<i64>::RING;
219    let data = [2, 3, 5, 7, 11, 13, 17, 19];
220    let mut actual = [0; 8];
221    let expected = [
222        3 * 5 * 7 * 11 * 13 * 17 * 19,
223        2 * 5 * 7 * 11 * 13 * 17 * 19,
224        2 * 3 * 7 * 11 * 13 * 17 * 19,
225        2 * 3 * 5 * 11 * 13 * 17 * 19,
226        2 * 3 * 5 * 7 * 13 * 17 * 19,
227        2 * 3 * 5 * 7 * 11 * 17 * 19,
228        2 * 3 * 5 * 7 * 11 * 13 * 19,
229        2 * 3 * 5 * 7 * 11 * 13 * 17
230    ];
231    product_except_one(&ring, (&data[..]).clone_els_by(|x| *x), &mut actual);
232    assert_eq!(expected, actual);
233
234    let data = [2, 3, 5, 7, 11, 13, 17];
235    let mut actual = [0; 7];
236    let expected = [
237        3 * 5 * 7 * 11 * 13 * 17,
238        2 * 5 * 7 * 11 * 13 * 17,
239        2 * 3 * 7 * 11 * 13 * 17,
240        2 * 3 * 5 * 11 * 13 * 17,
241        2 * 3 * 5 * 7 * 13 * 17,
242        2 * 3 * 5 * 7 * 11 * 17,
243        2 * 3 * 5 * 7 * 11 * 13
244    ];
245    product_except_one(&ring, (&data[..]).clone_els_by(|x| *x), &mut actual);
246    assert_eq!(expected, actual);
247
248    let data = [2, 3, 5, 7, 11, 13];
249    let mut actual = [0; 6];
250    let expected = [
251        3 * 5 * 7 * 11 * 13,
252        2 * 5 * 7 * 11 * 13,
253        2 * 3 * 7 * 11 * 13,
254        2 * 3 * 5 * 11 * 13,
255        2 * 3 * 5 * 7 * 13,
256        2 * 3 * 5 * 7 * 11
257    ];
258    product_except_one(&ring, (&data[..]).clone_els_by(|x| *x), &mut actual);
259    assert_eq!(expected, actual);
260}
261
262#[test]
263fn test_interpolate() {
264    let ring = StaticRing::<i64>::RING;
265    let poly_ring = DensePolyRing::new(ring, "X");
266    let poly = poly_ring.from_terms([(3, 0), (1, 1), (-1, 3), (2, 4), (1, 5)].into_iter());
267    let x = (0..6).map_fn(|x| x.try_into().unwrap());
268    let actual = interpolate(&poly_ring, x.clone(), x.map_fn(|x| poly_ring.evaluate(&poly, &x, &ring.identity())), Global).unwrap();
269    assert_el_eq!(&poly_ring, &poly, &actual);
270
271    let ring = Zn::new(29).as_field().ok().unwrap();
272    let poly_ring = DensePolyRing::new(ring, "X");
273    let x = (0..5).map_fn(|x| ring.int_hom().map(x as i32));
274    let y = (0..5).map_fn(|x| if x == 3 { ring.int_hom().map(6) } else { ring.zero() });
275    let poly = interpolate(&poly_ring, x.clone(), y.clone(), Global).unwrap();
276    for i in 0..5 {
277        assert_el_eq!(ring, y.at(i), poly_ring.evaluate(&poly, &x.at(i), &ring.identity()));
278    }
279}
280
281#[test]
282fn test_interpolate_multivariate() {
283    let ring = Zn::new(29).as_field().ok().unwrap();
284    let poly_ring: MultivariatePolyRingImpl<_> = MultivariatePolyRingImpl::new(ring, 2);
285
286    let interpolation_points = (0..2).map_fn(|_| (0..5).map_fn(|x| ring.int_hom().map(x as i32)));
287    let values = (0..25).map(|x| ring.int_hom().map(x & 1)).collect::<Vec<_>>();
288    let poly = interpolate_multivariate(&poly_ring, &interpolation_points, values, Global).unwrap();
289
290    for x in 0..5 {
291        for y in 0..5 {
292            let expected = (x * 5 + y) & 1;
293            assert_el_eq!(ring, ring.int_hom().map(expected), poly_ring.evaluate(&poly, [ring.int_hom().map(x), ring.int_hom().map(y)].into_clone_ring_els(&ring), &ring.identity()));
294        }
295    }
296
297    let poly_ring: MultivariatePolyRingImpl<_> = MultivariatePolyRingImpl::new(ring, 3);
298
299    let interpolation_points = (0..3).map_fn(|i| (0..(i + 2)).map_fn(|x| ring.int_hom().map(x as i32)));
300    let values = (0..24).map(|x| ring.int_hom().map(x / 2)).collect::<Vec<_>>();
301    let poly = interpolate_multivariate(&poly_ring, &interpolation_points, values, Global).unwrap();
302
303    for x in 0..2 {
304        for y in 0..3 {
305            for z in 0..4 {
306                let expected = (x * 12 + y * 4 + z) / 2;
307                assert_el_eq!(ring, ring.int_hom().map(expected), poly_ring.evaluate(&poly, [ring.int_hom().map(x), ring.int_hom().map(y), ring.int_hom().map(z)].into_clone_ring_els(&ring), &ring.identity()));
308            }
309        }
310    }
311}
312
313#[test]
314#[ignore]
315fn large_polynomial_interpolation() {
316    let field = Zn::new(65537).as_field().ok().unwrap();
317    let poly_ring = DensePolyRing::new(field, "X");
318    let hom = poly_ring.base_ring().can_hom(&StaticRing::<i64>::RING).unwrap();
319    let actual = interpolate(&poly_ring, (0..65536).map_fn(|x| hom.map(x as i64)), (0..65536).map_fn(|x| hom.map(x as i64)), Global).unwrap();
320    assert_el_eq!(&poly_ring, poly_ring.indeterminate(), actual);
321}
feanor_math/algorithms/interpolate.rs

feanor_math/algorithms/
interpolate.rs
