feanor_math/
integer.rs

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