feanor_math/reduce_lift/
poly_eval.rs

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

feanor_math/reduce_lift/
poly_eval.rs
