axiolid_exact/
conic.rs

1//! Conics: exact line/conic and conic/conic intersection.
2//!
3//! A conic is `A x^2 + B xy + C y^2 + D x + E y + F = 0` with exact dyadic
4//! coefficients; circles and ellipses built from `f64` data are exact.
5//!
6//! # Line meets conic
7//!
8//! Along `p(t) = from + t*d` the conic is a quadratic in `t`, so a hit is
9//! a [`Root2`] and everything [`crate::root`] decides applies: exact
10//! tangency, exact ordering along the line, across different conics.
11//!
12//! # Conic meets conic
13//!
14//! Up to four points whose coordinates are, in general, roots of a quartic
15//! that radicals cannot express usefully. They are computed the way CGAL's
16//! algebraic kernel does it:
17//!
18//! 1. Shear `u = x + k*y` for a small integer `k`, chosen so no two
19//!    intersection points share `u` (finitely many `k` are bad).
20//! 2. Eliminate `y`: the resultant `R(u)` has the points' `u` as roots.
21//! 3. The common root in `y` is rational in `u`: `y = N(u) / D(u)` from
22//!    the first subresultant, valid because the shear made `D(u) != 0` at
23//!    every real root.
24//! 4. Each point is a [`RealRoot`] `u0` plus those polynomials. Exact
25//!    coordinates as [`RealRoot`]s follow from a second resultant, and the
26//!    sign of any conic or line at the point is the sign of one polynomial
27//!    at `u0`.
28//!
29//! Tangency is exact: in sheared coordinates the intersection multiplicity
30//! at a point is the multiplicity of `u0` as a root of `R`.
31
32use axiolid_core::Point2;
33use axiolid_guarantees::Sign;
34use num_bigint::BigInt;
35
36use crate::arith::{sign_product, Arith};
37use crate::certify::{require_finite, ExactError};
38use crate::construct::{Branch, Line};
39use crate::dyadic::Dyadic;
40use crate::poly::{IntPoly, RealRoot};
41use crate::root::Root2;
42
43/// `A x^2 + B xy + C y^2 + D x + E y + F = 0`, coefficients exact.
44#[derive(Debug, Clone, PartialEq, Eq)]
45pub struct Conic {
46    coeffs: [Dyadic; 6],
47}
48
49fn d(value: f64) -> Dyadic {
50    Dyadic::from_f64(value)
51}
52
53impl Conic {
54    /// From the six coefficients `[A, B, C, D, E, F]`.
55    ///
56    /// Refuses non-finite input, and a quadratic part that is identically
57    /// zero (that is a line; use the line APIs).
58    pub fn from_coefficients(coeffs: [f64; 6]) -> Result<Self, ExactError> {
59        require_finite(&coeffs)?;
60        Self::from_dyadic(coeffs.map(d))
61    }
62
63    fn from_dyadic(coeffs: [Dyadic; 6]) -> Result<Self, ExactError> {
64        if coeffs[..3].iter().all(|c| c.sign() == Some(Sign::Zero)) {
65            return Err(ExactError::DegenerateConic);
66        }
67        Ok(Self { coeffs })
68    }
69
70    /// The circle about `centre` with `radius > 0`.
71    pub fn circle(centre: Point2, radius: f64) -> Result<Self, ExactError> {
72        require_finite(&[centre.x, centre.y, radius])?;
73        if radius <= 0.0 {
74            return Err(ExactError::NegativeRadius);
75        }
76        let (cx, cy, r) = (d(centre.x), d(centre.y), d(radius));
77        let two = d(2.0);
78        Self::from_dyadic([
79            d(1.0),
80            Dyadic::zero(),
81            d(1.0),
82            two.mul(&cx).neg(),
83            two.mul(&cy).neg(),
84            cx.square().add(&cy.square()).sub(&r.square()),
85        ])
86    }
87
88    /// The ellipse about `centre` with semi-axis `a` along `axis` and `b`
89    /// across it. `axis` need not be unit length: the implicit form is
90    /// scaled by `|axis|^2`, so no square root or division is needed.
91    pub fn ellipse(centre: Point2, axis: Point2, a: f64, b: f64) -> Result<Self, ExactError> {
92        require_finite(&[centre.x, centre.y, axis.x, axis.y, a, b])?;
93        if a <= 0.0 || b <= 0.0 {
94            return Err(ExactError::NegativeRadius);
95        }
96        if axis.x == 0.0 && axis.y == 0.0 {
97            return Err(ExactError::DegenerateLine);
98        }
99        // b^2 ((p-c).u)^2 + a^2 ((p-c).u_perp)^2 = a^2 b^2 |u|^2,
100        // u = (ux, uy), u_perp = (-uy, ux). With q = p - c:
101        //   (q.u)^2    = ux^2 qx^2 + 2 ux uy qx qy + uy^2 qy^2
102        //   (q.uperp)^2 = uy^2 qx^2 - 2 ux uy qx qy + ux^2 qy^2
103        let (ux, uy) = (d(axis.x), d(axis.y));
104        let (a2, b2) = (d(a).square(), d(b).square());
105        let two = d(2.0);
106        let qa = b2.mul(&ux.square()).add(&a2.mul(&uy.square()));
107        let qb = two.mul(&ux).mul(&uy).mul(&b2.sub(&a2));
108        let qc = b2.mul(&uy.square()).add(&a2.mul(&ux.square()));
109        let norm2 = ux.square().add(&uy.square());
110        let rhs = a2.mul(&b2).mul(&norm2);
111        Self::from_centred(qa, qb, qc, rhs, centre)
112    }
113
114    /// `qa X^2 + qb X Y + qc Y^2 = rhs` in `X = x - cx`, `Y = y - cy`.
115    fn from_centred(
116        qa: Dyadic,
117        qb: Dyadic,
118        qc: Dyadic,
119        rhs: Dyadic,
120        centre: Point2,
121    ) -> Result<Self, ExactError> {
122        let (cx, cy) = (d(centre.x), d(centre.y));
123        let two = d(2.0);
124        let dd = two.mul(&qa).mul(&cx).add(&qb.mul(&cy)).neg();
125        let ee = two.mul(&qc).mul(&cy).add(&qb.mul(&cx)).neg();
126        let ff = qa
127            .mul(&cx.square())
128            .add(&qb.mul(&cx).mul(&cy))
129            .add(&qc.mul(&cy.square()))
130            .sub(&rhs);
131        Self::from_dyadic([qa, qb, qc, dd, ee, ff])
132    }
133
134    /// The coefficients `[A, B, C, D, E, F]`.
135    #[must_use]
136    pub fn coefficients(&self) -> &[Dyadic; 6] {
137        &self.coeffs
138    }
139
140    /// Exact value at an `f64` point.
141    #[must_use]
142    pub fn eval(&self, p: Point2) -> Dyadic {
143        let (x, y) = (d(p.x), d(p.y));
144        let [a, b, c, dd, e, f] = &self.coeffs;
145        a.mul(&x.square())
146            .add(&b.mul(&x).mul(&y))
147            .add(&c.mul(&y.square()))
148            .add(&dd.mul(&x))
149            .add(&e.mul(&y))
150            .add(f)
151    }
152
153    /// Integer coefficients with the same zero set.
154    fn integer(&self) -> [BigInt; 6] {
155        let poly = IntPoly::from_dyadic(&self.coeffs);
156        let mut out: [BigInt; 6] = Default::default();
157        // `from_dyadic` scales all six by one power of two; it trims
158        // trailing zeros, which the padding restores.
159        for (slot, c) in out.iter_mut().zip(poly.coeffs()) {
160            *slot = c.clone();
161        }
162        out
163    }
164}
165
166// ------------------------------------------------------------- line x conic
167
168/// Where a line meets a conic.
169#[derive(Debug, Clone, PartialEq)]
170pub enum ConicLineHits {
171    /// The line misses the conic.
172    None,
173    /// The line is tangent: one hit of multiplicity two.
174    Tangent(ConicLineHit),
175    /// Two distinct hits, first then second along the line.
176    Secant(ConicLineHit, ConicLineHit),
177    /// Exactly one hit because the line is parallel to an asymptote (or
178    /// the axis of a parabola): the quadratic in `t` degenerates to linear.
179    Single(ConicLineHit),
180    /// The whole line lies on the conic (a degenerate conic).
181    OnConic,
182}
183
184/// One exact hit of a line with a conic: `t = (-beta +- sqrt(disc)) / alpha`
185/// on the line, or `t = -gamma / (2 beta)` in the linear case.
186#[derive(Debug, Clone, PartialEq)]
187pub struct ConicLineHit {
188    line: Line,
189    t: Root2<Dyadic>,
190}
191
192/// Quadratic `alpha t^2 + 2 beta t + gamma` of the conic along the line.
193fn along(line: Line, conic: &Conic) -> (Dyadic, Dyadic, Dyadic) {
194    let (fx, fy) = (d(line.from().x), d(line.from().y));
195    let (dx, dy) = (d(line.to().x).sub(&fx), d(line.to().y).sub(&fy));
196    let [a, b, c, dd, e, _] = &conic.coeffs;
197    let two = d(2.0);
198    let alpha = a
199        .mul(&dx.square())
200        .add(&b.mul(&dx).mul(&dy))
201        .add(&c.mul(&dy.square()));
202    // 2*beta = 2A fx dx + B (fx dy + fy dx) + 2C fy dy + D dx + E dy
203    let two_beta = two
204        .mul(a)
205        .mul(&fx)
206        .mul(&dx)
207        .add(&b.mul(&fx.mul(&dy).add(&fy.mul(&dx))))
208        .add(&two.mul(c).mul(&fy).mul(&dy))
209        .add(&dd.mul(&dx))
210        .add(&e.mul(&dy));
211    let gamma = conic.eval(line.from());
212    // Along the line the conic is alpha t^2 + two_beta t + gamma. Return
213    // (2 alpha, two_beta, 2 gamma) = (al, be, ga): the roots are then
214    // t = (-be +- sqrt(be^2 - al*ga)) / al with no halving anywhere.
215    (two.mul(&alpha), two_beta, two.mul(&gamma))
216}
217
218/// Exact intersection of a line with a conic.
219pub fn line_conic_hits(line: Line, conic: &Conic) -> Result<ConicLineHits, ExactError> {
220    let (al, be, ga) = along(line, conic);
221    let exact = |x: &Dyadic| x.sign().expect("exact");
222    if exact(&al) == Sign::Zero {
223        return Ok(match (exact(&be), exact(&ga)) {
224            (Sign::Zero, Sign::Zero) => ConicLineHits::OnConic,
225            (Sign::Zero, _) => ConicLineHits::None,
226            // Linear: be t + gamma = 0, i.e. t = -ga / (2 be).
227            _ => ConicLineHits::Single(ConicLineHit {
228                line,
229                t: Root2 {
230                    a: ga.neg(),
231                    b: Dyadic::zero(),
232                    c: Dyadic::zero(),
233                    d: d(2.0).mul(&be),
234                },
235            }),
236        });
237    }
238    let disc = be.square().sub(&al.mul(&ga));
239    let hit = |branch: Branch| ConicLineHit {
240        line,
241        t: Root2 {
242            a: be.neg(),
243            b: match branch {
244                Branch::Minus => d(-1.0),
245                Branch::Plus => d(1.0),
246            },
247            c: disc.clone(),
248            d: al.clone(),
249        },
250    };
251    Ok(match exact(&disc) {
252        Sign::Negative => ConicLineHits::None,
253        Sign::Zero => ConicLineHits::Tangent(hit(Branch::Minus)),
254        _ => {
255            // Order along the line: t = (-be -+ sqrt)/al, so the Minus
256            // branch comes first only when al > 0.
257            let (first, second) = if exact(&al) == Sign::Positive {
258                (Branch::Minus, Branch::Plus)
259            } else {
260                (Branch::Plus, Branch::Minus)
261            };
262            ConicLineHits::Secant(hit(first), hit(second))
263        }
264    })
265}
266
267impl ConicLineHit {
268    /// The line this hit lies on.
269    #[must_use]
270    pub const fn line(&self) -> Line {
271        self.line
272    }
273
274    /// Exact sign of `t - value`.
275    pub fn cmp_param(&self, value: f64) -> Result<Sign, ExactError> {
276        require_finite(&[value])?;
277        let point = Root2 {
278            a: d(value),
279            b: Dyadic::zero(),
280            c: Dyadic::zero(),
281            d: d(1.0),
282        };
283        self.t.cmp_sign(&point).ok_or(ExactError::Undefined)
284    }
285
286    /// Exact sign of `self - other` along their common line.
287    pub fn compare_along(&self, other: &Self) -> Result<Sign, ExactError> {
288        if self.line != other.line {
289            return Err(ExactError::DifferentLines);
290        }
291        self.t.cmp_sign(&other.t).ok_or(ExactError::Undefined)
292    }
293
294    /// An approximate point, for output only.
295    #[must_use]
296    pub fn approx_point(&self) -> Point2 {
297        let t = approx_root2(&self.t);
298        let (from, to) = (self.line.from(), self.line.to());
299        Point2::new(from.x + t * (to.x - from.x), from.y + t * (to.y - from.y))
300    }
301}
302
303fn approx_root2(r: &Root2<Dyadic>) -> f64 {
304    let root = r.c.to_f64().max(0.0).sqrt();
305    (r.a.to_f64() + r.b.to_f64() * root) / r.d.to_f64()
306}
307
308// ------------------------------------------------------------ conic x conic
309
310/// How two conics meet.
311#[derive(Debug, Clone, PartialEq, Eq)]
312pub enum ConicIntersection {
313    /// Finitely many points (possibly none), in increasing order of the
314    /// sheared coordinate; use [`ConicPoint::x`] to sort by `x`.
315    Points(Vec<ConicPoint>),
316    /// The conics share a component (a whole curve), so the intersection
317    /// is not a finite point set.
318    Overlapping,
319}
320
321/// One exact intersection point of two conics.
322#[derive(Debug, Clone, PartialEq, Eq)]
323pub struct ConicPoint {
324    u: RealRoot,
325    shear: i64,
326    /// y = y_num(u) / den(u), x = x_num(u) / den(u).
327    y_num: IntPoly,
328    x_num: IntPoly,
329    den: IntPoly,
330    tangent: bool,
331}
332
333// Small integer-polynomial helpers (coefficients lowest degree first).
334fn ip(c: &[i64]) -> IntPoly {
335    IntPoly::new(c.iter().map(|&v| BigInt::from(v)).collect())
336}
337
338fn pconst(c: &BigInt) -> IntPoly {
339    IntPoly::new(vec![c.clone()])
340}
341
342fn padd(a: &IntPoly, b: &IntPoly) -> IntPoly {
343    let n = a.coeffs().len().max(b.coeffs().len());
344    let get = |p: &IntPoly, i: usize| p.coeffs().get(i).cloned().unwrap_or_default();
345    IntPoly::new((0..n).map(|i| get(a, i) + get(b, i)).collect())
346}
347
348fn pneg(a: &IntPoly) -> IntPoly {
349    IntPoly::new(a.coeffs().iter().map(|c| -c).collect())
350}
351
352fn psub(a: &IntPoly, b: &IntPoly) -> IntPoly {
353    padd(a, &pneg(b))
354}
355
356fn pmul(a: &IntPoly, b: &IntPoly) -> IntPoly {
357    if a.is_zero() || b.is_zero() {
358        return IntPoly::new(vec![]);
359    }
360    let mut out = vec![BigInt::from(0); a.coeffs().len() + b.coeffs().len() - 1];
361    for (i, x) in a.coeffs().iter().enumerate() {
362        for (j, y) in b.coeffs().iter().enumerate() {
363            out[i + j] += x * y;
364        }
365    }
366    IntPoly::new(out)
367}
368
369/// The conic in sheared coordinates `u = x + k y`, as a quadratic in `y`
370/// with coefficients in `Z[u]`: `q2 y^2 + q1(u) y + q0(u)`.
371fn sheared(c: &[BigInt; 6], k: i64) -> [IntPoly; 3] {
372    let [a, b, cc, dd, e, f] = c;
373    let k = BigInt::from(k);
374    // x = u - k y:
375    // A(u - ky)^2 + B(u - ky)y + C y^2 + D(u - ky) + E y + F
376    // y^2: A k^2 - B k + C
377    // y^1: (B - 2Ak) u + (E - Dk)
378    // y^0: A u^2 + D u + F
379    let q2 = pconst(&(a * &k * &k - b * &k + cc));
380    let q1 = IntPoly::new(vec![e - dd * &k, b - BigInt::from(2) * a * &k]);
381    let q0 = IntPoly::new(vec![f.clone(), dd.clone(), a.clone()]);
382    [q2, q1, q0]
383}
384
385/// Integer determinant by Bareiss fraction-free elimination.
386fn determinant(mut m: Vec<Vec<BigInt>>) -> BigInt {
387    let n = m.len();
388    let mut sign = BigInt::from(1);
389    let mut prev = BigInt::from(1);
390    for k in 0..n {
391        if m[k][k].sign() == num_bigint::Sign::NoSign {
392            let Some(swap) = (k + 1..n).find(|&r| m[r][k].sign() != num_bigint::Sign::NoSign)
393            else {
394                return BigInt::from(0);
395            };
396            m.swap(k, swap);
397            sign = -sign;
398        }
399        for i in k + 1..n {
400            for j in k + 1..n {
401                let v = &m[i][j] * &m[k][k] - &m[i][k] * &m[k][j];
402                m[i][j] = v / &prev;
403            }
404        }
405        prev = m[k][k].clone();
406    }
407    sign * &m[n - 1][n - 1]
408}
409
410/// Resultant of two integer polynomials with formal degrees `n`, `m`.
411fn resultant(p: &IntPoly, n: usize, q: &IntPoly, m: usize) -> BigInt {
412    let size = n + m;
413    if size == 0 {
414        return BigInt::from(1);
415    }
416    let coef = |poly: &IntPoly, i: usize| poly.coeffs().get(i).cloned().unwrap_or_default();
417    let mut rows = Vec::with_capacity(size);
418    for r in 0..m {
419        let mut row = vec![BigInt::from(0); size];
420        for i in 0..=n {
421            row[r + i] = coef(p, n - i);
422        }
423        rows.push(row);
424    }
425    for r in 0..n {
426        let mut row = vec![BigInt::from(0); size];
427        for i in 0..=m {
428            row[r + i] = coef(q, m - i);
429        }
430        rows.push(row);
431    }
432    determinant(rows)
433}
434
435/// `Res_u(r(u), den(u) * Y - num(u))` as a polynomial in `Y`, by
436/// evaluation at `Y = 0..=deg r` and exact interpolation (scaled by
437/// `deg r!`, which does not move roots).
438fn eliminate(r: &IntPoly, num: &IntPoly, den: &IntPoly) -> IntPoly {
439    let n = r.degree().unwrap_or(0);
440    let m = num.degree().unwrap_or(0).max(den.degree().unwrap_or(0));
441    let values: Vec<BigInt> = (0..=n as i64)
442        .map(|y| {
443            let line = psub(&pmul(den, &ip(&[y])), num);
444            resultant(r, n, &line, m)
445        })
446        .collect();
447    // n! * M(Y) = sum_i (-1)^(n-i) C(n,i) v_i prod_{j != i} (Y - j)
448    let mut out = IntPoly::new(vec![]);
449    for (i, v) in values.iter().enumerate() {
450        let mut term = pconst(&(v * binomial(n, i)));
451        if (n - i) % 2 == 1 {
452            term = pneg(&term);
453        }
454        for j in 0..=n {
455            if j != i {
456                term = pmul(&term, &ip(&[-(j as i64), 1]));
457            }
458        }
459        out = padd(&out, &term);
460    }
461    out
462}
463
464fn binomial(n: usize, k: usize) -> BigInt {
465    let mut out = BigInt::from(1);
466    for i in 0..k {
467        out = out * BigInt::from(n - i) / BigInt::from(i + 1);
468    }
469    out
470}
471
472/// `sum c_i * p_i` for exact dyadic scalars, scaled to integers by one
473/// positive power of two (roots and signs unchanged).
474fn combine(terms: &[(Dyadic, &IntPoly)]) -> IntPoly {
475    let len = terms
476        .iter()
477        .map(|(_, p)| p.coeffs().len())
478        .max()
479        .unwrap_or(0);
480    let coeffs: Vec<Dyadic> = (0..len)
481        .map(|i| {
482            terms.iter().fold(Dyadic::zero(), |acc, (c, p)| {
483                let pi = p.coeffs().get(i).cloned().unwrap_or_default();
484                acc.add(&c.mul(&Dyadic::from_parts(pi, 0)))
485            })
486        })
487        .collect();
488    IntPoly::from_dyadic(&coeffs)
489}
490
491/// Exact intersection of two conics.
492pub fn conic_intersections(first: &Conic, second: &Conic) -> Result<ConicIntersection, ExactError> {
493    let (c1, c2) = (first.integer(), second.integer());
494    for k in 0..=MAX_SHEAR {
495        let [a2, a1, a0] = sheared(&c1, k);
496        let [b2, b1, b0] = sheared(&c2, k);
497        // Both must stay genuinely quadratic in y, or the resultant formula
498        // below is not the resultant.
499        if a2.is_zero() || b2.is_zero() {
500            continue;
501        }
502        // y-resultant of two quadratics:
503        // (a2 b0 - a0 b2)^2 - (a2 b1 - a1 b2)(a1 b0 - a0 b1)
504        let n = psub(&pmul(&a2, &b0), &pmul(&a0, &b2));
505        let den = psub(&pmul(&b2, &a1), &pmul(&a2, &b1));
506        let res = psub(
507            &pmul(&n, &n),
508            &pmul(
509                &psub(&pmul(&a2, &b1), &pmul(&a1, &b2)),
510                &psub(&pmul(&a1, &b0), &pmul(&a0, &b1)),
511            ),
512        );
513        if res.is_zero() {
514            return Ok(ConicIntersection::Overlapping);
515        }
516        if res.degree() == Some(0) {
517            return Ok(ConicIntersection::Points(Vec::new()));
518        }
519        // y = n(u) / den(u) needs den(u0) != 0 at every real root u0; a
520        // shared real root means two points share u (or a vertical
521        // tangent in these coordinates): try the next shear.
522        let sf = res.square_free();
523        let shared = sf.gcd(&den);
524        if shared.degree().unwrap_or(0) >= 1 && !shared.real_roots().is_empty() {
525            continue;
526        }
527        let doubled = res.gcd(&res.derivative());
528        // Complex roots shared with den must go before coordinates are
529        // eliminated: there num vanishes too, so den*Y - num would share
530        // that root for every Y and the eliminant would be identically
531        // zero. Real roots are unaffected (none is shared, checked above).
532        let sf = if shared.degree().unwrap_or(0) >= 1 {
533            sf.exact_div(&shared)
534        } else {
535            sf
536        };
537        let x_num = psub(&pmul(&ip(&[0, 1]), &den), &pmul(&ip(&[k]), &n));
538        let points = sf
539            .real_roots()
540            .into_iter()
541            .map(|u| {
542                let tangent =
543                    doubled.degree().unwrap_or(0) >= 1 && u.sign_of(&doubled) == Sign::Zero;
544                ConicPoint {
545                    u,
546                    shear: k,
547                    y_num: n.clone(),
548                    x_num: x_num.clone(),
549                    den: den.clone(),
550                    tangent,
551                }
552            })
553            .collect();
554        return Ok(ConicIntersection::Points(points));
555    }
556    Err(ExactError::DegenerateConic)
557}
558
559/// Shears tried before giving up. Each pair of points and each direction
560/// of asymptote rules out at most one or two, so 0..=16 always suffices
561/// for two proper conics.
562const MAX_SHEAR: i64 = 16;
563
564impl ConicPoint {
565    /// Whether the conics are tangent here (intersection multiplicity at
566    /// least two), decided exactly.
567    #[must_use]
568    pub fn is_tangent(&self) -> bool {
569        self.tangent
570    }
571
572    /// Exact sign of the conic `other` at this point: zero on it, and for
573    /// a circle or ellipse negative inside, positive outside.
574    #[must_use]
575    pub fn sign_of_conic(&self, other: &Conic) -> Sign {
576        // Multiply through by den^2 > 0: every term is a polynomial in u.
577        let [a, b, c, dd, e, f] = other.integer();
578        let (xn, yn, den) = (&self.x_num, &self.y_num, &self.den);
579        let terms = [
580            pmul(&pconst(&a), &pmul(xn, xn)),
581            pmul(&pconst(&b), &pmul(xn, yn)),
582            pmul(&pconst(&c), &pmul(yn, yn)),
583            pmul(&pconst(&dd), &pmul(xn, den)),
584            pmul(&pconst(&e), &pmul(yn, den)),
585            pmul(&pconst(&f), &pmul(den, den)),
586        ];
587        let total = terms
588            .iter()
589            .fold(IntPoly::new(vec![]), |acc, t| padd(&acc, t));
590        if total.is_zero() {
591            return Sign::Zero;
592        }
593        self.u.sign_of(&total)
594    }
595
596    /// Exact side of the directed line `a -> b`: positive to the left.
597    #[must_use]
598    pub fn side_of_line(&self, a: Point2, b: Point2) -> Sign {
599        // Times den(u0):  dx*(y_num - ay*den) - dy*(x_num - ax*den),
600        // assembled with exact dyadic coefficients so the scaling to
601        // integers is one common positive factor.
602        let (dx, dy) = (d(b.x).sub(&d(a.x)), d(b.y).sub(&d(a.y)));
603        let shift = dx.mul(&d(a.y)).sub(&dy.mul(&d(a.x)));
604        let value = combine(&[
605            (dx, &self.y_num),
606            (dy.neg(), &self.x_num),
607            (shift.neg(), &self.den),
608        ]);
609        if value.is_zero() {
610            return Sign::Zero;
611        }
612        sign_product(self.u.sign_of(&value), self.u.sign_of(&self.den))
613    }
614
615    /// The exact `x` coordinate, as a root of an integer polynomial.
616    #[must_use]
617    pub fn x(&self) -> RealRoot {
618        self.coordinate(&self.x_num)
619    }
620
621    /// The exact `y` coordinate, as a root of an integer polynomial.
622    #[must_use]
623    pub fn y(&self) -> RealRoot {
624        self.coordinate(&self.y_num)
625    }
626
627    /// The root `num(u0) / den(u0)` among the roots of the eliminant.
628    fn coordinate(&self, num: &IntPoly) -> RealRoot {
629        let eliminant = eliminate(self.u.poly(), num, &self.den);
630        let sd = self.u.sign_of(&self.den);
631        // value - y0 has the sign of (num - y0 den)(u0) * sign(den(u0)).
632        let side = |y0: &Dyadic| {
633            let shifted = IntPoly::from_dyadic(
634                &(0..num.coeffs().len().max(self.den.coeffs().len()))
635                    .map(|i| {
636                        let n = num.coeffs().get(i).cloned().unwrap_or_default();
637                        let dd = self.den.coeffs().get(i).cloned().unwrap_or_default();
638                        Dyadic::from_parts(n, 0).sub(&y0.mul(&Dyadic::from_parts(dd, 0)))
639                    })
640                    .collect::<Vec<_>>(),
641            );
642            if shifted.is_zero() {
643                return Sign::Zero;
644            }
645            sign_product(self.u.sign_of(&shifted), sd)
646        };
647        eliminant
648            .real_roots()
649            .into_iter()
650            .find(|candidate| {
651                let (lo, hi) = candidate.bounds();
652                if candidate.is_exact() {
653                    side(lo) == Sign::Zero
654                } else {
655                    side(lo) == Sign::Positive && side(hi) == Sign::Negative
656                }
657            })
658            .expect("the coordinate is a root of its eliminant")
659    }
660
661    /// An approximate point, for output only.
662    #[must_use]
663    pub fn approx(&self) -> Point2 {
664        Point2::new(self.x().approx(), self.y().approx())
665    }
666
667    /// The shear used internally (diagnostics).
668    #[must_use]
669    pub fn shear(&self) -> i64 {
670        self.shear
671    }
672}
axiolid_exact/conic.rs

axiolid_exact/
conic.rs
