feanor_math/algorithms/
int_bisect.rs

1use crate::ring::*;
2use crate::integer::*;
3use crate::ordered::*;
4
5///
6/// Finds some integer `left <= n < right` such that `f(n) <= 0` and `f(n + 1) > 0`, given
7/// that `f(left) <= 0` and `f(right) > 0`.
8/// 
9/// If we consider a continuous extension of `f` to the real numbers, this means that
10/// the function finds `floor(x)` for some root `x` of `f` between `left` and `right`.
11/// 
12pub fn bisect_floor<R, F>(ZZ: R, left: El<R>, right: El<R>, mut func: F) -> El<R>
13    where R: RingStore, R::Type: IntegerRing, F: FnMut(&El<R>) -> El<R>
14{
15    assert!(ZZ.is_lt(&left, &right));
16    let mut l = left;
17    let mut r = right;
18    assert!(!ZZ.is_pos(&func(&l)));
19    assert!(ZZ.is_pos(&func(&r)));
20    loop {
21        let mut mid = ZZ.add_ref(&l, &r);
22        ZZ.euclidean_div_pow_2(&mut mid, 1);
23
24        if ZZ.eq_el(&mid, &l) || ZZ.eq_el(&mid, &r) {
25            return l;
26        } else if ZZ.is_pos(&func(&mid)) {
27            r = mid;
28        } else {
29            l = mid;
30        }
31    }
32}
33
34///
35/// Given a function `f` with `lim_{x -> -inf} f(x) = -inf` and `lim_{x -> inf} f(x) = inf`,
36/// finds some integer `n` such that `f(n) <= 0` and `f(n + 1) > 0`.
37/// 
38/// Good performance is only guaranteed for increasing functions `f`. In this case, the
39/// solution is unique, and the function terminates `O(log|approx - solution|)` steps.
40/// 
41/// If we consider a continuous extension of `f` to the real numbers, this means that
42/// the function finds `floor(x)` for some root `x` of `f`.
43/// 
44pub fn find_root_floor<R, F>(ZZ: R, approx: El<R>, mut func: F) -> El<R>
45    where R: RingStore, R::Type: IntegerRing, F: FnMut(&El<R>) -> El<R>
46{
47    let mut left = ZZ.clone_el(&approx);
48    let mut step = ZZ.one();
49    while ZZ.is_pos(&func(&left)) {
50        ZZ.sub_assign_ref(&mut left, &step);
51        ZZ.mul_pow_2(&mut step, 1);
52    }
53    step = ZZ.one();
54    let mut right = approx;
55    while !ZZ.is_pos(&func(&right)) {
56        ZZ.add_assign_ref(&mut right, &step);
57        ZZ.mul_pow_2(&mut step, 1);
58    }
59    return bisect_floor(ZZ, left, right, func);
60}
61
62///
63/// Computes the largest integer `x` such that `x^root <= n`, where `n` is a positive integer.
64/// 
65/// # Integer overflow
66/// 
67/// This function is designed to avoid integer overflow if possible without a huge performance
68/// benefit. In particular, in some cases the naive way would require evaluating `(x + 1)^root`
69/// to ensure this is really larger than `n` - however, this might lead to integer overflow as in
70/// ```should_panic
71/// // compute the 62nd root of 2^62
72/// let result: i64 = 2;
73/// assert!(result.pow(62) <= (1 << 62) && (result + 1).pow(62) > (1 << 62));
74/// ```
75/// These cases can somewhat be avoided by first doing a size estimate.
76/// ```rust
77/// # use feanor_math::ring::*;
78/// # use feanor_math::algorithms::int_bisect::*;
79/// # use feanor_math::primitive_int::*;
80/// let ZZ = StaticRing::<i64>::RING;
81/// assert_eq!(2, root_floor(&ZZ, 1 << 62, 62));
82/// ```
83/// I some edge cases this is not enough, and so integer overflow can occur.
84/// ```should_panic
85/// # use feanor_math::ring::*;
86/// # use feanor_math::algorithms::int_bisect::*;
87/// # use feanor_math::primitive_int::*;
88/// let ZZ = StaticRing::<i64>::RING;
89/// assert_eq!(26, root_floor(&ZZ, ZZ.pow(5, 26), 26));
90/// ```
91/// 
92pub fn root_floor<R>(ZZ: R, n: El<R>, root: usize) -> El<R>
93    where R: RingStore, R::Type: IntegerRing
94{
95    assert!(root > 0);
96    if root == 1 {
97        return n;
98    }
99    if ZZ.is_zero(&n) {
100        return ZZ.zero();
101    }
102    assert!(!ZZ.is_neg(&n));
103    let log2_ceil_n = ZZ.abs_log2_ceil(&n).unwrap();
104
105    return find_root_floor(
106        &ZZ, 
107        ZZ.from_float_approx(ZZ.to_float_approx(&n).powf(1. / root as f64)).unwrap_or(ZZ.zero()), 
108        |x| {
109            let x_pow_root_half = ZZ.pow(ZZ.clone_el(x), root / 2);
110            // we first make a size estimate, mainly to avoid situations (high `root`) where we get avoidable integer overflow (also might improve performance)
111            if ZZ.abs_log2_ceil(x).unwrap_or(0) >= 2 && (ZZ.abs_log2_ceil(&x_pow_root_half).unwrap() - 1) * 2 >= log2_ceil_n {
112                if ZZ.is_neg(x) {
113                    return ZZ.negate(ZZ.clone_el(&n));
114                } else {
115                    return ZZ.clone_el(&n);
116                }
117            } else {
118                let x_pow_root = if root % 2 == 0 {
119                    ZZ.pow(x_pow_root_half, 2)
120                } else {
121                    ZZ.mul_ref_snd(ZZ.pow(x_pow_root_half, 2), x)
122                };
123                return ZZ.sub_ref_snd(x_pow_root, &n);
124            }
125        }
126    );
127}
128
129#[cfg(test)]
130use crate::primitive_int::StaticRing;
131
132#[test]
133fn test_bisect_floor() {
134    assert_eq!(0, bisect_floor(&StaticRing::<i64>::RING, 0, 10, |x| if *x == 0 { 0 } else { 1 }));
135    assert_eq!(9, bisect_floor(&StaticRing::<i64>::RING, 0, 10, |x| if *x == 10 { 1 } else { 0 }));
136    assert_eq!(-15, bisect_floor(&StaticRing::<i64>::RING, -20, -10, |x| *x + 15));
137}
138
139#[test]
140fn test_root_floor() {
141    assert_eq!(4, root_floor(&StaticRing::<i64>::RING, 16, 2));
142    assert_eq!(3, root_floor(&StaticRing::<i64>::RING, 27, 3));
143    assert_eq!(4, root_floor(&StaticRing::<i64>::RING, 17, 2));
144    assert_eq!(3, root_floor(&StaticRing::<i64>::RING, 28, 3));
145    assert_eq!(4, root_floor(&StaticRing::<i64>::RING, 24, 2));
146    assert_eq!(3, root_floor(&StaticRing::<i64>::RING, 63, 3));
147    assert_eq!(5, root_floor(&StaticRing::<i64>::RING, StaticRing::<i64>::RING.pow(5, 25), 25));
148    assert_eq!(4, root_floor(&StaticRing::<i64>::RING, StaticRing::<i64>::RING.pow(5, 25), 26));
149    assert_eq!(4, root_floor(&StaticRing::<i64>::RING, StaticRing::<i64>::RING.pow(5, 25), 27));
150    assert_eq!(4, root_floor(&StaticRing::<i64>::RING, StaticRing::<i64>::RING.pow(5, 25), 28));
151}
feanor_math/algorithms/int_bisect.rs

feanor_math/algorithms/
int_bisect.rs
