axiolid_exact/
poly.rs

1//! Integer polynomials and exact real-root isolation.
2//!
3//! Conic intersections lead to a quartic whose roots are nested radicals
4//! at best and not radicals at all in general (casus irreducibilis). So an
5//! exact root is represented the way CGAL's `Algebraic_kernel_d` does it:
6//! a square-free integer polynomial plus a dyadic interval that contains
7//! exactly one of its roots. Every comparison is then decided exactly:
8//! by Sturm counts, by refining the intervals, or by a gcd when two roots
9//! might be equal.
10//!
11//! Coefficients are `BigInt`. A polynomial with dyadic coefficients is
12//! scaled by a power of two first ([`IntPoly::from_dyadic`]), which does
13//! not move its roots.
14
15use std::cmp::Ordering;
16
17use num_bigint::BigInt;
18use num_integer::Integer;
19
20use axiolid_guarantees::Sign;
21
22use crate::arith::Arith;
23use crate::dyadic::Dyadic;
24
25/// A polynomial with integer coefficients, lowest degree first, no
26/// trailing zeros (the zero polynomial is empty).
27#[derive(Debug, Clone, PartialEq, Eq, Hash)]
28pub struct IntPoly {
29    coeffs: Vec<BigInt>,
30}
31
32impl IntPoly {
33    /// From coefficients, lowest degree first. Trailing zeros are dropped.
34    #[must_use]
35    pub fn new(mut coeffs: Vec<BigInt>) -> Self {
36        while coeffs
37            .last()
38            .is_some_and(|c| c.sign() == num_bigint::Sign::NoSign)
39        {
40            coeffs.pop();
41        }
42        Self { coeffs }
43    }
44
45    /// From dyadic coefficients, lowest degree first, scaled by a power of
46    /// two to integers. The roots are those of the dyadic polynomial.
47    #[must_use]
48    pub fn from_dyadic(coeffs: &[Dyadic]) -> Self {
49        let lowest = coeffs
50            .iter()
51            .filter(|c| c.sign() != Some(Sign::Zero))
52            .map(Dyadic::exponent)
53            .min()
54            .unwrap_or(0);
55        Self::new(
56            coeffs
57                .iter()
58                .map(|c| {
59                    if c.sign() == Some(Sign::Zero) {
60                        BigInt::from(0)
61                    } else {
62                        c.mantissa() << (c.exponent() - lowest) as u64
63                    }
64                })
65                .collect(),
66        )
67    }
68
69    /// Coefficients, lowest degree first.
70    #[must_use]
71    pub fn coeffs(&self) -> &[BigInt] {
72        &self.coeffs
73    }
74
75    /// Degree, or `None` for the zero polynomial.
76    #[must_use]
77    pub fn degree(&self) -> Option<usize> {
78        self.coeffs.len().checked_sub(1)
79    }
80
81    /// Whether this is the zero polynomial.
82    #[must_use]
83    pub fn is_zero(&self) -> bool {
84        self.coeffs.is_empty()
85    }
86
87    fn lead(&self) -> &BigInt {
88        self.coeffs.last().expect("non-zero polynomial")
89    }
90
91    /// Exact value at a dyadic point.
92    #[must_use]
93    pub fn eval(&self, x: &Dyadic) -> Dyadic {
94        // Horner, exact.
95        let mut acc = Dyadic::zero();
96        for c in self.coeffs.iter().rev() {
97            acc = acc.mul(x).add(&Dyadic::from_parts(c.clone(), 0));
98        }
99        acc
100    }
101
102    /// Exact sign at a dyadic point.
103    #[must_use]
104    pub fn sign_at(&self, x: &Dyadic) -> Sign {
105        self.eval(x).sign().expect("exact")
106    }
107
108    /// The derivative.
109    #[must_use]
110    pub fn derivative(&self) -> Self {
111        Self::new(
112            self.coeffs
113                .iter()
114                .enumerate()
115                .skip(1)
116                .map(|(i, c)| c * BigInt::from(i))
117                .collect(),
118        )
119    }
120
121    fn content(&self) -> BigInt {
122        self.coeffs.iter().fold(BigInt::from(0), |g, c| g.gcd(c))
123    }
124
125    /// Divided by the gcd of its coefficients, leading coefficient
126    /// positive. Keeps coefficient growth in Sturm chains in check.
127    #[must_use]
128    pub fn primitive(&self) -> Self {
129        if self.is_zero() {
130            return self.clone();
131        }
132        let mut g = self.content();
133        if self.lead().sign() == num_bigint::Sign::Minus {
134            g = -g;
135        }
136        Self::new(self.coeffs.iter().map(|c| c / &g).collect())
137    }
138
139    /// Divided by the absolute gcd of its coefficients: same signs, same
140    /// roots, smaller numbers. Sturm chains need this rather than
141    /// [`Self::primitive`], which may flip the sign.
142    fn scaled_down(&self) -> Self {
143        if self.is_zero() {
144            return self.clone();
145        }
146        let g = self.content();
147        let g = if g.sign() == num_bigint::Sign::Minus {
148            -g
149        } else {
150            g
151        };
152        Self::new(self.coeffs.iter().map(|c| c / &g).collect())
153    }
154
155    /// Pseudo-remainder of `self` by `divisor`: the remainder of
156    /// `lead(divisor)^k * self` by `divisor`, with `k` just large enough to
157    /// keep all arithmetic in the integers. Positive leading factor, so the
158    /// sign structure a Sturm chain needs is preserved.
159    fn pseudo_rem(&self, divisor: &Self) -> Self {
160        let ld = divisor.lead().clone();
161        let dd = divisor.degree().expect("non-zero divisor");
162        let mut r = self.coeffs.clone();
163        // lead(divisor) squared is positive, so multiply by |lead| as needed
164        // while keeping the sign right: using ld^2 per step when ld < 0
165        // would also work, but scaling by |ld| suffices because we only
166        // eliminate the top term and the factor is positive.
167        let scale = if ld.sign() == num_bigint::Sign::Minus {
168            -ld.clone()
169        } else {
170            ld.clone()
171        };
172        while r.len() > dd && !r.is_empty() {
173            let top = r.last().cloned().unwrap_or_default();
174            if top.sign() == num_bigint::Sign::NoSign {
175                r.pop();
176                continue;
177            }
178            let shift = r.len() - 1 - dd;
179            // r := scale * r - (top * sign(ld)) * x^shift * divisor
180            let factor = if ld.sign() == num_bigint::Sign::Minus {
181                -top
182            } else {
183                top
184            };
185            for c in r.iter_mut() {
186                *c *= &scale;
187            }
188            for (i, dc) in divisor.coeffs.iter().enumerate() {
189                r[shift + i] -= &factor * dc;
190            }
191            r.pop();
192        }
193        Self::new(r)
194    }
195
196    fn gcd_poly(&self, other: &Self) -> Self {
197        let (mut a, mut b) = (self.primitive(), other.primitive());
198        while !b.is_zero() {
199            let r = a.pseudo_rem(&b).primitive();
200            a = b;
201            b = r;
202        }
203        a.primitive()
204    }
205
206    pub(crate) fn exact_div(&self, divisor: &Self) -> Self {
207        // Polynomial long division known to be exact over the rationals;
208        // done over the integers after scaling, then made primitive.
209        let dd = divisor.degree().expect("non-zero divisor");
210        let Some(n) = self.degree() else {
211            return self.clone();
212        };
213        if n < dd {
214            return Self::new(vec![]);
215        }
216        let ld = divisor.lead().clone();
217        let mut r = self.coeffs.clone();
218        let mut q = vec![BigInt::from(0); n - dd + 1];
219        // Scale the dividend so each step divides exactly.
220        let steps = (n - dd + 1) as u32;
221        let scale = num_traits_pow(&ld, steps);
222        for c in r.iter_mut() {
223            *c *= &scale;
224        }
225        for k in (0..=n - dd).rev() {
226            let top = r[k + dd].clone();
227            let coef = &top / &ld;
228            debug_assert_eq!(&coef * &ld, top, "exact division");
229            for (i, dc) in divisor.coeffs.iter().enumerate() {
230                r[k + i] -= &coef * dc;
231            }
232            q[k] = coef;
233        }
234        Self::new(q).primitive()
235    }
236
237    /// The square-free part: same real roots, each with multiplicity one.
238    #[must_use]
239    pub fn square_free(&self) -> Self {
240        if self.degree().unwrap_or(0) < 1 {
241            return self.primitive();
242        }
243        let g = self.gcd_poly(&self.derivative());
244        if g.degree() == Some(0) {
245            self.primitive()
246        } else {
247            self.exact_div(&g)
248        }
249    }
250
251    /// Greatest common divisor, primitive with positive leading
252    /// coefficient.
253    #[must_use]
254    pub fn gcd(&self, other: &Self) -> Self {
255        if self.is_zero() {
256            return other.primitive();
257        }
258        if other.is_zero() {
259            return self.primitive();
260        }
261        self.gcd_poly(other)
262    }
263
264    /// The Sturm chain of a square-free polynomial.
265    fn sturm(&self) -> Vec<Self> {
266        let mut chain = vec![self.clone(), self.derivative()];
267        loop {
268            let n = chain.len();
269            if chain[n - 1].is_zero() {
270                chain.pop();
271                break;
272            }
273            // -rem(p_{k-2}, p_{k-1}); a positive scaling keeps signs.
274            let r = chain[n - 2].pseudo_rem(&chain[n - 1]);
275            if r.is_zero() {
276                break;
277            }
278            // Sturm needs exactly -rem up to a POSITIVE factor.
279            let neg = Self::new(r.scaled_down().coeffs.iter().map(|c| -c).collect());
280            chain.push(neg);
281        }
282        chain
283    }
284
285    /// A bound `B` (a power of two) with every real root in `(-B, B)`.
286    fn root_bound(&self) -> Dyadic {
287        // Cauchy: 1 + max |c_i / c_n|, rounded up to a power of two.
288        let lead_bits = self.lead().bits();
289        let max_bits = self.coeffs.iter().map(BigInt::bits).max().unwrap_or(0);
290        let shift = max_bits.saturating_sub(lead_bits) + 2;
291        Dyadic::from_parts(BigInt::from(1), shift as i64)
292    }
293
294    /// Isolate every real root of this polynomial.
295    ///
296    /// Returns roots in increasing order, each as a [`RealRoot`] whose
297    /// interval contains exactly that root. Multiple roots are reported
298    /// once. The zero polynomial has no isolated roots (it vanishes
299    /// everywhere); callers must handle it before asking.
300    #[must_use]
301    pub fn real_roots(&self) -> Vec<RealRoot> {
302        if self.degree().unwrap_or(0) < 1 {
303            return Vec::new();
304        }
305        let sf = self.square_free();
306        let chain = sf.sturm();
307        let bound = sf.root_bound();
308        let mut out = Vec::new();
309        isolate(&sf, &chain, bound.neg(), bound, &mut out);
310        out
311    }
312}
313
314fn num_traits_pow(base: &BigInt, exp: u32) -> BigInt {
315    let mut out = BigInt::from(1);
316    for _ in 0..exp {
317        out *= base;
318    }
319    out
320}
321
322/// Sign changes of a Sturm chain at `x`, zeros skipped.
323fn variations(chain: &[IntPoly], x: &Dyadic) -> usize {
324    let mut count = 0;
325    let mut last = Sign::Zero;
326    for p in chain {
327        let s = p.sign_at(x);
328        if s == Sign::Zero {
329            continue;
330        }
331        if last != Sign::Zero && s != last {
332            count += 1;
333        }
334        last = s;
335    }
336    count
337}
338
339/// Roots in `(lo, hi]` by Sturm's theorem.
340fn roots_in(chain: &[IntPoly], lo: &Dyadic, hi: &Dyadic) -> usize {
341    variations(chain, lo) - variations(chain, hi)
342}
343
344fn midpoint(lo: &Dyadic, hi: &Dyadic) -> Dyadic {
345    lo.add(hi).mul(&Dyadic::from_parts(BigInt::from(1), -1))
346}
347
348fn isolate(poly: &IntPoly, chain: &[IntPoly], lo: Dyadic, hi: Dyadic, out: &mut Vec<RealRoot>) {
349    let count = roots_in(chain, &lo, &hi);
350    // An open interval must not start at a root: `lo` can be the previous
351    // interval's right end, and a root there. Bisect until it does not.
352    let lo_is_root = poly.sign_at(&lo) == Sign::Zero;
353    match count {
354        0 => {}
355        1 if !lo_is_root || poly.sign_at(&hi) == Sign::Zero => {
356            // Exactly one root in (lo, hi]. If it is hi itself, record it as
357            // an exact point, so intervals stay open at both ends otherwise.
358            if poly.sign_at(&hi) == Sign::Zero {
359                out.push(RealRoot::exact(poly.clone(), hi));
360            } else {
361                out.push(RealRoot {
362                    poly: poly.clone(),
363                    lo,
364                    hi,
365                });
366            }
367        }
368        _ => {
369            let mid = midpoint(&lo, &hi);
370            isolate(poly, chain, lo, mid.clone(), out);
371            isolate(poly, chain, mid, hi, out);
372        }
373    }
374}
375
376/// One real root of a square-free integer polynomial.
377///
378/// Either an exact dyadic point (`lo == hi`), or the unique root in the
379/// open interval `(lo, hi)` with the polynomial non-zero at both ends.
380#[derive(Debug, Clone, PartialEq, Eq)]
381pub struct RealRoot {
382    poly: IntPoly,
383    lo: Dyadic,
384    hi: Dyadic,
385}
386
387impl RealRoot {
388    fn exact(poly: IntPoly, at: Dyadic) -> Self {
389        Self {
390            poly,
391            lo: at.clone(),
392            hi: at,
393        }
394    }
395
396    /// The defining square-free polynomial.
397    #[must_use]
398    pub fn poly(&self) -> &IntPoly {
399        &self.poly
400    }
401
402    /// Current isolating interval `[lo, hi]`.
403    #[must_use]
404    pub fn bounds(&self) -> (&Dyadic, &Dyadic) {
405        (&self.lo, &self.hi)
406    }
407
408    /// Whether the root is a known dyadic value.
409    #[must_use]
410    pub fn is_exact(&self) -> bool {
411        self.lo == self.hi
412    }
413
414    /// Halve the isolating interval once (or pin the root exactly).
415    pub fn refine(&mut self) {
416        if self.is_exact() {
417            return;
418        }
419        let mid = midpoint(&self.lo, &self.hi);
420        let s_mid = self.poly.sign_at(&mid);
421        if s_mid == Sign::Zero {
422            self.lo = mid.clone();
423            self.hi = mid;
424            return;
425        }
426        if s_mid == self.poly.sign_at(&self.lo) {
427            self.lo = mid;
428        } else {
429            self.hi = mid;
430        }
431    }
432
433    /// Refine until the interval is no wider than `width`, or the root is
434    /// exact.
435    pub fn refine_to_width(&mut self, width: &Dyadic) {
436        while !self.is_exact() && self.hi.sub(&self.lo).sub(width).sign() == Some(Sign::Positive) {
437            self.refine();
438        }
439    }
440
441    /// A double near the root (for output only; never for decisions).
442    #[must_use]
443    pub fn approx(&self) -> f64 {
444        let mut copy = self.clone();
445        // 2^-60 relative is well below f64 resolution for any root the
446        // bound admits; exact roots return immediately.
447        for _ in 0..200 {
448            if copy.is_exact() {
449                break;
450            }
451            let lo = copy.lo.to_f64();
452            let hi = copy.hi.to_f64();
453            if lo == hi || hi.next_down() <= lo {
454                break;
455            }
456            copy.refine();
457        }
458        let (lo, hi) = (copy.lo.to_f64(), copy.hi.to_f64());
459        lo + (hi - lo) / 2.0
460    }
461
462    /// Exact sign of `root - x` for a dyadic `x`.
463    #[must_use]
464    pub fn cmp_dyadic(&self, x: &Dyadic) -> Sign {
465        if self.is_exact() {
466            return self.lo.sub(x).sign().expect("exact");
467        }
468        // x outside the open interval decides immediately.
469        if x.sub(&self.lo).sign() != Some(Sign::Positive) {
470            return Sign::Positive;
471        }
472        if x.sub(&self.hi).sign() != Some(Sign::Negative) {
473            return Sign::Negative;
474        }
475        // x strictly inside: the polynomial's sign at x relative to its
476        // sign at lo says which side the root is on.
477        let s = self.poly.sign_at(x);
478        if s == Sign::Zero {
479            return Sign::Zero;
480        }
481        if s == self.poly.sign_at(&self.lo) {
482            // No sign change between lo and x: root is above x.
483            Sign::Positive
484        } else {
485            Sign::Negative
486        }
487    }
488
489    /// Exact sign of the polynomial `q` at this root.
490    ///
491    /// Zero is decided through `gcd(poly, q)`: `q` vanishes at the root iff
492    /// the gcd has a root in the isolating interval, which holds no other
493    /// root of `poly`. Otherwise the interval is refined until `q` has no
494    /// root in it, and `q`'s sign there is its sign at the root.
495    #[must_use]
496    pub fn sign_of(&self, q: &IntPoly) -> Sign {
497        if q.is_zero() {
498            return Sign::Zero;
499        }
500        if self.is_exact() {
501            return q.sign_at(&self.lo);
502        }
503        let g = self.poly.gcd(q);
504        if g.degree().unwrap_or(0) >= 1 && roots_in(&g.sturm(), &self.lo, &self.hi) > 0 {
505            return Sign::Zero;
506        }
507        let chain = q.square_free().sturm();
508        let mut r = self.clone();
509        loop {
510            if q.sign_at(&r.lo) != Sign::Zero && roots_in(&chain, &r.lo, &r.hi) == 0 {
511                return q.sign_at(&r.hi);
512            }
513            r.refine();
514            if r.is_exact() {
515                return q.sign_at(&r.lo);
516            }
517        }
518    }
519
520    /// Exact sign of `self - other`.
521    ///
522    /// Distinct roots are separated by refinement, which terminates because
523    /// they differ. Equality is decided exactly first: `g = gcd` of the two
524    /// polynomials. A root of `g` in the overlap of the two isolating
525    /// intervals is a root of each polynomial in each interval, and each
526    /// interval holds only one, so it is both roots. No such root means
527    /// they differ.
528    #[must_use]
529    pub fn cmp_root(&self, other: &Self) -> Sign {
530        if self.is_exact() {
531            return other.cmp_dyadic(&self.lo).flip();
532        }
533        if other.is_exact() {
534            return self.cmp_dyadic(&other.lo);
535        }
536        if let Some(sign) = disjoint(self, other) {
537            return sign;
538        }
539        let g = self.poly.gcd(&other.poly);
540        if g.degree().unwrap_or(0) >= 1 {
541            // Overlap is open at both ends, and neither end is a root of
542            // g: each end is an end of one interval, where that
543            // interval's polynomial -- a multiple of g -- is non-zero.
544            let lo = max_dyadic(&self.lo, &other.lo);
545            let hi = min_dyadic(&self.hi, &other.hi);
546            if roots_in(&g.sturm(), &lo, &hi) > 0 {
547                return Sign::Zero;
548            }
549        }
550        let (mut a, mut b) = (self.clone(), other.clone());
551        loop {
552            a.refine();
553            b.refine();
554            if a.is_exact() || b.is_exact() {
555                return a.cmp_root(&b);
556            }
557            if let Some(sign) = disjoint(&a, &b) {
558                return sign;
559            }
560        }
561    }
562}
563
564fn disjoint(a: &RealRoot, b: &RealRoot) -> Option<Sign> {
565    if a.hi.sub(&b.lo).sign() != Some(Sign::Positive) {
566        return Some(Sign::Negative);
567    }
568    if b.hi.sub(&a.lo).sign() != Some(Sign::Positive) {
569        return Some(Sign::Positive);
570    }
571    None
572}
573
574fn max_dyadic(a: &Dyadic, b: &Dyadic) -> Dyadic {
575    if a.sub(b).sign() == Some(Sign::Negative) {
576        b.clone()
577    } else {
578        a.clone()
579    }
580}
581
582fn min_dyadic(a: &Dyadic, b: &Dyadic) -> Dyadic {
583    if a.sub(b).sign() == Some(Sign::Positive) {
584        b.clone()
585    } else {
586        a.clone()
587    }
588}
589
590impl PartialOrd for RealRoot {
591    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
592        Some(match self.cmp_root(other) {
593            Sign::Negative => Ordering::Less,
594            Sign::Positive => Ordering::Greater,
595            _ => Ordering::Equal,
596        })
597    }
598}
axiolid_exact/poly.rs

axiolid_exact/
poly.rs
