feanor_math/
integer.rs

1use std::fmt::Debug;
2
3use crate::reduce_lift::poly_factor_gcd::IntegerPolyGCDRing;
4use crate::reduce_lift::poly_eval::EvalPolyLocallyRing;
5use crate::divisibility::*;
6use crate::ring::*;
7use crate::homomorphism::*;
8use crate::pid::*;
9use crate::ordered::*;
10
11///
12/// Type alias for the current default used big integer ring implementation.
13/// 
14/// The type this points to may change when features or other compilation parameters
15/// change.
16///
17#[cfg(feature = "mpir")]
18pub type BigIntRing = crate::rings::mpir::MPZ;
19///
20/// Type alias for the current default used big integer ring implementation.
21/// 
22/// The type this points to may change when features or other compilation parameters
23/// change.
24/// 
25#[cfg(not(feature = "mpir"))]
26pub type BigIntRing = crate::rings::rust_bigint::RustBigintRing;
27///
28/// Type alias for the current default used big integer ring implementation.
29/// 
30/// The type this points to may change when features or other compilation parameters
31/// change.
32/// 
33#[cfg(feature = "mpir")]
34pub type BigIntRingBase = crate::rings::mpir::MPZBase;
35///
36/// Type alias for the current default used big integer ring implementation.
37/// 
38/// The type this points to may change when features or other compilation parameters
39/// change.
40/// 
41#[cfg(not(feature = "mpir"))]
42pub type BigIntRingBase = crate::rings::rust_bigint::RustBigintRingBase;
43
44///
45/// Trait for rings that are isomorphic to the ring of integers `ZZ = { ..., -2, -1, 0, 1, 2, ... }`.
46/// 
47/// Some of the functionality in this trait refers to the binary expansion of
48/// a positive integer. While this is not really general, it is often required
49/// for fast operations with integers.
50/// 
51/// As an additional requirement, the euclidean division (i.e. [`EuclideanRing::euclidean_div_rem()`] and
52/// [`IntegerRing::euclidean_div_pow_2()`]) are additionally expected to round towards zero.
53/// 
54/// Currently [`IntegerRing`] is a subtrait of the unstable traits [`IntegerPolyGCDRing`] and,
55/// [`EvalPolyLocallyRing`] so it is at the moment impossible to implement [`IntegerRing`] for a
56/// custom ring type without enabling unstable features. Sorry.
57/// 
58pub trait IntegerRing: Domain + EuclideanRing + OrderedRing + HashableElRing + IntegerPolyGCDRing + EvalPolyLocallyRing + Debug {
59
60    ///
61    /// Computes a float value that is "close" to the given integer.
62    /// 
63    /// However, no guarantees are made on how close it must be, in particular,
64    /// this function may also always return `0.` (this is just an example - 
65    /// it's not a good idea).
66    /// 
67    /// Some use cases include:
68    ///  - Estimating control parameters (e.g. bounds for prime numbers
69    ///    during factoring algorithms)
70    ///  - First performing some operation on floating point numbers, and
71    ///    then refining it to integers.
72    /// 
73    /// Note that a high-quality implementation of this function can vastly improve
74    /// performance in some cases, e.g. of [`crate::algorithms::int_bisect::root_floor()`] or 
75    /// [`crate::algorithms::lll::exact::lll()`].
76    /// 
77    /// # Example
78    /// 
79    /// ```rust
80    /// # use feanor_math::ring::*;
81    /// # use feanor_math::integer::*;
82    /// let ZZ = BigIntRing::RING;
83    /// let x = ZZ.power_of_two(1023);
84    /// assert!(ZZ.to_float_approx(&x) > 2f64.powi(1023) * 0.99999);
85    /// assert!(ZZ.to_float_approx(&x) < 2f64.powi(1023) * 1.000001);
86    /// ```
87    /// If the value is too large for the exponent of a `f64`, infinity is returned.
88    /// ```rust
89    /// # use feanor_math::ring::*;
90    /// # use feanor_math::integer::*;
91    /// let ZZ = BigIntRing::RING;
92    /// let x = ZZ.power_of_two(1024);
93    /// assert!(ZZ.to_float_approx(&x).is_infinite());
94    /// ```
95    /// 
96    fn to_float_approx(&self, value: &Self::Element) -> f64;
97
98    ///
99    /// Computes a value that is "close" to the given float. However, no guarantees
100    /// are made on the definition of close, in particular, this does not have to be
101    /// the closest integer to the given float, and cannot be used to compute rounding.
102    /// It is also implementation-defined when to return `None`, although this is usually
103    /// the case on infinity and NaN.
104    /// 
105    /// For information when to use (or not use) this, see its counterpart [`IntegerRing::to_float_approx()`].
106    /// 
107    fn from_float_approx(&self, value: f64) -> Option<Self::Element>;
108
109    ///
110    /// Return whether the `i`-th bit in the two-complements representation of `abs(value)`
111    /// is `1`.
112    /// 
113    /// # Example
114    /// ```rust
115    /// # use feanor_math::primitive_int::*;
116    /// # use feanor_math::integer::*;
117    /// # use feanor_math::ring::*;
118    /// assert_eq!(false, StaticRing::<i32>::RING.abs_is_bit_set(&4, 1));
119    /// assert_eq!(true, StaticRing::<i32>::RING.abs_is_bit_set(&4, 2));
120    /// assert_eq!(true, StaticRing::<i32>::RING.abs_is_bit_set(&-4, 2));
121    /// ```
122    /// 
123    fn abs_is_bit_set(&self, value: &Self::Element, i: usize) -> bool;
124
125    ///
126    /// Returns the index of the highest set bit in the two-complements representation of `abs(value)`,
127    /// or `None` if the value is zero.
128    /// 
129    /// # Example
130    /// ```rust
131    /// # use feanor_math::primitive_int::*;
132    /// # use feanor_math::integer::*;
133    /// # use feanor_math::ring::*;
134    /// assert_eq!(None, StaticRing::<i32>::RING.abs_highest_set_bit(&0));
135    /// assert_eq!(Some(0), StaticRing::<i32>::RING.abs_highest_set_bit(&-1));
136    /// assert_eq!(Some(2), StaticRing::<i32>::RING.abs_highest_set_bit(&4));
137    /// ```
138    /// 
139    fn abs_highest_set_bit(&self, value: &Self::Element) -> Option<usize>;
140
141    ///
142    /// Returns the index of the lowest set bit in the two-complements representation of `abs(value)`,
143    /// or `None` if the value is zero.
144    /// 
145    /// # Example
146    /// ```rust
147    /// # use feanor_math::primitive_int::*;
148    /// # use feanor_math::integer::*;
149    /// # use feanor_math::ring::*;
150    /// assert_eq!(None, StaticRing::<i32>::RING.abs_lowest_set_bit(&0));
151    /// assert_eq!(Some(0), StaticRing::<i32>::RING.abs_lowest_set_bit(&1));
152    /// assert_eq!(Some(0), StaticRing::<i32>::RING.abs_lowest_set_bit(&-3));
153    /// assert_eq!(Some(2), StaticRing::<i32>::RING.abs_lowest_set_bit(&4));
154    /// ```
155    /// 
156    fn abs_lowest_set_bit(&self, value: &Self::Element) -> Option<usize>;
157
158    ///
159    /// Computes the euclidean division by a power of two, always rounding to zero (note that
160    /// euclidean division requires that `|remainder| < |divisor|`, and thus would otherwise
161    /// leave multiple possible results).
162    /// 
163    /// # Example
164    /// ```rust
165    /// # use feanor_math::primitive_int::*;
166    /// # use feanor_math::integer::*;
167    /// # use feanor_math::ring::*;
168    /// let mut value = -7;
169    /// StaticRing::<i32>::RING.euclidean_div_pow_2(&mut value, 1);
170    /// assert_eq!(-3, value);
171    /// ```
172    /// 
173    fn euclidean_div_pow_2(&self, value: &mut Self::Element, power: usize);
174
175    ///
176    /// Multiplies the element by a power of two.
177    /// 
178    fn mul_pow_2(&self, value: &mut Self::Element, power: usize);
179
180    ///
181    /// Computes a uniformly random integer in `[0, 2^log_bound_exclusive - 1]`, assuming that
182    /// `rng` provides uniformly random values in the whole range of `u64`.
183    /// 
184    fn get_uniformly_random_bits<G: FnMut() -> u64>(&self, log2_bound_exclusive: usize, rng: G) -> Self::Element;
185
186    ///
187    /// Computes the rounded division, i.e. rounding to the closest integer.
188    /// In the case of a tie (i.e. `round(0.5)`), we round towards `+/- infinity`.
189    /// 
190    /// # Example
191    /// ```rust
192    /// # use feanor_math::primitive_int::*;
193    /// # use feanor_math::integer::*;
194    /// # use feanor_math::ring::*;
195    /// assert_eq!(2, StaticRing::<i32>::RING.rounded_div(7, &3));
196    /// assert_eq!(-2, StaticRing::<i32>::RING.rounded_div(-7, &3));
197    /// assert_eq!(-2, StaticRing::<i32>::RING.rounded_div(7, &-3));
198    /// assert_eq!(2, StaticRing::<i32>::RING.rounded_div(-7, &-3));
199    /// 
200    /// assert_eq!(3, StaticRing::<i32>::RING.rounded_div(8, &3));
201    /// assert_eq!(-3, StaticRing::<i32>::RING.rounded_div(-8, &3));
202    /// assert_eq!(-3, StaticRing::<i32>::RING.rounded_div(8, &-3));
203    /// assert_eq!(3, StaticRing::<i32>::RING.rounded_div(-8, &-3));
204    /// 
205    /// assert_eq!(4, StaticRing::<i32>::RING.rounded_div(7, &2));
206    /// assert_eq!(-4, StaticRing::<i32>::RING.rounded_div(-7, &2));
207    /// assert_eq!(-4, StaticRing::<i32>::RING.rounded_div(7, &-2));
208    /// assert_eq!(4, StaticRing::<i32>::RING.rounded_div(-7, &-2));
209    /// ```
210    /// 
211    fn rounded_div(&self, lhs: Self::Element, rhs: &Self::Element) -> Self::Element {
212        let mut rhs_half = self.abs(self.clone_el(rhs));
213        self.euclidean_div_pow_2(&mut rhs_half, 1);
214        if self.is_neg(&lhs) {
215            return self.euclidean_div(self.sub(lhs, rhs_half), rhs);
216        } else {
217            return self.euclidean_div(self.add(lhs, rhs_half), rhs);
218        }
219    }
220
221    ///
222    /// Computes the division `lhs / rhs`, rounding towards `+ infinity`.
223    /// 
224    /// In particular, if `rhs` is positive, this gives the smallest
225    /// integer `quo` such that `quo * rhs >= lhs`. On the other hand, if
226    /// `rhs` is negative, this computes the largest integer `quo` such that
227    /// `quo * rhs <= lhs`.
228    /// 
229    /// # Example
230    /// ```rust
231    /// # use feanor_math::primitive_int::*;
232    /// # use feanor_math::integer::*;
233    /// # use feanor_math::ring::*;
234    /// assert_eq!(3, StaticRing::<i32>::RING.ceil_div(7, &3));
235    /// assert_eq!(-2, StaticRing::<i32>::RING.ceil_div(-7, &3));
236    /// assert_eq!(-2, StaticRing::<i32>::RING.ceil_div(7, &-3));
237    /// assert_eq!(3, StaticRing::<i32>::RING.ceil_div(-7, &-3));
238    /// ```
239    /// 
240    fn ceil_div(&self, lhs: Self::Element, rhs: &Self::Element) -> Self::Element {
241        assert!(!self.is_zero(rhs));
242        if self.is_zero(&lhs) {
243            return self.zero();
244        }
245        let one = self.one();
246        return match (self.is_pos(&lhs), self.is_pos(rhs)) {
247            (true, true) => self.add(self.euclidean_div(self.sub_ref_snd(lhs, &one), rhs), one),
248            (false, true) => self.euclidean_div(lhs, rhs),
249            (true, false) => self.euclidean_div(lhs, rhs),
250            (false, false) => self.add(self.euclidean_div(self.add_ref_snd(lhs, &one), rhs), one)
251        };
252    }
253
254    ///
255    /// Computes the division `lhs / rhs`, rounding towards `- infinity`.
256    /// 
257    /// In particular, if `rhs` is positive, this gives the largest
258    /// integer `quo` such that `quo * rhs <= lhs`. On the other hand, if
259    /// `rhs` is negative, this computes the smallest integer `quo` such that
260    /// `quo * rhs >= lhs`.
261    /// 
262    /// # Example
263    /// ```rust
264    /// # use feanor_math::primitive_int::*;
265    /// # use feanor_math::integer::*;
266    /// # use feanor_math::ring::*;
267    /// assert_eq!(2, StaticRing::<i32>::RING.floor_div(7, &3));
268    /// assert_eq!(-3, StaticRing::<i32>::RING.floor_div(-7, &3));
269    /// assert_eq!(-3, StaticRing::<i32>::RING.floor_div(7, &-3));
270    /// assert_eq!(2, StaticRing::<i32>::RING.floor_div(-7, &-3));
271    /// ```
272    /// 
273    fn floor_div(&self, lhs: Self::Element, rhs: &Self::Element) -> Self::Element {
274        self.negate(self.ceil_div(self.negate(lhs), rhs))
275    }
276
277    ///
278    /// Returns the value `2^power` in this integer ring.
279    /// 
280    fn power_of_two(&self, power: usize) -> Self::Element {
281        let mut result = self.one();
282        self.mul_pow_2(&mut result, power);
283        return result;
284    }
285
286    ///
287    /// Returns `n` such that this ring can represent at least `[-2^n, ..., 2^n - 1]`.
288    /// Returning `None` means that the size of representable integers is unbounded.
289    /// 
290    fn representable_bits(&self) -> Option<usize>;
291
292    ///
293    /// Parses the given string as a number.
294    /// 
295    /// Returns `Err(())` if it is not a valid number w.r.t. base, i.e. if the string
296    /// is not a sequence of digit characters, optionally beginning with `+` or `-`. 
297    /// To denote digits larger than `9`, the same characters as in [`u64::from_str_radix()`]
298    /// should be used.
299    /// 
300    fn parse(&self, string: &str, base: u32) -> Result<Self::Element, ()> {
301        generic_impls::parse(RingRef::new(self), string, base)
302    }
303}
304
305impl<I, J> CanHomFrom<I> for J
306    where I: ?Sized + IntegerRing, J: ?Sized + IntegerRing
307{
308    type Homomorphism = ();
309
310    fn has_canonical_hom(&self, _: &I) -> Option<Self::Homomorphism> {
311        Some(())
312    }
313
314    fn map_in(&self, from: &I, el: <I as RingBase>::Element, _: &Self::Homomorphism) -> Self::Element {
315        int_cast(el, &RingRef::new(self), &RingRef::new(from))
316    }
317
318    default fn map_in_ref(&self, from: &I, el: &<I as RingBase>::Element, hom: &Self::Homomorphism) -> Self::Element {
319        <J as CanHomFrom<I>>::map_in(self, from, from.clone_el(el), hom)
320    }
321}
322
323impl<I, J> CanIsoFromTo<I> for J
324    where I: ?Sized + IntegerRing, J: ?Sized + IntegerRing
325{
326    type Isomorphism = ();
327
328    fn has_canonical_iso(&self, _: &I) -> Option<Self::Isomorphism> {
329        Some(())
330    }
331
332    fn map_out(&self, from: &I, el: Self::Element, _: &Self::Isomorphism) -> <I as RingBase>::Element {
333        int_cast(el, &RingRef::new(from), &RingRef::new(self))
334    }
335}
336
337///
338/// Helper trait to simplify conversion between ints.
339/// 
340/// More concretely, `IntCast` defines a conversion between two
341/// integer rings, and is default-implemented for all integer rings
342/// using a double-and-and technique. All implementors of integer
343/// rings are encouraged to provide specializations for improved performance.
344/// 
345/// # Why yet another conversion trait?
346/// 
347/// It is a common requirement to convert between arbitrary (i.e. generic)
348/// integer rings. To achieve this, we require a blanket implementation
349/// anyway.
350/// 
351/// Now it would be possible to just provide a blanket implementation of
352/// [`CanHomFrom`] and specialize it for all integer rings. However, specialization
353/// with default types is currently a pain in the ass. Furthermore, this trait is simpler.
354/// 
355pub trait IntCast<F: ?Sized + IntegerRing>: IntegerRing {
356
357    ///
358    /// Maps the given integer into this ring.
359    /// 
360    /// For the difference to [`RingStore::coerce()`] or [`RingStore::can_hom()`],
361    /// see the documentation of [`IntCast`].
362    /// 
363    fn cast(&self, from: &F, value: F::Element) -> Self::Element;
364}
365
366impl<F: ?Sized + IntegerRing, T: ?Sized + IntegerRing> IntCast<F> for T {
367
368    default fn cast(&self, from: &F, value: F::Element) -> Self::Element {
369        generic_impls::map_from_integer_ring(RingRef::new(from), RingRef::new(self), value)
370    }
371}
372
373///
374/// Conversion of elements between two rings representing the integers `ZZ`.
375/// 
376/// The underlying conversion functionality is the same as provided by [`IntCast`], and
377/// indirectly also by [`CanHomFrom`] and [`CanIsoFromTo`].
378/// 
379/// # Example
380/// 
381/// ```rust
382/// # use feanor_math::ring::*;
383/// # use feanor_math::integer::*;
384/// # use feanor_math::primitive_int::*;
385/// # use feanor_math::assert_el_eq;
386/// let ZZi64 = StaticRing::<i64>::RING;
387/// let ZZbig = BigIntRing::RING;
388/// let ZZi8 = StaticRing::<i8>::RING;
389/// assert_eq!(7, int_cast(7, ZZi64, ZZi8));
390/// assert_eq!(65536, int_cast(ZZbig.power_of_two(16), ZZi64, ZZbig));
391/// assert_el_eq!(ZZbig, ZZbig.power_of_two(16), int_cast(65536, ZZbig, ZZi64));
392///  ```
393/// 
394pub fn int_cast<T: RingStore, F: RingStore>(value: El<F>, to: T, from: F) -> El<T>
395    where T::Type: IntegerRing, F::Type: IntegerRing
396{
397    <T::Type as IntCast<F::Type>>::cast(to.get_ring(), from.get_ring(), value)
398}
399
400///
401/// Computes the binomial coefficient of `n` and `k`, defined as `n(n - 1)...(n - k + 1)/k!`.
402/// 
403/// The above definition works for any `n` and `k >= 0`. If `k < 0`, we define the binomial coefficient
404/// to be zero. This function will not overflow, if the integer rings supports number up to 
405/// `binomial(n, k) * k`.
406/// 
407/// # Example
408/// ```rust
409/// # use feanor_math::ring::*;
410/// # use feanor_math::integer::*;
411/// # use feanor_math::iters::*;
412/// # use feanor_math::primitive_int::*;
413/// // the binomial coefficient is equal to the number of combinations of fixed size
414/// assert_eq!(
415///     binomial(10, &3, StaticRing::<i64>::RING) as usize,
416///     multiset_combinations(&[1; 10], 3, |_| ()).count()
417/// );
418/// ```
419/// 
420#[stability::unstable(feature = "enable")]
421pub fn binomial<I>(n: El<I>, k: &El<I>, ring: I) -> El<I>
422    where I: RingStore + Copy,
423        I::Type: IntegerRing
424{
425    if ring.is_neg(&n) {
426        let mut result = binomial(ring.sub(ring.sub_ref_fst(&k, n), ring.one()), k, ring);
427        if !ring.is_even(k) {
428            ring.negate_inplace(&mut result);
429        }
430        return result;
431    } else if ring.is_neg(k) || ring.is_gt(k, &n) {
432        return ring.zero();
433    } else {
434        // this formula works always, and is guaranteed not to overflow if k <= n/2 and `binomial(n, k) * k` 
435        // fits into an integer; thus distinguish this case that k > n/2
436        let n_neg_k = ring.sub_ref(&n, &k);
437        if ring.is_lt(&n_neg_k, k) {
438            return binomial(n, &n_neg_k, ring);
439        }
440        let mut result = ring.one();
441        let mut i = ring.one();
442        while ring.is_leq(&i, &k) {
443            ring.mul_assign(&mut result, ring.sub_ref_snd(ring.add_ref_fst(&n, ring.one()), &i));
444            result = ring.checked_div(&result, &i).unwrap();
445            ring.add_assign(&mut i, ring.one());
446        }
447        return result;
448    }
449}
450
451#[stability::unstable(feature = "enable")]
452pub fn factorial<I>(n: &El<I>, ring: I) -> El<I>
453    where I: RingStore + Copy,
454        I::Type: IntegerRing
455{
456    let mut current = ring.zero();
457    let one = ring.one();
458    return ring.prod((0..).map_while(|_| {
459        if ring.is_lt(&current, &n) {
460            ring.add_assign_ref(&mut current, &one);
461            return Some(ring.clone_el(&current));
462        } else {
463            return None;
464        }
465    }));
466}
467
468///
469/// Trait for [`RingStore`]s that store [`IntegerRing`]s. Mainly used
470/// to provide a convenient interface to the `IntegerRing`-functions.
471/// 
472pub trait IntegerRingStore: RingStore
473    where Self::Type: IntegerRing
474{
475    delegate!{ IntegerRing, fn to_float_approx(&self, value: &El<Self>) -> f64 }
476    delegate!{ IntegerRing, fn from_float_approx(&self, value: f64) -> Option<El<Self>> }
477    delegate!{ IntegerRing, fn abs_is_bit_set(&self, value: &El<Self>, i: usize) -> bool }
478    delegate!{ IntegerRing, fn abs_highest_set_bit(&self, value: &El<Self>) -> Option<usize> }
479    delegate!{ IntegerRing, fn abs_lowest_set_bit(&self, value: &El<Self>) -> Option<usize> }
480    delegate!{ IntegerRing, fn euclidean_div_pow_2(&self, value: &mut El<Self>, power: usize) -> () }
481    delegate!{ IntegerRing, fn mul_pow_2(&self, value: &mut El<Self>, power: usize) -> () }
482    delegate!{ IntegerRing, fn power_of_two(&self, power: usize) -> El<Self> }
483    delegate!{ IntegerRing, fn rounded_div(&self, lhs: El<Self>, rhs: &El<Self>) -> El<Self> }
484    delegate!{ IntegerRing, fn floor_div(&self, lhs: El<Self>, rhs: &El<Self>) -> El<Self> }
485    delegate!{ IntegerRing, fn ceil_div(&self, lhs: El<Self>, rhs: &El<Self>) -> El<Self> }
486    delegate!{ IntegerRing, fn parse(&self, string: &str, base: u32) -> Result<El<Self>, ()> }
487
488    ///
489    /// Using the randomness of the given rng, samples a uniformly random integer
490    /// from the set `{ 0, 1, ..., bound_exclusive - 1 }`.
491    /// 
492    /// Uses rejection sampling on top of [`IntegerRing::get_uniformly_random_bits()`].
493    /// 
494    fn get_uniformly_random<G: FnMut() -> u64>(&self, bound_exclusive: &El<Self>, mut rng: G) -> El<Self> {
495        assert!(self.is_gt(bound_exclusive, &self.zero()));
496        let log2_ceil_bound = self.abs_highest_set_bit(bound_exclusive).unwrap() + 1;
497        let mut result = self.get_ring().get_uniformly_random_bits(log2_ceil_bound, || rng());
498        while self.is_geq(&result, bound_exclusive) {
499            result = self.get_ring().get_uniformly_random_bits(log2_ceil_bound, || rng());
500        }
501        return result;
502    }
503
504    ///
505    /// Computes `floor(log2(abs(value)))`, and returns `None` if the argument is 0.
506    /// 
507    /// This is equivalent to [`IntegerRingStore::abs_highest_set_bit`].
508    /// 
509    fn abs_log2_floor(&self, value: &El<Self>) -> Option<usize> {
510        self.abs_highest_set_bit(value)
511    }
512
513    ///
514    /// Computes `ceil(log2(abs(value)))`, and returns `None` if the argument is 0.
515    /// 
516    fn abs_log2_ceil(&self, value: &El<Self>) -> Option<usize> {
517        let highest_bit = self.abs_highest_set_bit(value)?;
518        if self.abs_lowest_set_bit(value).unwrap() == highest_bit {
519            return Some(highest_bit);
520        } else {
521            return Some(highest_bit + 1);
522        }
523    }
524
525    ///
526    /// Returns true if the given integer is divisible by 2.
527    /// 
528    fn is_even(&self, value: &El<Self>) -> bool {
529        !self.abs_is_bit_set(value, 0)
530    }
531
532    ///
533    /// Returns true if the given integer is not divisible by 2.
534    /// 
535    fn is_odd(&self, value: &El<Self>) -> bool {
536        !self.is_even(value)
537    }
538
539    ///
540    /// Assumes the given integer is even, and computes its quotient by 2.
541    /// 
542    fn half_exact(&self, mut value: El<Self>) -> El<Self> {
543        assert!(self.is_even(&value));
544        self.euclidean_div_pow_2(&mut value, 1);
545        return value;
546    }
547}
548
549impl<R> IntegerRingStore for R
550    where R: RingStore,
551        R::Type: IntegerRing
552{}
553
554///
555/// Implementations of [`IntegerRing`] functionality that make sense for
556/// any [`IntegerRing`], and thus can be used as default if no better implementation
557/// can be provided for a concrete [`IntegerRing`].
558/// 
559pub mod generic_impls {
560    use crate::ring::*;
561    use crate::primitive_int::*;
562    use super::*;
563    
564    #[stability::unstable(feature = "enable")]
565    pub fn map_from_integer_ring<I, R>(from: I, to: R, mut x: El<I>) -> El<R>
566        where I: RingStore,
567            I::Type: IntegerRing,
568            R: RingStore
569    {
570        let basis = to.int_hom().map(1 << 16);
571        let is_neg = if from.is_neg(&x) {
572            from.negate_inplace(&mut x);
573            true
574        } else {
575            false
576        };
577        let mut current = to.zero();
578        let mut current_pow = to.one();
579        while !from.is_zero(&x) {
580            let mut quo = from.clone_el(&x);
581            from.euclidean_div_pow_2(&mut quo, 16);
582            let mut rem = from.clone_el(&quo);
583            from.mul_pow_2(&mut rem, 16);
584            from.sub_self_assign(&mut rem, x);
585            let rem = int_cast(rem, StaticRing::<i32>::RING, &from);
586            to.add_assign(&mut current, to.mul_ref_snd(to.int_hom().map(rem), &current_pow));
587            x = quo;
588            to.mul_assign_ref(&mut current_pow, &basis);
589        }
590        if is_neg {
591            return to.negate(current);
592        } else {
593            return current;
594        }
595    }
596
597    #[stability::unstable(feature = "enable")]
598    pub fn parse<I>(ring: I, string: &str, base: u32) -> Result<El<I>, ()>
599        where I: RingStore, I::Type: IntegerRing
600    {
601        let (negative, rest) = if string.chars().next() == Some('-') {
602            (true, string.split_at(1).1)
603        } else if string.chars().next() == Some('+') {
604            (false, string.split_at(1).1)
605        } else {
606            (false, string)
607        };
608        assert!(base >= 2);
609
610        let bits_per_chunk = u32::BITS as usize;
611        assert!(bits_per_chunk <= i64::BITS as usize - 2);
612        let chunk_size = (bits_per_chunk as f32 / (base as f32).log2()).floor() as usize;
613        let it = <str as AsRef<[u8]>>::as_ref(rest)
614            .rchunks(chunk_size)
615            .rev()
616            .map(std::str::from_utf8)
617            .map(|chunk| chunk.map_err(|_| ()))
618            .map(|chunk| chunk.and_then(|n| 
619                u64::from_str_radix(n, base).map_err(|_| ()))
620            );
621        let chunk_base = ring.pow(int_cast(base as i64, &ring, StaticRing::<i64>::RING), chunk_size);
622        let result = it.fold(Ok(ring.zero()), |current, next| {
623            current.and_then(|mut current| next.map(|next| {
624                ring.mul_assign_ref(&mut current, &chunk_base);
625                ring.add(current, int_cast(next as i64, &ring, StaticRing::<i64>::RING))
626            }))
627        });
628        if negative {
629            return result.map(|result| ring.negate(result));
630        } else {
631            return result;
632        }
633    }
634}
635
636#[cfg(test)]
637use crate::primitive_int::*;
638#[cfg(test)]
639use generic_impls::map_from_integer_ring;
640
641#[allow(missing_docs)]
642#[cfg(any(test, feature = "generic_tests"))]
643pub mod generic_tests {
644
645    use crate::ring::El;
646    use super::*;
647        
648    pub fn test_integer_get_uniformly_random<R: RingStore>(ring: R) 
649        where R::Type: IntegerRing
650    {
651        for b in [15, 16] {
652            let bound = ring.int_hom().map(b);
653            let mut rng = oorandom::Rand64::new(1);
654            let elements: Vec<El<R>> = (0..1000).map(|_| ring.get_uniformly_random(&bound, || rng.rand_u64())).collect();
655            for i in 0..b {
656                assert!(elements.iter().any(|x| ring.eq_el(x, &ring.int_hom().map(i))))
657            }
658            for x in &elements {
659                assert!(ring.is_lt(x, &bound));
660            }
661        }
662    }
663
664    pub fn test_integer_axioms<R: IntegerRingStore, I: Iterator<Item = El<R>>>(ring: R, edge_case_elements: I) 
665        where R::Type: IntegerRing
666    {
667        let elements = edge_case_elements.collect::<Vec<_>>();
668
669        // test abs_highest_set_bit on standard values
670        assert_eq!(None, ring.abs_highest_set_bit(&ring.int_hom().map(0)));
671        assert_eq!(Some(0), ring.abs_highest_set_bit(&ring.int_hom().map(1)));
672        assert_eq!(Some(1), ring.abs_highest_set_bit(&ring.int_hom().map(2)));
673
674        // generic test of mul_pow_2 resp. euclidean_div_pow_2
675        for a in &elements {
676            let mut ceil_pow_2 = ring.int_hom().map(2);
677            ring.mul_pow_2(&mut ceil_pow_2, ring.abs_highest_set_bit(a).unwrap_or(0));
678            assert!(ring.is_lt(a, &ceil_pow_2));
679            assert!(ring.is_lt(&ring.negate(ring.clone_el(a)), &ceil_pow_2));
680            
681            for i in 0..ring.abs_highest_set_bit(a).unwrap_or(0) {
682                let mut pow_2 = ring.one();
683                ring.mul_pow_2(&mut pow_2, i);
684                let mut b = ring.clone_el(a);
685                ring.mul_pow_2(&mut b, i);
686                assert_el_eq!(ring, b, ring.mul(ring.clone_el(a), ring.clone_el(&pow_2)));
687                ring.euclidean_div_pow_2(&mut b, i);
688                assert_el_eq!(ring, b, a);
689                ring.euclidean_div_pow_2(&mut b, i);
690                assert_el_eq!(ring, b, ring.euclidean_div(ring.clone_el(a), &pow_2));
691            }
692        }
693
694        // test euclidean div round to zero
695        let d = ring.int_hom().map(8);
696        for k in -10..=10 {
697            let mut a = ring.int_hom().map(k);
698            assert_el_eq!(ring, ring.int_hom().map(k / 8), ring.euclidean_div(ring.clone_el(&a), &d));
699            ring.euclidean_div_pow_2(&mut a, 3);
700            assert_el_eq!(ring, ring.int_hom().map(k / 8), a);
701        }
702        let d = ring.int_hom().map(-8);
703        for k in -10..=10 {
704            let a = ring.int_hom().map(k);
705            assert_el_eq!(ring, ring.int_hom().map(k / -8), ring.euclidean_div(ring.clone_el(&a), &d));
706        }
707
708        // test rounded_div
709        assert_el_eq!(ring, ring.int_hom().map(2), ring.rounded_div(ring.int_hom().map(7), &ring.int_hom().map(3)));
710        assert_el_eq!(ring, ring.int_hom().map(-2), ring.rounded_div(ring.int_hom().map(-7), &ring.int_hom().map(3)));
711        assert_el_eq!(ring, ring.int_hom().map(-2), ring.rounded_div(ring.int_hom().map(7), &ring.int_hom().map(-3)));
712        assert_el_eq!(ring, ring.int_hom().map(2), ring.rounded_div(ring.int_hom().map(-7), &ring.int_hom().map(-3)));
713
714        assert_el_eq!(ring, ring.int_hom().map(3), ring.rounded_div(ring.int_hom().map(8), &ring.int_hom().map(3)));
715        assert_el_eq!(ring, ring.int_hom().map(-3), ring.rounded_div(ring.int_hom().map(-8), &ring.int_hom().map(3)));
716        assert_el_eq!(ring, ring.int_hom().map(-3), ring.rounded_div(ring.int_hom().map(8), &ring.int_hom().map(-3)));
717        assert_el_eq!(ring, ring.int_hom().map(3), ring.rounded_div(ring.int_hom().map(-8), &ring.int_hom().map(-3)));
718
719        assert_el_eq!(ring, ring.int_hom().map(4), ring.rounded_div(ring.int_hom().map(7), &ring.int_hom().map(2)));
720        assert_el_eq!(ring, ring.int_hom().map(-4), ring.rounded_div(ring.int_hom().map(-7), &ring.int_hom().map(2)));
721        assert_el_eq!(ring, ring.int_hom().map(-4), ring.rounded_div(ring.int_hom().map(7), &ring.int_hom().map(-2)));
722        assert_el_eq!(ring, ring.int_hom().map(4), ring.rounded_div(ring.int_hom().map(-7), &ring.int_hom().map(-2)));
723    }
724}
725
726#[test]
727fn test_int_div_assumption() {
728    assert_eq!(-1, -10 / 8);
729    assert_eq!(-1, 10 / -8);
730    assert_eq!(1, 10 / 8);
731    assert_eq!(1, -10 / -8);
732}
733
734#[test]
735fn test_rounded_div() {
736    let ZZ = StaticRing::<i32>::RING;
737    assert_el_eq!(ZZ, 3, ZZ.rounded_div(20, &7));
738    assert_el_eq!(ZZ, -3, ZZ.rounded_div(-20, &7));
739    assert_el_eq!(ZZ, -3, ZZ.rounded_div(20, &-7));
740    assert_el_eq!(ZZ, 3, ZZ.rounded_div(-20, &-7));
741    assert_el_eq!(ZZ, 3, ZZ.rounded_div(22, &7));
742    assert_el_eq!(ZZ, -3, ZZ.rounded_div(-22, &7));
743    assert_el_eq!(ZZ, -3, ZZ.rounded_div(22, &-7));
744    assert_el_eq!(ZZ, 3, ZZ.rounded_div(-22, &-7));
745}
746
747#[test]
748fn test_binomial() {
749    let ZZ = StaticRing::<i32>::RING;
750    assert_eq!(0, binomial(-4, &-1, ZZ));
751    assert_eq!(1, binomial(-4, &0, ZZ));
752    assert_eq!(-4, binomial(-4, &1, ZZ));
753    assert_eq!(10, binomial(-4, &2, ZZ));
754    assert_eq!(-20, binomial(-4, &3, ZZ));
755    assert_eq!(35, binomial(-4, &4, ZZ));
756    assert_eq!(-56, binomial(-4, &5, ZZ));
757
758    assert_eq!(0, binomial(3, &-1, ZZ));
759    assert_eq!(1, binomial(3, &0, ZZ));
760    assert_eq!(3, binomial(3, &1, ZZ));
761    assert_eq!(3, binomial(3, &2, ZZ));
762    assert_eq!(1, binomial(3, &3, ZZ));
763    assert_eq!(0, binomial(3, &4, ZZ));
764    
765    assert_eq!(0, binomial(8, &-1, ZZ));
766    assert_eq!(1, binomial(8, &0, ZZ));
767    assert_eq!(8, binomial(8, &1, ZZ));
768    assert_eq!(28, binomial(8, &2, ZZ));
769    assert_eq!(56, binomial(8, &3, ZZ));
770    assert_eq!(70, binomial(8, &4, ZZ));
771
772    // a naive computation would overflow
773    assert_eq!(145422675, binomial(30, &14, ZZ));
774}
775
776#[test]
777fn test_factorial() {
778    let ZZ = StaticRing::<i32>::RING;
779    assert_eq!(1, factorial(&0, ZZ));
780    assert_eq!(1, factorial(&1, ZZ));
781    assert_eq!(2, factorial(&2, ZZ));
782    assert_eq!(6, factorial(&3, ZZ));
783    assert_eq!(24, factorial(&4, ZZ));
784}
785
786#[test]
787fn test_ceil_floor_div() {
788    let ZZ = StaticRing::<i32>::RING;
789    for rhs in [-10, -3, -2, -1, 1, 2, 3, 10] {
790        for lhs in [-10, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 10] {
791            let result = ZZ.ceil_div(lhs, &rhs);
792            assert_eq!(i32::div_ceil(lhs, rhs), result);
793            assert_eq!((lhs as f64 / rhs as f64).ceil() as i32, result);
794
795            let result = ZZ.floor_div(lhs, &rhs);
796            assert_eq!(i32::div_floor(lhs, rhs), result);
797            assert_eq!((lhs as f64 / rhs as f64).floor() as i32, result);
798        }
799    }
800}
801
802#[test]
803fn test_parse() {
804    let ZZbig = BigIntRing::RING;
805    assert_el_eq!(&ZZbig, &ZZbig.int_hom().map(3), ZZbig.parse("3", 10).unwrap());
806    assert_el_eq!(&ZZbig, &ZZbig.power_of_two(100), ZZbig.parse("1267650600228229401496703205376", 10).unwrap());
807    assert_el_eq!(&ZZbig, &ZZbig.power_of_two(100), ZZbig.parse("+1267650600228229401496703205376", 10).unwrap());
808    assert_el_eq!(&ZZbig, &ZZbig.negate(ZZbig.power_of_two(100)), ZZbig.parse("-1267650600228229401496703205376", 10).unwrap());
809    assert_el_eq!(&ZZbig, &ZZbig.mul(ZZbig.pow(ZZbig.int_hom().map(11), 26), ZZbig.int_hom().map(10)), ZZbig.parse("a00000000000000000000000000", 11).unwrap());
810    assert!(ZZbig.parse("238597a", 10).is_err());
811}
812
813#[test]
814fn test_map_from_integer() {
815    let ZZbig = BigIntRing::RING;
816    assert_el_eq!(&ZZbig, ZZbig.parse("0", 10).unwrap(), map_from_integer_ring(ZZbig, ZZbig, ZZbig.parse("0", 10).unwrap()));
817    assert_el_eq!(&ZZbig, ZZbig.parse("-1", 10).unwrap(), map_from_integer_ring(ZZbig, ZZbig, ZZbig.parse("-1", 10).unwrap()));
818    assert_el_eq!(&ZZbig, ZZbig.parse("80934517340527398320561", 10).unwrap(), map_from_integer_ring(ZZbig, ZZbig, ZZbig.parse("80934517340527398320561", 10).unwrap()));
819    assert_el_eq!(&ZZbig, ZZbig.parse("-4861235897136598713465987138952643", 10).unwrap(), map_from_integer_ring(ZZbig, ZZbig, ZZbig.parse("-4861235897136598713465987138952643", 10).unwrap()));
820}
feanor_math/integer.rs

feanor_math/
integer.rs
