sirius/solver/poly/
interval.rs

1//! Interval arithmetic with [`Poly`] endpoints. Allows attempted interval tracking, given non-negative variables, for e.g. loop bounds.
2//!
3//! | operation | result                              | requires             |
4//! | --------- | ----------------------------------- | -------------------- |
5//! | $x + y$   | $[a + c, b + d]$                    |                      |
6//! | $-x$      | $[-b, -a]$                          |                      |
7//! | $px$      | $[pa, pb]$                          | $p \geq 0$           |
8//! | $xy$      | $[ac, bd]$                          | $a \geq 0, c \geq 0$ |
9//! | $x^n$     | $[a^n, b^n]$                        | $n$ odd              |
10//! | $x^n$     | $[a^n, b^n]$                        | $n$ even, $a \geq 0$ |
11//! | $x^n$     | $[0, a^n + b^n]$                    | $n$ even             |
12//!
13//! ```
14//! use sirius::solver::poly::{poly, Var, interval::Interval};
15//!
16//! let nonneg = |_: Var| true;
17//! let i = Interval::new(poly!(0), poly!(a - 1));
18//! let j = Interval::new(poly!(0), poly!(b - 1));
19//!
20//! let range = poly!(i*b + j)
21//!     .substitute_interval('j' as Var, &j, &nonneg)
22//!     .and_then(|r| r.substitute('i' as Var, &i, &nonneg))
23//!     .unwrap();
24//! assert_eq!(range, Interval::new(poly!(0), poly!(a*b - 1)));
25//! ```
26
27use super::Poly;
28use super::coef::Coef;
29use super::mono::{Mono, Pow, Var};
30
31/// A closed range $[\text{lo}, \text{hi}]$ with [`Poly`] endpoints.
32#[derive(Clone, PartialEq, Eq, Hash, Debug)]
33pub struct Interval {
34    pub lo: Poly,
35    pub hi: Poly,
36}
37
38impl Interval {
39    pub fn new(lo: Poly, hi: Poly) -> Self {
40        Self { lo, hi }
41    }
42
43    /// The range of a value known exactly, `[p, p]`.
44    pub fn exact(p: Poly) -> Self {
45        Self {
46            lo: p.clone(),
47            hi: p,
48        }
49    }
50
51    pub fn is_exact(&self) -> bool {
52        self.lo == self.hi
53    }
54
55    /// Every value in the range is non-negative, given non-negative variables.
56    pub fn is_nonneg(&self, nonneg: &dyn Fn(Var) -> bool) -> bool {
57        self.lo.always_nonneg(nonneg)
58    }
59
60    /// Every value in the range is non-positive, given non-negative variables.
61    pub fn is_nonpos(&self, nonneg: &dyn Fn(Var) -> bool) -> bool {
62        self.hi.neg().always_nonneg(nonneg)
63    }
64
65    pub fn add(&self, rhs: &Self) -> Self {
66        Self {
67            lo: self.lo.add(&rhs.lo),
68            hi: self.hi.add(&rhs.hi),
69        }
70    }
71
72    pub fn neg(&self) -> Self {
73        Self {
74            lo: self.hi.neg(),
75            hi: self.lo.neg(),
76        }
77    }
78
79    pub fn sub(&self, rhs: &Self) -> Self {
80        self.add(&rhs.neg())
81    }
82
83    pub fn mul_scalar<T: Into<Coef>>(&self, c: T) -> Self {
84        let c = c.into();
85        if c.is_positive() || c.is_zero() {
86            Self {
87                lo: self.lo.mul_scalar(c),
88                hi: self.hi.mul_scalar(c),
89            }
90        } else {
91            Self {
92                lo: self.hi.mul_scalar(c),
93                hi: self.lo.mul_scalar(c),
94            }
95        }
96    }
97
98    /// Multiply by a polynomial whose sign is known.
99    pub fn scale(&self, p: &Poly, nonneg: &dyn Fn(Var) -> bool) -> Option<Self> {
100        if p.always_nonneg(nonneg) {
101            Some(Self {
102                lo: self.lo.mul(p),
103                hi: self.hi.mul(p),
104            })
105        } else if p.neg().always_nonneg(nonneg) {
106            Some(Self {
107                lo: self.hi.mul(p),
108                hi: self.lo.mul(p),
109            })
110        } else {
111            None
112        }
113    }
114
115    /// The product of two ranges, when the sign of each is known. The product of two exact
116    /// values is exact regardless.
117    pub fn mul(&self, rhs: &Self, nonneg: &dyn Fn(Var) -> bool) -> Option<Self> {
118        if self.is_exact() {
119            return if rhs.is_exact() {
120                Some(Self::exact(self.lo.mul(&rhs.lo)))
121            } else {
122                rhs.scale(&self.lo, nonneg)
123            };
124        }
125        if rhs.is_exact() {
126            return self.scale(&rhs.lo, nonneg);
127        }
128
129        // x in [a, b], y in [c, d]; each case is the all-non-negative one under a sign flip
130        let (a, b, c, d) = (&self.lo, &self.hi, &rhs.lo, &rhs.hi);
131        let (lo, hi) = match (
132            self.is_nonneg(nonneg),
133            self.is_nonpos(nonneg),
134            rhs.is_nonneg(nonneg),
135            rhs.is_nonpos(nonneg),
136        ) {
137            (true, _, true, _) => (a.mul(c), b.mul(d)),
138            (true, _, _, true) => (b.mul(c), a.mul(d)),
139            (_, true, true, _) => (a.mul(d), b.mul(c)),
140            (_, true, _, true) => (b.mul(d), a.mul(c)),
141            _ => return None,
142        };
143        Some(Self { lo, hi })
144    }
145
146    pub fn pow(&self, e: Pow, nonneg: &dyn Fn(Var) -> bool) -> Self {
147        if e == 0 {
148            return Self::exact(Poly::constant(1));
149        }
150        if e == 1 || self.is_exact() {
151            return Self {
152                lo: self.lo.pow(e),
153                hi: self.hi.pow(e),
154            };
155        }
156        // odd powers are monotone everywhere; even ones only away from zero
157        let (lo, hi) = if e % 2 == 1 || self.is_nonneg(nonneg) {
158            (self.lo.pow(e), self.hi.pow(e))
159        } else if self.is_nonpos(nonneg) {
160            (self.hi.pow(e), self.lo.pow(e))
161        } else {
162            // |x| <= max(|lo|, |hi|), and max(|lo|, |hi|)^e <= lo^e + hi^e for even e
163            (Poly::zero(), self.lo.pow(e).add(&self.hi.pow(e)))
164        };
165        Self { lo, hi }
166    }
167
168    /// The range as `v` also ranges over `range`: the endpoints widen outward.
169    pub fn substitute(
170        &self,
171        v: Var,
172        range: &Interval,
173        nonneg: &dyn Fn(Var) -> bool,
174    ) -> Option<Self> {
175        Some(Self {
176            lo: self.lo.substitute_interval(v, range, nonneg)?.lo,
177            hi: self.hi.substitute_interval(v, range, nonneg)?.hi,
178        })
179    }
180}
181
182impl Poly {
183    /// The range of `self` as `v` ranges over `range`, every other variable held fixed. Each term
184    /// `c * m * v^e` is bounded as `(c * m) * range^e`, so the sign of `c * m` has to be known.
185    pub fn substitute_interval(
186        &self,
187        v: Var,
188        range: &Interval,
189        nonneg: &dyn Fn(Var) -> bool,
190    ) -> Option<Interval> {
191        let mut acc = Interval::exact(Poly::zero());
192        for (c, m) in self.terms() {
193            let e = m.degree_in(v);
194            let term = if e == 0 {
195                Interval::exact(Poly::term(c, m))
196            } else {
197                let rest = m
198                    .div(&Mono::new([(v, e)]))
199                    .expect("monomial divides by its own factor");
200                range.pow(e, nonneg).scale(&Poly::term(c, rest), nonneg)?
201            };
202            acc = acc.add(&term);
203        }
204        Some(acc)
205    }
206}
207
208#[cfg(test)]
209mod tests {
210    use super::{Interval, Poly, Var};
211    use crate::solver::poly::poly;
212
213    fn v(name: char) -> Var {
214        name as Var
215    }
216
217    fn all(_: Var) -> bool {
218        true
219    }
220
221    /// Everything but `a`.
222    fn all_but_a(x: Var) -> bool {
223        x != v('a')
224    }
225
226    fn c(k: i128) -> Poly {
227        Poly::constant(k)
228    }
229
230    fn iv(lo: i128, hi: i128) -> Interval {
231        Interval::new(c(lo), c(hi))
232    }
233
234    fn none(_: Var) -> bool {
235        false
236    }
237
238    #[test]
239    fn row_major_access() {
240        let i = Interval::new(poly!(0), poly!(a - 1));
241        let j = Interval::new(poly!(0), poly!(b - 1));
242        let k = Interval::new(poly!(0), poly!(c - 1));
243
244        let index = poly!(i * b * c + j * c + k);
245        let range = index
246            .substitute_interval(v('k'), &k, &all)
247            .and_then(|r| r.substitute(v('j'), &j, &all))
248            .and_then(|r| r.substitute(v('i'), &i, &all))
249            .unwrap();
250        assert_eq!(range, Interval::new(poly!(0), poly!(a * b * c - 1)));
251
252        // a reversed access
253        let range = poly!(a * b - 1 - i * b - j)
254            .substitute_interval(v('j'), &j, &all)
255            .and_then(|r| r.substitute(v('i'), &i, &all))
256            .unwrap();
257        assert_eq!(range, Interval::new(poly!(0), poly!(a * b - 1)));
258    }
259
260    #[test]
261    fn correlated_bounds_survive_one_substitution() {
262        // for i from 0 to N: for j from i to N: j - i
263        let j = Interval::new(poly!(i), poly!(n - 1));
264        let range = poly!(j - i).substitute_interval(v('j'), &j, &all).unwrap();
265        assert_eq!(range, Interval::new(poly!(0), poly!(n - 1 - i)));
266    }
267
268    #[test]
269    fn substitution_needs_a_coefficient_sign() {
270        let i = Interval::new(poly!(0), poly!(n - 1));
271        // `a` is not known non-negative, so `a*i` cannot be bounded
272        assert_eq!(
273            poly!(a * i).substitute_interval(v('i'), &i, &all_but_a),
274            None
275        );
276        // but `a^2*i` can
277        assert_eq!(
278            poly!(a ^ 2 * i).substitute_interval(v('i'), &i, &all_but_a),
279            Some(Interval::new(poly!(0), poly!(a ^ 2 * n - a ^ 2)))
280        );
281        // a variable absent from the polynomial changes nothing
282        assert_eq!(
283            poly!(a + 1).substitute_interval(v('i'), &i, &all_but_a),
284            Some(Interval::exact(poly!(a + 1)))
285        );
286    }
287
288    #[test]
289    fn products_by_sign() {
290        // both non-negative
291        assert_eq!(iv(1, 2).mul(&iv(3, 4), &none), Some(iv(3, 8)));
292        // one non-positive
293        assert_eq!(iv(1, 2).mul(&iv(-4, -3), &none), Some(iv(-8, -3)));
294        assert_eq!(iv(-4, -3).mul(&iv(1, 2), &none), Some(iv(-8, -3)));
295        // both non-positive
296        assert_eq!(iv(-2, -1).mul(&iv(-4, -3), &none), Some(iv(3, 8)));
297        // straddling zero
298        assert_eq!(iv(-1, 2).mul(&iv(3, 4), &none), None);
299        // exact factors scale
300        assert_eq!(
301            Interval::exact(c(-2)).mul(&iv(-1, 3), &none),
302            Some(iv(-6, 2))
303        );
304        assert_eq!(
305            iv(-1, 3).mul(&Interval::exact(c(-2)), &none),
306            Some(iv(-6, 2))
307        );
308        assert_eq!(
309            Interval::exact(c(-2)).mul(&Interval::exact(c(-3)), &none),
310            Some(Interval::exact(c(6)))
311        );
312
313        // symbolic: `n` is non-negative, `a` is not
314        let n = Interval::exact(poly!(n));
315        assert_eq!(
316            n.mul(&Interval::new(poly!(a), poly!(a + 1)), &all_but_a),
317            Some(Interval::new(poly!(a * n), poly!(a * n + n)))
318        );
319        assert_eq!(
320            Interval::exact(poly!(a)).mul(&Interval::new(poly!(0), poly!(n)), &all_but_a),
321            None
322        );
323    }
324
325    #[test]
326    fn powers_by_sign() {
327        assert_eq!(iv(2, 3).pow(2, &none), iv(4, 9));
328        assert_eq!(iv(-3, -2).pow(2, &none), iv(4, 9));
329        assert_eq!(iv(-3, -2).pow(3, &none), iv(-27, -8));
330        assert_eq!(iv(-3, 2).pow(3, &none), iv(-27, 8));
331        // straddling zero: a loose but sound bound
332        assert_eq!(iv(-3, 2).pow(2, &none), iv(0, 13));
333        assert_eq!(iv(-3, 2).pow(0, &none), Interval::exact(c(1)));
334        assert_eq!(
335            Interval::exact(poly!(a)).pow(2, &all_but_a),
336            Interval::exact(poly!(a ^ 2))
337        );
338    }
339
340    /// Random polynomials over nested symbolic ranges, checked at random points.
341    #[test]
342    fn substitution_fuzz() {
343        use rand::prelude::*;
344
345        let mut rng = SmallRng::seed_from_u64(7);
346
347        // `n`, `m` are non-negative typevars; `a` is a rigid var of either sign
348        let rigid = [v('n'), v('m'), v('a')];
349        let nonneg_base = |x: Var| x == v('n') || x == v('m');
350
351        fn random_poly(rng: &mut SmallRng, vars: &[Var], terms: usize, max_pow: u32) -> Poly {
352            let mut p = Poly::zero();
353            for _ in 0..rng.gen_range(0..=terms) {
354                let mut t = Poly::constant(rng.gen_range(-4..=4));
355                for &x in vars {
356                    t = t.mul(&Poly::var(x, rng.gen_range(0..=max_pow)));
357                }
358                p = p.add(&t);
359            }
360            p
361        }
362
363        let mut bounded = 0;
364        let mut checked = 0;
365
366        for _ in 0..3000 {
367            // i ranges over polys in the rigid vars; j over those plus i
368            let i_lo = random_poly(&mut rng, &rigid, 2, 1);
369            let i_hi = i_lo.add(&random_poly(&mut rng, &rigid, 2, 1));
370            let i_range = Interval::new(i_lo, i_hi);
371
372            let with_i = [v('n'), v('m'), v('a'), v('i')];
373            let j_lo = random_poly(&mut rng, &with_i, 2, 1);
374            let j_hi = j_lo.add(&random_poly(&mut rng, &with_i, 2, 1));
375            let j_range = Interval::new(j_lo, j_hi);
376
377            // a var is non-negative only when its lower bound is
378            let i_nonneg = i_range.lo.always_nonneg(&nonneg_base);
379            let nonneg_i = |x: Var| nonneg_base(x) || (x == v('i') && i_nonneg);
380            let j_nonneg = j_range.lo.always_nonneg(&nonneg_i);
381            let nonneg = |x: Var| nonneg_i(x) || (x == v('j') && j_nonneg);
382
383            let all_vars = [v('n'), v('m'), v('a'), v('i'), v('j')];
384            let p = random_poly(&mut rng, &all_vars, 4, 2);
385
386            // innermost first, as the typechecker does
387            let Some(range) = p
388                .substitute_interval(v('j'), &j_range, &nonneg)
389                .and_then(|r| r.substitute(v('i'), &i_range, &nonneg))
390            else {
391                continue;
392            };
393            bounded += 1;
394
395            for &x in range.lo.vars().iter().chain(range.hi.vars().iter()) {
396                assert!(rigid.contains(&x), "range mentions a substituted var");
397            }
398
399            for _ in 0..20 {
400                let n = rng.gen_range(0..=5);
401                let m = rng.gen_range(0..=5);
402                let a = rng.gen_range(-5..=5);
403                let at_rigid = |x: Var| match x {
404                    _ if x == v('n') => n,
405                    _ if x == v('m') => m,
406                    _ if x == v('a') => a,
407                    _ => panic!("not rigid"),
408                };
409                let (i_lo, i_hi) = (
410                    i_range.lo.eval(at_rigid).get(),
411                    i_range.hi.eval(at_rigid).get(),
412                );
413                if i_lo > i_hi {
414                    continue;
415                }
416                let i = rng.gen_range(i_lo..=i_hi);
417                assert!(!i_nonneg || i >= 0);
418                let at_i = |x: Var| if x == v('i') { i } else { at_rigid(x) };
419                let (j_lo, j_hi) = (j_range.lo.eval(at_i).get(), j_range.hi.eval(at_i).get());
420                if j_lo > j_hi {
421                    continue;
422                }
423                let j = rng.gen_range(j_lo..=j_hi);
424                assert!(!j_nonneg || j >= 0);
425                let at_j = |x: Var| if x == v('j') { j } else { at_i(x) };
426
427                let value = p.eval(at_j).get();
428                let (lo, hi) = (range.lo.eval(at_rigid).get(), range.hi.eval(at_rigid).get());
429                assert!(
430                    lo <= value && value <= hi,
431                    "{p:?} = {value} at n={n} m={m} a={a} i={i} j={j} escapes {range:?} = [{lo}, {hi}] \
432                     with i in {i_range:?}, j in {j_range:?}"
433                );
434                checked += 1;
435            }
436        }
437
438        assert!(bounded > 500, "only {bounded} polynomials were bounded");
439        assert!(checked > 5000, "only {checked} points were checked");
440    }
441
442    /// Products and powers of ranges with symbolic endpoints, checked at random points.
443    #[test]
444    fn product_fuzz() {
445        use rand::prelude::*;
446
447        let mut rng = SmallRng::seed_from_u64(11);
448        let rigid = [v('n'), v('a')];
449        let nonneg = |x: Var| x == v('n');
450
451        fn random_poly(rng: &mut SmallRng, vars: &[Var]) -> Poly {
452            let mut p = Poly::zero();
453            for _ in 0..rng.gen_range(0..=2) {
454                let mut t = Poly::constant(rng.gen_range(-3..=3));
455                for &x in vars {
456                    t = t.mul(&Poly::var(x, rng.gen_range(0..=1)));
457                }
458                p = p.add(&t);
459            }
460            p
461        }
462
463        let mut products = 0;
464        let mut powers = 0;
465
466        for _ in 0..3000 {
467            let x = Interval::new(random_poly(&mut rng, &rigid), random_poly(&mut rng, &rigid));
468            let y = Interval::new(random_poly(&mut rng, &rigid), random_poly(&mut rng, &rigid));
469            let e = rng.gen_range(0..=4);
470
471            let product = x.mul(&y, &nonneg);
472            let power = x.pow(e, &nonneg);
473            products += product.is_some() as usize;
474
475            for _ in 0..20 {
476                let n = rng.gen_range(0..=4);
477                let a = rng.gen_range(-4..=4);
478                let at = |var: Var| if var == v('n') { n } else { a };
479                let (x_lo, x_hi) = (x.lo.eval(at).get(), x.hi.eval(at).get());
480                let (y_lo, y_hi) = (y.lo.eval(at).get(), y.hi.eval(at).get());
481                if x_lo > x_hi {
482                    continue;
483                }
484                let xv = rng.gen_range(x_lo..=x_hi);
485
486                let (lo, hi) = (power.lo.eval(at).get(), power.hi.eval(at).get());
487                let value = xv.pow(e);
488                assert!(
489                    lo <= value && value <= hi,
490                    "{xv}^{e} escapes {power:?} = [{lo}, {hi}] with x in {x:?} at n={n} a={a}"
491                );
492                powers += 1;
493
494                let Some(product) = &product else {
495                    continue;
496                };
497                if y_lo > y_hi {
498                    continue;
499                }
500                let yv = rng.gen_range(y_lo..=y_hi);
501                let (lo, hi) = (product.lo.eval(at).get(), product.hi.eval(at).get());
502                assert!(
503                    lo <= xv * yv && xv * yv <= hi,
504                    "{xv}*{yv} escapes {product:?} = [{lo}, {hi}] with x in {x:?}, y in {y:?} \
505                     at n={n} a={a}"
506                );
507            }
508        }
509
510        assert!(products > 300, "only {products} products were bounded");
511        assert!(powers > 5000, "only {powers} powers were checked");
512    }
513}
sirius/solver/poly/interval.rs

sirius/solver/poly/
interval.rs
