feanor_math/
compute_locally.rs

1
2use crate::algorithms::miller_rabin::prev_prime;
3use crate::divisibility::{DivisibilityRing, Domain};
4use crate::homomorphism::*;
5use crate::integer::{IntegerRing, IntegerRingStore};
6use crate::pid::PrincipalIdealRing;
7use crate::primitive_int::StaticRing;
8use crate::ring::*;
9use crate::rings::field::{AsField, AsFieldBase};
10use crate::rings::zn::{zn_64, zn_rns, ZnRingStore};
11use crate::seq::VectorView;
12
13///
14/// Trait for rings that can be temporarily replaced by an extension when we need more points,
15/// e.g. for interpolation.
16/// 
17/// Note that a trivial implementation is possible for every ring of characteristic 0, since
18/// these already have infinitely many points whose pairwise differences are non-zero-divisors. 
19/// Such an implementation can be added to new types using the macro [`impl_interpolation_base_ring_char_zero!`].
20/// 
21#[stability::unstable(feature = "enable")]
22pub trait InterpolationBaseRing: DivisibilityRing {
23
24    ///
25    /// Restricting this here to be `DivisibilityRing + PrincipalIdealRing + Domain`
26    /// is necessary, because of a compiler bug, see also
27    /// [`crate::compute_locally::EvaluatePolyLocallyRing::LocalRingBase`]
28    /// 
29    type ExtendedRingBase<'a>: ?Sized + DivisibilityRing + PrincipalIdealRing + Domain
30        where Self: 'a;
31
32    type ExtendedRing<'a>: RingStore<Type = Self::ExtendedRingBase<'a>> + Clone
33        where Self: 'a;
34
35    fn in_base<'a, S>(&self, ext_ring: S, el: El<S>) -> Option<Self::Element>
36        where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>;
37
38    fn in_extension<'a, S>(&self, ext_ring: S, el: Self::Element) -> El<S>
39        where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>;
40
41    ///
42    /// Returns `count` points such that the difference between any two of them
43    /// is a non-zero-divisor.
44    /// 
45    /// Any two calls must give elements in the same order.
46    /// 
47    fn interpolation_points<'a>(&'a self, count: usize) -> (Self::ExtendedRing<'a>, Vec<El<Self::ExtendedRing<'a>>>);
48}
49
50///
51/// [`RingStore`] for [`InterpolationBaseRing`].
52/// 
53#[stability::unstable(feature = "enable")]
54pub trait InterpolationBaseRingStore: RingStore
55    where Self::Type: InterpolationBaseRing
56{}
57
58impl<R> InterpolationBaseRingStore for R
59    where R: RingStore, R::Type: InterpolationBaseRing
60{}
61
62///
63/// The inclusion map `R -> S` for a ring `R` and one of its extensions `S`
64/// as given by [`InterpolationBaseRing`].
65/// 
66#[stability::unstable(feature = "enable")]
67pub struct ToExtRingMap<'a, R>
68    where R: ?Sized + InterpolationBaseRing
69{
70    ring: RingRef<'a, R>,
71    ext_ring: R::ExtendedRing<'a>
72}
73
74impl<'a, R> ToExtRingMap<'a, R>
75    where R: ?Sized + InterpolationBaseRing
76{
77    #[stability::unstable(feature = "enable")]
78    pub fn for_interpolation(ring: &'a R, point_count: usize) -> (Self, Vec<El<R::ExtendedRing<'a>>>) {
79        let (ext_ring, points) = ring.interpolation_points(point_count);
80        return (Self { ring: RingRef::new(ring), ext_ring: ext_ring }, points);
81    }
82
83    #[stability::unstable(feature = "enable")]
84    pub fn as_base_ring_el(&self, el: El<R::ExtendedRing<'a>>) -> R::Element {
85        self.ring.get_ring().in_base(&self.ext_ring, el).unwrap()
86    }
87}
88
89impl<'a, R> Homomorphism<R, R::ExtendedRingBase<'a>> for ToExtRingMap<'a, R>
90    where R: ?Sized + InterpolationBaseRing
91{
92    type CodomainStore = R::ExtendedRing<'a>;
93    type DomainStore = RingRef<'a, R>;
94
95    fn codomain<'b>(&'b self) -> &'b Self::CodomainStore {
96        &self.ext_ring
97    }
98
99    fn domain<'b>(&'b self) -> &'b Self::DomainStore {
100        &self.ring
101    }
102
103    fn map(&self, x: <R as RingBase>::Element) -> <R::ExtendedRingBase<'a> as RingBase>::Element {
104        self.ring.get_ring().in_extension(&self.ext_ring, x)
105    }
106}
107
108///
109/// Trait for rings that support performing computations locally. 
110/// 
111/// Note that here (and in `feanor-math` generally), the term "local" is used to refer to algorithms
112/// that work modulo prime ideals (or their powers), which is different from the mathematical concept
113/// of localization.
114/// 
115/// More concretely, a ring `R` implementing this trait should be endowed with a "pseudo norm"
116/// ```text
117///   |.|: R  ->  [0, ∞)
118/// ```
119/// i.e. a symmetric, sub-additive, sub-multiplicative map.
120/// Furthermore, for any bound `b`, the ring should be able to provide prime ideals
121/// `p1, ..., pk` together with the rings `Ri = R / pi`, such that the restricted
122/// reduction map
123/// ```text
124///   { x in R | |x| <= b }  ->  R1 x ... x Rk
125/// ```
126/// is injective.
127/// This means that a computation can be performed in the simpler ring `R1 x ... x Rk`,
128/// and - assuming the result is of pseudo-norm `<= b`, mapped back to `R`.
129/// 
130/// The standard use case is the evaluation of a multivariate polynomial `f(X1, ..., Xm)`
131/// over this ring. The trait is designed to enable the following approach:
132///  - Given ring elements `a1, ..., am`, compute an upper bound `B` on `|f(a1, ..., am)|`.
133///    The values `|ai|` are given by [`EvaluatePolyLocallyRing::pseudo_norm()`].
134///  - Get a sufficient number of prime ideals, using [`EvaluatePolyLocallyRing::local_computation()`] 
135///  - Compute `f(a1 mod pi, ..., am mod pi) mod pi` for each prime `pi` within the ring given by 
136///    [`EvaluatePolyLocallyRing::local_ring_at()`]. The reductions `ai mod pj` are given by
137///    [`EvaluatePolyLocallyRing::reduce()`].
138///  - Recombine the results to an element of `R` by using [`EvaluatePolyLocallyRing::lift_combine()`].
139/// 
140#[stability::unstable(feature = "enable")]
141pub trait EvaluatePolyLocallyRing: RingBase {
142    
143    ///
144    /// The proper way would be to define this with two lifetime parameters `'ring` and `'data`,
145    /// to allow it to reference both the ring itself and the current `LocalComputationData`.
146    /// However, when doing this, I ran into the compiler bug
147    /// [https://github.com/rust-lang/rust/issues/100013](https://github.com/rust-lang/rust/issues/100013).
148    /// 
149    /// This is also the reason why we restrict this type here to be [`PrincipalIdealRing`], because
150    /// unfortunately, the a constraint `for<'a> SomeRing::LocalRingBase<'a>: PrincipalIdealRing` again
151    /// triggers the bug in some settings.
152    /// 
153    type LocalRingBase<'ring>: ?Sized + PrincipalIdealRing + Domain
154        where Self: 'ring;
155
156    type LocalRing<'ring>: RingStore<Type = Self::LocalRingBase<'ring>>
157        where Self: 'ring;
158
159    type LocalComputationData<'ring>
160        where Self: 'ring;
161
162    ///
163    /// Computes (an upper bound of) the natural logarithm of the pseudo norm of a ring element.
164    /// 
165    /// The pseudo norm should be
166    ///  - symmetric, i.e. `|-x| = |x|`,
167    ///  - sub-additive, i.e. `|x + y| <= |x| + |y|`
168    ///  - sub-multiplicative, i.e. `|xy| <= |x| |y|`
169    /// and this function should give `ln|x|`
170    /// 
171    fn ln_pseudo_norm(&self, el: &Self::Element) -> f64;
172
173    ///
174    /// Sets up the context for a new polynomial evaluation, whose output
175    /// should have pseudo norm less than the given bound.
176    /// 
177    fn local_computation<'ring>(&'ring self, ln_pseudo_norm_bound: f64) -> Self::LocalComputationData<'ring>;
178
179    ///
180    /// Returns the number `k` of local rings that are required
181    /// to get the correct result of the given computation.
182    /// 
183    fn local_ring_count<'ring>(&self, computation: &Self::LocalComputationData<'ring>) -> usize
184        where Self: 'ring;
185
186    ///
187    /// Returns the `i`-th local ring belonging to the given computation.
188    /// 
189    fn local_ring_at<'ring>(&self, computation: &Self::LocalComputationData<'ring>, i: usize) -> Self::LocalRing<'ring>
190        where Self: 'ring;
191
192    ///
193    /// Computes the map `R -> R1 x ... x Rk`, i.e. maps the given element into each of
194    /// the local rings.
195    /// 
196    fn reduce<'ring>(&self, computation: &Self::LocalComputationData<'ring>, el: &Self::Element) -> Vec<<Self::LocalRingBase<'ring> as RingBase>::Element>
197        where Self: 'ring;
198
199    ///
200    /// Computes a preimage under the map `R -> R1 x ... x Rk`, i.e. a ring element `x` that reduces
201    /// to each of the given local rings under the map [`EvaluatePolyLocallyRing::reduce()`].
202    /// 
203    /// The result should have pseudo-norm bounded by the bound given when the computation
204    /// was initialized, via [`EvaluatePolyLocallyRing::local_computation()`].
205    /// 
206    fn lift_combine<'ring>(&self, computation: &Self::LocalComputationData<'ring>, el: &[<Self::LocalRingBase<'ring> as RingBase>::Element]) -> Self::Element
207        where Self: 'ring;
208}
209
210impl<I> EvaluatePolyLocallyRing for I
211    where I: IntegerRing
212{
213    type LocalComputationData<'ring> = zn_rns::Zn<AsField<zn_64::Zn>, RingRef<'ring, Self>>
214        where Self: 'ring;
215
216    type LocalRing<'ring> = AsField<zn_64::Zn>
217        where Self: 'ring;
218
219    type LocalRingBase<'ring> = AsFieldBase<zn_64::Zn>
220        where Self: 'ring;
221
222    fn ln_pseudo_norm(&self, el: &Self::Element) -> f64 {
223        RingRef::new(self).abs_log2_ceil(el).unwrap_or(0) as f64 * 2f64.ln()
224    }
225
226    fn local_computation<'ring>(&'ring self, ln_pseudo_norm_bound: f64) -> Self::LocalComputationData<'ring> {
227        let mut primes = Vec::new();
228        let mut ln_current = 0.;
229        let mut current_value = (1 << 62) / 9;
230        while ln_current < ln_pseudo_norm_bound + 1. {
231            current_value = prev_prime(StaticRing::<i64>::RING, current_value).unwrap();
232            if current_value < (1 << 32) {
233                panic!("not enough primes");
234            }
235            primes.push(current_value);
236            ln_current += (current_value as f64).ln();
237        }
238        return zn_rns::Zn::new(
239            primes.into_iter().map(|p| AsField::from(AsFieldBase::promise_is_perfect_field(zn_64::Zn::new(p as u64)))).collect(),
240            RingRef::new(self)
241        );
242    }
243
244    fn local_ring_at<'ring>(&self, computation: &Self::LocalComputationData<'ring>, i: usize) -> Self::LocalRing<'ring>
245        where Self: 'ring
246    {
247        computation.at(i).clone()
248    }
249
250    fn local_ring_count<'ring>(&self, computation: &Self::LocalComputationData<'ring>) -> usize
251        where Self: 'ring
252    {
253        computation.len()
254    }
255
256    fn reduce<'ring>(&self, computation: &Self::LocalComputationData<'ring>, el: &Self::Element) -> Vec<<Self::LocalRingBase<'ring> as RingBase>::Element>
257        where Self: 'ring
258    {
259        computation.get_congruence(&computation.coerce(RingValue::from_ref(self), self.clone_el(el))).as_iter().map(|x| *x).collect()
260    }
261
262    fn lift_combine<'ring>(&self, computation: &Self::LocalComputationData<'ring>, el: &[<Self::LocalRingBase<'ring> as RingBase>::Element]) -> Self::Element
263        where Self: 'ring
264    {
265        computation.smallest_lift(computation.from_congruence(el.iter().copied()))
266    }
267}
268
269///
270/// The map `R -> R/p` for a ring `R` and one of its local quotients `R/p` as
271/// given by [`EvaluatePolyLocallyRing`].
272/// 
273#[stability::unstable(feature = "enable")]
274pub struct EvaluatePolyLocallyReductionMap<'ring, 'data, R>
275    where R: 'ring + ?Sized + EvaluatePolyLocallyRing, 'ring: 'data
276{
277    ring: RingRef<'data, R>,
278    data: &'data R::LocalComputationData<'ring>,
279    local_ring: R::LocalRing<'ring>,
280    index: usize
281}
282
283impl<'ring, 'data, R> EvaluatePolyLocallyReductionMap<'ring, 'data, R>
284    where R: 'ring +?Sized + EvaluatePolyLocallyRing, 'ring: 'data
285{
286    #[stability::unstable(feature = "enable")]
287    pub fn new(ring: &'data R, data: &'data R::LocalComputationData<'ring>, index: usize) -> Self {
288        Self { ring: RingRef::new(ring), data: data, local_ring: ring.local_ring_at(data, index), index: index }
289    }
290}
291
292impl<'ring, 'data, R> Homomorphism<R, R::LocalRingBase<'ring>> for EvaluatePolyLocallyReductionMap<'ring, 'data, R>
293    where R: 'ring +?Sized + EvaluatePolyLocallyRing, 'ring: 'data
294{
295    type CodomainStore = R::LocalRing<'ring>;
296    type DomainStore = RingRef<'data, R>;
297
298    fn codomain<'b>(&'b self) -> &'b Self::CodomainStore {
299        &self.local_ring
300    }
301
302    fn domain<'b>(&'b self) -> &'b Self::DomainStore {
303        &self.ring
304    }
305
306    fn map(&self, x: <R as RingBase>::Element) -> <R::LocalRingBase<'ring> as RingBase>::Element {
307        let ring_ref: &'data R = self.ring.into();
308        let mut reductions: Vec<<R::LocalRingBase<'ring> as RingBase>::Element> = ring_ref.reduce(self.data, &x);
309        return reductions.swap_remove(self.index);
310    }
311}
312
313#[macro_export]
314macro_rules! impl_interpolation_base_ring_char_zero {
315    (<{$($gen_args:tt)*}> InterpolationBaseRing for $self_type:ty where $($constraints:tt)*) => {
316        impl<$($gen_args)*> $crate::compute_locally::InterpolationBaseRing for $self_type where $($constraints)* {
317                
318            type ExtendedRing<'a> = RingRef<'a, Self>
319                where Self: 'a;
320
321            type ExtendedRingBase<'a> = Self
322                where Self: 'a;
323
324            fn in_base<'a, S>(&self, _ext_ring: S, el: El<S>) -> Option<Self::Element>
325                where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>
326            {
327                Some(el)
328            }
329
330            fn in_extension<'a, S>(&self, _ext_ring: S, el: Self::Element) -> El<S>
331                where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>
332            {
333                el
334            }
335
336            fn interpolation_points<'a>(&'a self, count: usize) -> (Self::ExtendedRing<'a>, Vec<El<Self::ExtendedRing<'a>>>) {
337                let ZZbig = $crate::integer::BigIntRing::RING;
338                assert!(ZZbig.is_zero(&self.characteristic(&ZZbig).unwrap()));
339                let ring = $crate::ring::RingRef::new(self);
340                (ring, (0..count).map(|n| <_ as $crate::homomorphism::Homomorphism<_, _>>::map(&ring.int_hom(), n as i32)).collect())
341            }
342        }
343    };
344}
feanor_math/compute_locally.rs

feanor_math/
compute_locally.rs
