feanor_math/reduce_lift/
poly_eval.rs

1use std::fmt::Debug;
2
3use crate::algorithms::linsolve::LinSolveRing;
4use crate::algorithms::resultant::ComputeResultantRing;
5use crate::divisibility::{DivisibilityRing, Domain};
6use crate::field::*;
7use crate::homomorphism::*;
8use crate::pid::PrincipalIdealRing;
9use crate::ring::*;
10use crate::rings::finite::FiniteRing;
11use crate::specialization::FiniteRingSpecializable;
12
13/// Trait for rings that can be temporarily replaced by an extension when we need more points,
14/// e.g. for interpolation.
15///
16/// Note that a trivial implementation is possible for every ring of characteristic 0, since
17/// these already have infinitely many points whose pairwise differences are non-zero-divisors.
18/// Such an implementation can be added to new types using the macro
19/// [`impl_interpolation_base_ring_char_zero!`].
20#[stability::unstable(feature = "enable")]
21pub trait InterpolationBaseRing: DivisibilityRing {
22    /// The type of the extension ring we can switch to to get more points.
23    ///
24    /// For the reason why there are so many quite specific trait bounds here:
25    /// See the doc of [`EvalPolyLocallyRing::LocalRingBase`].
26    type ExtendedRingBase<'a>: ?Sized + PrincipalIdealRing + Domain + LinSolveRing + ComputeResultantRing + Debug
27    where
28        Self: 'a;
29
30    type ExtendedRing<'a>: RingStore<Type = Self::ExtendedRingBase<'a>> + Clone
31    where
32        Self: 'a;
33
34    fn in_base<'a, S>(&self, ext_ring: S, el: El<S>) -> Option<Self::Element>
35    where
36        Self: 'a,
37        S: RingStore<Type = Self::ExtendedRingBase<'a>>;
38
39    fn in_extension<'a, S>(&self, ext_ring: S, el: Self::Element) -> El<S>
40    where
41        Self: 'a,
42        S: RingStore<Type = Self::ExtendedRingBase<'a>>;
43
44    /// Returns `count` points such that the difference between any two of them
45    /// is a non-zero-divisor.
46    ///
47    /// Any two calls must give elements in the same order.
48    fn interpolation_points<'a>(&'a self, count: usize) -> (Self::ExtendedRing<'a>, Vec<El<Self::ExtendedRing<'a>>>);
49}
50
51/// [`RingStore`] for [`InterpolationBaseRing`].
52#[stability::unstable(feature = "enable")]
53pub trait InterpolationBaseRingStore: RingStore
54where
55    Self::Type: InterpolationBaseRing,
56{
57}
58
59impl<R> InterpolationBaseRingStore for R
60where
61    R: RingStore,
62    R::Type: InterpolationBaseRing,
63{
64}
65
66/// The inclusion map `R -> S` for a ring `R` and one of its extensions `S`
67/// as given by [`InterpolationBaseRing`].
68#[stability::unstable(feature = "enable")]
69pub struct ToExtRingMap<'a, R>
70where
71    R: ?Sized + InterpolationBaseRing,
72{
73    ring: RingRef<'a, R>,
74    ext_ring: R::ExtendedRing<'a>,
75}
76
77impl<'a, R> ToExtRingMap<'a, R>
78where
79    R: ?Sized + InterpolationBaseRing,
80{
81    #[stability::unstable(feature = "enable")]
82    pub fn for_interpolation(ring: &'a R, point_count: usize) -> (Self, Vec<El<R::ExtendedRing<'a>>>) {
83        let (ext_ring, points) = ring.interpolation_points(point_count);
84        return (
85            Self {
86                ring: RingRef::new(ring),
87                ext_ring,
88            },
89            points,
90        );
91    }
92
93    #[stability::unstable(feature = "enable")]
94    pub fn as_base_ring_el(&self, el: El<R::ExtendedRing<'a>>) -> R::Element {
95        self.ring.get_ring().in_base(&self.ext_ring, el).unwrap()
96    }
97}
98
99impl<'a, R> Homomorphism<R, R::ExtendedRingBase<'a>> for ToExtRingMap<'a, R>
100where
101    R: ?Sized + InterpolationBaseRing,
102{
103    type CodomainStore = R::ExtendedRing<'a>;
104    type DomainStore = RingRef<'a, R>;
105
106    fn codomain(&self) -> &Self::CodomainStore { &self.ext_ring }
107
108    fn domain(&self) -> &Self::DomainStore { &self.ring }
109
110    fn map(&self, x: <R as RingBase>::Element) -> <R::ExtendedRingBase<'a> as RingBase>::Element {
111        self.ring.get_ring().in_extension(&self.ext_ring, x)
112    }
113}
114
115/// Trait for rings that support performing computations locally.
116///
117/// Note that here (and in `feanor-math` generally), the term "local" is used to refer to algorithms
118/// that work modulo prime ideals (or their powers), which is different from the mathematical
119/// concept of localization.
120///
121/// More concretely, a ring `R` implementing this trait should be endowed with a "pseudo norm"
122/// ```text
123///   |.|: R  ->  [0, ∞)
124/// ```
125/// i.e. a symmetric, sub-additive, sub-multiplicative map.
126/// Furthermore, for any bound `b`, the ring should be able to provide prime ideals
127/// `p1, ..., pk` together with the rings `Ri = R / pi`, such that the restricted
128/// reduction map
129/// ```text
130///   { x in R | |x| <= b }  ->  R1 x ... x Rk
131/// ```
132/// is injective.
133/// This means that a computation can be performed in the simpler ring `R1 x ... x Rk`,
134/// and - assuming the result is of pseudo-norm `<= b`, mapped back to `R`.
135///
136/// The standard use case is the evaluation of a multivariate polynomial `f(X1, ..., Xm)`
137/// over this ring. The trait is designed to enable the following approach:
138///  - Given ring elements `a1, ..., am`, compute an upper bound `B` on `|f(a1, ..., am)|`. The
139///    values `|ai|` are given by [`EvalPolyLocallyRing::pseudo_norm()`].
140///  - Get a sufficient number of prime ideals, using [`EvalPolyLocallyRing::local_computation()`]
141///  - Compute `f(a1 mod pi, ..., am mod pi) mod pi` for each prime `pi` within the ring given by
142///    [`EvalPolyLocallyRing::local_ring_at()`]. The reductions `ai mod pj` are given by
143///    [`EvalPolyLocallyRing::reduce()`].
144///  - Recombine the results to an element of `R` by using [`EvalPolyLocallyRing::lift_combine()`].
145///
146/// # Relationship with [`crate::reduce_lift::poly_factor_gcd::PolyGCDLocallyDomain`]
147///
148/// There are generally two ways of computing something via a reduce-modulo-primes-then-lift
149/// approach. Either one can take many different prime ideals, or one can take a large power
150/// of a single/a small amount of prime ideals.
151///
152/// This trait is for the former approach, which is especially suitable when the computation to
153/// perform can be written as a polynomial evaluation. In particular, this applicable to
154/// determinants, resultant, and (with some caveats) solving linear systems.
155///
156/// On the other hand, when factoring polynomials or computing their gcds, it is common to instead
157/// rely on Hensel lifting to compute the result modulo a large power of a single prime, or very
158/// few primes. This approach is formalized by
159/// [`crate::reduce_lift::poly_factor_gcd::PolyGCDLocallyDomain`].
160///
161/// # Type-level recursion in feanor-math
162///
163/// [`EvalPolyLocallyRing`] and [`PolyGCDLocallyDomain`] are the two traits in feanor-math which
164/// use type-level recursion with blanket implementations. The idea is simple: If our ring is an
165/// [`EvalPolyLocallyRing`], whose quotients are again [`EvalPolyLocallyRing`], whose quotients are
166/// again [`EvalPolyLocallyRing`] and so on, ending with a finite field. Then, for all kinds of
167/// operations which are actually just polynomial evaluations (resultants, determinants,
168/// unique-solution linear systems, ...) we already have all the information for an efficient
169/// algorithm. However, providing this through blanket implementations is not that simple. We now
170/// explain how it is done in feanor-math.
171///
172/// The proper solution:
173/// ```compile_fail,E0391
174/// #![feature(min_specialization)]
175/// // Some special property, in practice this might be LinSolveRing, ComputeResultantRing, FactorPolyRing, ...
176/// pub trait SpecialProperty {
177///
178///     fn foo<'a>(&'a self);
179/// }
180///
181/// impl<T> SpecialProperty for T
182///     where T: Recursive,
183///         for<'a> T::PrevCase<'a>: SpecialProperty
184/// {
185///     default fn foo<'a>(&'a self) {
186///         println!("recursion");
187///         <T::PrevCase<'a> as SpecialProperty>::foo(&self.map_down());
188///     }
189/// }
190///
191/// // Newtype to allow `&'a T: SpecialProperty` without conflicting with above blanket impl
192/// pub struct RefWrapper<'a, T>(&'a T);
193///
194/// impl<'a, T: Recursive> SpecialProperty for RefWrapper<'a, T> {
195///
196///     fn foo<'b>(&'b self) {
197///         self.0.foo()
198///     }
199/// }
200///
201/// // The main recursion type
202/// pub trait Recursive {
203///
204///     type PrevCase<'a> where Self: 'a;
205///
206///     fn map_down<'a>(&'a self) -> Self::PrevCase<'a>;
207/// }
208///
209/// #[derive(Clone, Copy)]
210/// pub struct BaseCase;
211///
212/// impl Recursive for BaseCase {
213///
214///     type PrevCase<'a> = Self where Self: 'a;
215///
216///     fn map_down<'a>(&'a self) -> Self::PrevCase<'a> { *self }
217/// }
218///
219/// impl SpecialProperty for BaseCase {
220///
221///     fn foo<'a>(&'a self) {
222///         println!("BaseCaseSpecial")
223///     }
224/// }
225///
226/// pub struct RecursiveCase<T: Recursive>(T);
227///
228/// impl<T: Recursive> Recursive for RecursiveCase<T> {
229///
230///     type PrevCase<'a> = RefWrapper<'a, T> where Self: 'a;
231///
232///     fn map_down<'a>(&'a self) -> Self::PrevCase<'a> {
233///         RefWrapper(&self.0)
234///     }
235/// }
236///
237/// fn run_type_recursion<'b, T>(t: &'b T)
238///     where T: 'b + Recursive,
239///         for<'a> T::PrevCase<'a>: SpecialProperty
240/// {
241///     t.foo()
242/// }
243///
244/// run_type_recursion(&BaseCase);
245/// run_type_recursion(&RecursiveCase(BaseCase));
246/// ```
247/// Unfortunately, this currently fails with an error
248/// ```text
249/// error[E0391]: cycle detected when computing whether impls specialize one another
250///   --> src/main.rs:9:1
251///    |
252/// 9  | / impl<T> SpecialProperty for T
253/// 10 | |     where T: Recursive,
254/// 11 | |         for<'a> T::PrevCase<'a>: SpecialProperty
255///    | |________________________________________________^
256///    |
257///    = note: ...which requires evaluating trait selection obligation `BaseCase: SpecialProperty`...
258///    = note: ...which again requires computing whether impls specialize one another, completing the cycle
259/// ```
260/// The only solution I found is to put the constraint by `SpecialProperty` directly into
261/// the declaration `type PrevCase<'a>: SpecialProperty where Self: 'a`. This works, but
262/// unfortunately makes things much less generic than we would want to. In particular,
263/// If we have a `NonSpecialBaseCase` in addition to `BaseCase`, and we want it to not
264/// implement `SpecialProperty`, then we cannot use `RecursiveCase` as well. On the other hand,
265/// if the compiler would accept the first variant, this would work out fine - the only effect
266/// of `NonSpecialBaseCase: !SpecialProperty` would then be that also
267/// `RecursiveCase<NonSpecialBaseCase>: !SpecialProperty`, which of course makes sense.
268///
269/// However, for now we stay with this workaround, which then looks as follows:
270/// ```
271/// #![feature(min_specialization)]
272/// // Some special property, in practice this might be LinSolveRing, ComputeResultantRing, FactorPolyRing, ...
273/// pub trait SpecialProperty {
274///
275///     fn foo<'a>(&'a self);
276/// }
277///
278/// impl<T> SpecialProperty for T
279///     where T: Recursive
280/// {
281///     default fn foo<'a>(&'a self) {
282///         println!("recursion");
283///         <T::PrevCase<'a> as SpecialProperty>::foo(&self.map_down());
284///     }
285/// }
286///
287/// // Newtype to allow `&'a T: SpecialProperty` without conflicting with above blanket impl
288/// pub struct RefWrapper<'a, T: SpecialProperty>(&'a T);
289///
290/// impl<'a, T: Recursive> SpecialProperty for RefWrapper<'a, T> {
291///
292///     fn foo<'b>(&'b self) {
293///         self.0.foo()
294///     }
295/// }
296///
297/// // The main recursion type
298/// pub trait Recursive {
299///
300///     type PrevCase<'a>: SpecialProperty where Self: 'a;
301///
302///     fn map_down<'a>(&'a self) -> Self::PrevCase<'a>;
303/// }
304///
305/// #[derive(Clone, Copy)]
306/// pub struct BaseCase;
307///
308/// impl Recursive for BaseCase {
309///
310///     type PrevCase<'a> = Self where Self: 'a;
311///
312///     fn map_down<'a>(&'a self) -> Self::PrevCase<'a> { *self }
313/// }
314///
315/// impl SpecialProperty for BaseCase {
316///
317///     fn foo<'a>(&'a self) {
318///         println!("BaseCaseSpecial")
319///     }
320/// }
321///
322/// pub struct RecursiveCase<T: Recursive + SpecialProperty>(T);
323///
324/// impl<T: Recursive> Recursive for RecursiveCase<T> {
325///
326///     type PrevCase<'a> = RefWrapper<'a, T> where Self: 'a;
327///
328///     fn map_down<'a>(&'a self) -> Self::PrevCase<'a> {
329///         RefWrapper(&self.0)
330///     }
331/// }
332///
333/// fn run_type_recursion<'b, T>(t: &'b T)
334///     where T: 'b + Recursive
335/// {
336///     t.foo()
337/// }
338///
339/// run_type_recursion(&BaseCase);
340/// run_type_recursion(&RecursiveCase(BaseCase));
341/// ```
342///
343/// Another advantage is that every function
344#[stability::unstable(feature = "enable")]
345pub trait EvalPolyLocallyRing: RingBase + FiniteRingSpecializable {
346    /// The type of the ring we get once quotienting by a prime ideal.
347    ///
348    /// The proper solution would be to just require `LocalRingBase<'ring>: ?Sized + RingBase`.
349    /// Algorithms which require this type to provide additional functionality (like being
350    /// a [`ComputeResultantRing`]) should then add a constraint
351    /// ```ignore
352    ///   for<'ring> R::LocalRingBase<'ring>: ComputeResultantRing
353    /// ```
354    /// However, this causes problems, more concretely an error
355    /// "cycle detected when computing whether impls specialize one another".
356    /// More on that in the trait-level doc.
357    ///
358    /// Hence, we have to already bound `LocalRingBase` by all the traits that are reasonably
359    /// required by some algorithms. This clearly makes it a lot less generic than it should be,
360    /// but so far it worked out without too much trouble.
361    type LocalRingBase<'ring>: ?Sized + PrincipalIdealRing + Domain + LinSolveRing + ComputeResultantRing + Debug
362    where
363        Self: 'ring;
364
365    type LocalRing<'ring>: RingStore<Type = Self::LocalRingBase<'ring>>
366    where
367        Self: 'ring;
368
369    /// A collection of prime ideals of the ring, and additionally any data required to reconstruct
370    /// a small ring element from its projections onto each prime ideal.
371    type LocalComputationData<'ring>
372    where
373        Self: 'ring;
374
375    /// Computes (an upper bound of) the natural logarithm of the pseudo norm of a ring element.
376    ///
377    /// The pseudo norm should be
378    ///  - symmetric, i.e. `|-x| = |x|`,
379    ///  - sub-additive, i.e. `|x + y| <= |x| + |y|`
380    ///  - sub-multiplicative, i.e. `|xy| <= |x| |y|`
381    ///
382    /// and this function should give `ln|x|`
383    fn ln_pseudo_norm(&self, el: &Self::Element) -> f64;
384
385    /// Sets up the context for a new polynomial evaluation, whose output
386    /// should have pseudo norm less than the given bound.
387    fn local_computation<'ring>(&'ring self, ln_pseudo_norm_bound: f64) -> Self::LocalComputationData<'ring>;
388
389    /// Returns the number `k` of local rings that are required
390    /// to get the correct result of the given computation.
391    fn local_ring_count<'ring>(&self, computation: &Self::LocalComputationData<'ring>) -> usize
392    where
393        Self: 'ring;
394
395    /// Returns the `i`-th local ring belonging to the given computation.
396    fn local_ring_at<'ring>(&self, computation: &Self::LocalComputationData<'ring>, i: usize) -> Self::LocalRing<'ring>
397    where
398        Self: 'ring;
399
400    /// Computes the map `R -> R1 x ... x Rk`, i.e. maps the given element into each of
401    /// the local rings.
402    fn reduce<'ring>(
403        &self,
404        computation: &Self::LocalComputationData<'ring>,
405        el: &Self::Element,
406    ) -> Vec<<Self::LocalRingBase<'ring> as RingBase>::Element>
407    where
408        Self: 'ring;
409
410    /// Computes a preimage under the map `R -> R1 x ... x Rk`, i.e. a ring element `x` that reduces
411    /// to each of the given local rings under the map [`EvalPolyLocallyRing::reduce()`].
412    ///
413    /// The result should have pseudo-norm bounded by the bound given when the computation
414    /// was initialized, via [`EvalPolyLocallyRing::local_computation()`].
415    fn lift_combine<'ring>(
416        &self,
417        computation: &Self::LocalComputationData<'ring>,
418        el: &[<Self::LocalRingBase<'ring> as RingBase>::Element],
419    ) -> Self::Element
420    where
421        Self: 'ring;
422}
423
424impl<R> EvalPolyLocallyRing for R
425where
426    R: ?Sized + FiniteRing + Field + Debug + SelfIso,
427{
428    type LocalComputationData<'ring>
429        = RingRef<'ring, Self>
430    where
431        Self: 'ring;
432
433    type LocalRing<'ring>
434        = RingRef<'ring, Self>
435    where
436        Self: 'ring;
437
438    type LocalRingBase<'ring>
439        = Self
440    where
441        Self: 'ring;
442
443    fn ln_pseudo_norm(&self, _el: &Self::Element) -> f64 { 0.0 }
444
445    fn local_computation<'ring>(&'ring self, _ln_pseudo_norm_bound: f64) -> Self::LocalComputationData<'ring> {
446        RingRef::new(self)
447    }
448
449    fn local_ring_at<'ring>(&self, computation: &Self::LocalComputationData<'ring>, _i: usize) -> Self::LocalRing<'ring>
450    where
451        Self: 'ring,
452    {
453        *computation
454    }
455
456    fn local_ring_count<'ring>(&self, _computation: &Self::LocalComputationData<'ring>) -> usize
457    where
458        Self: 'ring,
459    {
460        1
461    }
462
463    fn reduce<'ring>(
464        &self,
465        _computation: &Self::LocalComputationData<'ring>,
466        el: &Self::Element,
467    ) -> Vec<<Self::LocalRingBase<'ring> as RingBase>::Element>
468    where
469        Self: 'ring,
470    {
471        vec![self.clone_el(el)]
472    }
473
474    fn lift_combine<'ring>(
475        &self,
476        _computation: &Self::LocalComputationData<'ring>,
477        el: &[<Self::LocalRingBase<'ring> as RingBase>::Element],
478    ) -> Self::Element
479    where
480        Self: 'ring,
481    {
482        assert_eq!(1, el.len());
483        return self.clone_el(&el[0]);
484    }
485}
486
487/// The map `R -> R/p` for a ring `R` and one of its local quotients `R/p` as
488/// given by [`EvalPolyLocallyRing`].
489#[stability::unstable(feature = "enable")]
490pub struct EvaluatePolyLocallyReductionMap<'ring, 'data, R>
491where
492    R: 'ring + ?Sized + EvalPolyLocallyRing,
493    'ring: 'data,
494{
495    ring: RingRef<'data, R>,
496    data: &'data R::LocalComputationData<'ring>,
497    local_ring: R::LocalRing<'ring>,
498    index: usize,
499}
500
501impl<'ring, 'data, R> EvaluatePolyLocallyReductionMap<'ring, 'data, R>
502where
503    R: 'ring + ?Sized + EvalPolyLocallyRing,
504    'ring: 'data,
505{
506    #[stability::unstable(feature = "enable")]
507    pub fn new(ring: &'data R, data: &'data R::LocalComputationData<'ring>, index: usize) -> Self {
508        Self {
509            ring: RingRef::new(ring),
510            data,
511            local_ring: ring.local_ring_at(data, index),
512            index,
513        }
514    }
515}
516
517impl<'ring, 'data, R> Homomorphism<R, R::LocalRingBase<'ring>> for EvaluatePolyLocallyReductionMap<'ring, 'data, R>
518where
519    R: 'ring + ?Sized + EvalPolyLocallyRing,
520    'ring: 'data,
521{
522    type CodomainStore = R::LocalRing<'ring>;
523    type DomainStore = RingRef<'data, R>;
524
525    fn codomain(&self) -> &Self::CodomainStore { &self.local_ring }
526
527    fn domain(&self) -> &Self::DomainStore { &self.ring }
528
529    fn map(&self, x: <R as RingBase>::Element) -> <R::LocalRingBase<'ring> as RingBase>::Element {
530        let ring_ref: &'data R = self.ring.into();
531        let mut reductions: Vec<<R::LocalRingBase<'ring> as RingBase>::Element> = ring_ref.reduce(self.data, &x);
532        return reductions.swap_remove(self.index);
533    }
534}
535
536/// Generates an implementation of [`crate::reduce_lift::poly_eval::InterpolationBaseRing`]
537/// for a ring of characteristic zero. Not that in this case, the only sensible implementation
538/// is trivial, since the ring itself has enough elements for any interpolation task.
539///
540/// # Example
541/// ```rust
542/// # use std::fmt::Debug;
543/// # use feanor_math::ring::*;
544/// # use feanor_math::delegate::*;
545/// # use feanor_math::reduce_lift::poly_eval::*;
546/// # use feanor_math::primitive_int::*;
547/// # use feanor_math::algorithms::resultant::*;
548/// # use feanor_math::divisibility::*;
549/// # use feanor_math::homomorphism::Homomorphism;
550/// # use feanor_math::rings::poly::*;
551/// # use feanor_math::rings::poly::dense_poly::DensePolyRing;
552/// # use feanor_math::pid::*;
553/// # use feanor_math::impl_interpolation_base_ring_char_zero;
554/// // we wrap a `RingBase` here for simplicity, but in practice a wrapper should
555/// // always store a `RingStore` instead
556/// #[derive(PartialEq, Debug)]
557/// struct MyRingWrapper<R: RingBase>(R);
558/// impl<R: RingBase> DelegateRing for MyRingWrapper<R> {
559///     type Element = R::Element;
560///     type Base = R;
561///     fn get_delegate(&self) -> &Self::Base { &self.0 }
562///     fn delegate(&self, x: R::Element) -> R::Element { x }
563///     fn rev_delegate(&self, x: R::Element) -> R::Element { x }
564///     fn delegate_ref<'a>(&self, x: &'a R::Element) -> &'a R::Element { x }
565///     fn delegate_mut<'a>(&self, x: &'a mut R::Element) -> &'a mut R::Element { x }
566/// }
567/// impl<R: RingBase> DelegateRingImplEuclideanRing for MyRingWrapper<R> {}
568/// impl<R: RingBase> DelegateRingImplFiniteRing for MyRingWrapper<R> {}
569/// impl<R: Domain> Domain for MyRingWrapper<R> {}
570/// impl_interpolation_base_ring_char_zero!{ <{ R }> InterpolationBaseRing for MyRingWrapper<R> where R: PrincipalIdealRing + Domain + ComputeResultantRing + Debug }
571///
572/// // now we can use `InterpolationBaseRing`-functionality
573/// let ring = MyRingWrapper(StaticRing::<i64>::RING.into());
574/// let (embedding, points) = ToExtRingMap::for_interpolation(&ring, 3);
575/// assert_eq!(0, points[0]);
576/// assert_eq!(1, points[1]);
577/// assert_eq!(2, points[2]);
578///
579/// // There is a problem here, described in EvalPolyLocallyRing::LocalRingBase.
580/// // Short version: we need to manually impl ComputeResultantRing
581/// impl<R: ComputeResultantRing> ComputeResultantRing for MyRingWrapper<R> {
582///     fn resultant<P>(poly_ring: P, f: El<P>, g: El<P>) -> Self::Element
583///         where P: RingStore + Copy,
584///             P::Type: PolyRing,
585///             <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>
586///     {
587///         let new_poly_ring = DensePolyRing::new(RingRef::new(poly_ring.base_ring().get_ring().get_delegate()), "X");
588///         let hom = new_poly_ring.lifted_hom(&poly_ring, UnwrapHom::new(poly_ring.base_ring(), new_poly_ring.base_ring()));
589///         poly_ring.base_ring().get_ring().rev_delegate(R::resultant(&new_poly_ring, hom.map(f), hom.map(g)))
590///     }
591/// }
592/// ```
593#[macro_export]
594macro_rules! impl_interpolation_base_ring_char_zero {
595    (<{$($gen_args:tt)*}> InterpolationBaseRing for $self_type:ty where $($constraints:tt)*) => {
596        impl<$($gen_args)*> $crate::reduce_lift::poly_eval::InterpolationBaseRing for $self_type where $($constraints)* {
597
598            type ExtendedRing<'a> = RingRef<'a, Self>
599                where Self: 'a;
600
601            type ExtendedRingBase<'a> = Self
602                where Self: 'a;
603
604            fn in_base<'a, S>(&self, _ext_ring: S, el: El<S>) -> Option<Self::Element>
605                where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>
606            {
607                Some(el)
608            }
609
610            fn in_extension<'a, S>(&self, _ext_ring: S, el: Self::Element) -> El<S>
611                where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>
612            {
613                el
614            }
615
616            fn interpolation_points<'a>(&'a self, count: usize) -> (Self::ExtendedRing<'a>, Vec<El<Self::ExtendedRing<'a>>>) {
617                let ZZbig = $crate::integer::BigIntRing::RING;
618                assert!(ZZbig.is_zero(&self.characteristic(&ZZbig).unwrap()));
619                let ring = $crate::ring::RingRef::new(self);
620                (ring, (0..count).map(|n| <_ as $crate::homomorphism::Homomorphism<_, _>>::map(&ring.int_hom(), n.try_into().unwrap())).collect())
621            }
622        }
623    };
624    (InterpolationBaseRing for $self_type:ty) => {
625        impl_interpolation_base_ring_char_zero!{ <{}> InterpolationBaseRing for $self_type where }
626    }
627}
628
629/// Implements [`EvalPolyLocallyRing`] for an integer ring.
630///
631/// This uses a default implementation, where the prime ideals are given by the largest prime
632/// numbers such that the corresponding residue field can be implemented using
633/// [`crate::rings::zn::zn_64::Zn`]. This should be suitable in almost all scenarios.
634///
635/// The syntax is the same as for other impl-macros, see e.g.
636/// [`crate::impl_interpolation_base_ring_char_zero!`].
637#[macro_export]
638macro_rules! impl_eval_poly_locally_for_ZZ {
639    (EvalPolyLocallyRing for $int_ring_type:ty) => {
640        impl_eval_poly_locally_for_ZZ!{ <{}> EvalPolyLocallyRing for $int_ring_type where }
641    };
642    (<{$($gen_args:tt)*}> EvalPolyLocallyRing for $int_ring_type:ty where $($constraints:tt)*) => {
643
644        impl<$($gen_args)*> $crate::reduce_lift::poly_eval::EvalPolyLocallyRing for $int_ring_type
645            where $($constraints)*
646        {
647            type LocalComputationData<'ring> = $crate::rings::zn::zn_rns::Zn<$crate::rings::field::AsField<$crate::rings::zn::zn_64::Zn>, RingRef<'ring, Self>>
648                where Self: 'ring;
649
650            type LocalRing<'ring> = $crate::rings::field::AsField<$crate::rings::zn::zn_64::Zn>
651                where Self: 'ring;
652
653            type LocalRingBase<'ring> = $crate::rings::field::AsFieldBase<$crate::rings::zn::zn_64::Zn>
654                where Self: 'ring;
655
656            fn ln_pseudo_norm(&self, el: &Self::Element) -> f64 {
657                RingRef::new(self).abs_log2_ceil(el).unwrap_or(0) as f64 * 2f64.ln()
658            }
659
660            fn local_computation<'ring>(&'ring self, ln_pseudo_norm_bound: f64) -> Self::LocalComputationData<'ring> {
661                let mut primes = Vec::new();
662                let mut ln_current = 0.;
663                let mut prime_it = //$crate::reduce_lift::primelist::LARGE_PRIMES.iter().copied().chain
664                ((0..).scan((1 << 62) / 9, |current, _| {
665                    *current = $crate::algorithms::miller_rabin::prev_prime(StaticRing::<i64>::RING, *current).unwrap();
666                    if *current < (1 << 32) {
667                        panic!("not enough primes");
668                    }
669                    return Some($crate::rings::zn::zn_64::Zn::new(*current as u64));
670                }));
671                while ln_current < ln_pseudo_norm_bound + 1. {
672                    let Fp = prime_it.next().unwrap();
673                    ln_current += (*$crate::rings::zn::ZnRingStore::modulus(&Fp) as f64).ln();
674                    primes.push(Fp);
675                }
676                return $crate::rings::zn::zn_rns::Zn::new(
677                    primes.into_iter().map(|Fp| $crate::rings::field::AsField::from($crate::rings::field::AsFieldBase::promise_is_perfect_field(Fp))).collect(),
678                    RingRef::new(self)
679                );
680            }
681
682            fn local_ring_at<'ring>(&self, computation: &Self::LocalComputationData<'ring>, i: usize) -> Self::LocalRing<'ring>
683                where Self: 'ring
684            {
685                <_ as $crate::seq::VectorView<_>>::at(computation, i).clone()
686            }
687
688            fn local_ring_count<'ring>(&self, computation: &Self::LocalComputationData<'ring>) -> usize
689                where Self: 'ring
690            {
691                <_ as $crate::seq::VectorView<_>>::len(computation)
692            }
693
694            fn reduce<'ring>(&self, computation: &Self::LocalComputationData<'ring>, el: &Self::Element) -> Vec<<Self::LocalRingBase<'ring> as RingBase>::Element>
695                where Self: 'ring
696            {
697                <_ as $crate::seq::VectorView<_>>::as_iter(&computation.get_congruence(&computation.coerce(RingValue::from_ref(self), self.clone_el(el)))).map(|x| *x).collect()
698            }
699
700            fn lift_combine<'ring>(&self, computation: &Self::LocalComputationData<'ring>, el: &[<Self::LocalRingBase<'ring> as RingBase>::Element]) -> Self::Element
701                where Self: 'ring
702            {
703                <_ as $crate::rings::zn::ZnRingStore>::smallest_lift(computation, computation.from_congruence(el.iter().copied()))
704            }
705        }
706    };
707}
feanor_math/reduce_lift/poly_eval.rs

feanor_math/reduce_lift/
poly_eval.rs
