feanor_math/algorithms/
interpolate.rs

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

feanor_math/algorithms/
interpolate.rs
