axiolid_exact/
construct.rs

1//! Exact constructions over `f64` inputs, decided without division.
2//!
3//! Each construction's result is kept symbolically -- as its inputs plus a
4//! recipe -- and every question about it is a [`SignExpr`], answered by the
5//! interval filter or, when that cannot decide, exactly. Nothing is rounded
6//! until a caller asks for an approximate value for output: the `approx_*`
7//! methods are for output only, and a decision is always asked as a sign
8//! question.
9//!
10//! Where an answer is provable from structure, it is not evaluated: the two
11//! hits of one circle are ordered by [`Branch`], because subtracting two
12//! equal intervals can never certify zero.
13
14use axiolid_core::Point2;
15use axiolid_guarantees::{Certified, Sign};
16
17use crate::arith::{sign_product, Arith};
18use crate::certify::{certify, filter, require_finite, ExactError, SignExpr};
19use crate::root::{sign_root, Root2};
20
21/// The line through two distinct points, parameterised `from + t*(to - from)`.
22#[derive(Debug, Clone, Copy, PartialEq)]
23pub struct Line {
24    from: Point2,
25    to: Point2,
26}
27
28impl Line {
29    /// The line through `from` and `to`.
30    pub fn new(from: Point2, to: Point2) -> Result<Self, ExactError> {
31        require_finite(&[from.x, from.y, to.x, to.y])?;
32        // Exact comparison on purpose: any two distinct floats define a line.
33        if from.x == to.x && from.y == to.y {
34            return Err(ExactError::DegenerateLine);
35        }
36        Ok(Self { from, to })
37    }
38
39    /// Parameter 0.
40    #[must_use]
41    pub const fn from(self) -> Point2 {
42        self.from
43    }
44
45    /// Parameter 1.
46    #[must_use]
47    pub const fn to(self) -> Point2 {
48        self.to
49    }
50}
51
52/// A circle; a zero radius is allowed and behaves as a point.
53#[derive(Debug, Clone, Copy, PartialEq)]
54pub struct Circle {
55    centre: Point2,
56    radius: f64,
57}
58
59impl Circle {
60    /// The circle about `centre` with `radius`.
61    pub fn new(centre: Point2, radius: f64) -> Result<Self, ExactError> {
62        require_finite(&[centre.x, centre.y, radius])?;
63        if radius < 0.0 {
64            return Err(ExactError::NegativeRadius);
65        }
66        Ok(Self { centre, radius })
67    }
68
69    /// Centre.
70    #[must_use]
71    pub const fn centre(self) -> Point2 {
72        self.centre
73    }
74
75    /// Radius.
76    #[must_use]
77    pub const fn radius(self) -> f64 {
78        self.radius
79    }
80}
81
82fn v<T: Arith>(value: f64) -> T {
83    T::from_f64(value)
84}
85
86/// `orient2d(a, b, p)` as a polynomial.
87fn orient<T: Arith>(a: Point2, b: Point2, px: &T, py: &T) -> T {
88    let (ax, ay) = (v::<T>(a.x), v::<T>(a.y));
89    let bax = v::<T>(b.x).sub(&ax);
90    let bay = v::<T>(b.y).sub(&ay);
91    bax.mul(&py.sub(&ay)).sub(&bay.mul(&px.sub(&ax)))
92}
93
94// ---------------------------------------------------------------- crossings
95
96/// Where lines `first` and `second` cross, `(Nx/D, Ny/D)`, homogeneously.
97fn crossing<T: Arith>(first: Line, second: Line) -> (T, T, T) {
98    let (x1, y1) = (v::<T>(first.from.x), v::<T>(first.from.y));
99    let (x2, y2) = (v::<T>(first.to.x), v::<T>(first.to.y));
100    let (x3, y3) = (v::<T>(second.from.x), v::<T>(second.from.y));
101    let (x4, y4) = (v::<T>(second.to.x), v::<T>(second.to.y));
102    let d = x1
103        .sub(&x2)
104        .mul(&y3.sub(&y4))
105        .sub(&y1.sub(&y2).mul(&x3.sub(&x4)));
106    let c1 = x1.mul(&y2).sub(&y1.mul(&x2));
107    let c2 = x3.mul(&y4).sub(&y3.mul(&x4));
108    let nx = c1.mul(&x3.sub(&x4)).sub(&x1.sub(&x2).mul(&c2));
109    let ny = c1.mul(&y3.sub(&y4)).sub(&y1.sub(&y2).mul(&c2));
110    (nx, ny, d)
111}
112
113struct CrossingDenominator(Line, Line);
114
115impl SignExpr for CrossingDenominator {
116    fn sign_in<T: Arith>(&self) -> Option<Sign> {
117        crossing::<T>(self.0, self.1).2.sign()
118    }
119}
120
121struct CrossingOrientation {
122    first: Line,
123    second: Line,
124    a: Point2,
125    b: Point2,
126}
127
128impl SignExpr for CrossingOrientation {
129    fn sign_in<T: Arith>(&self) -> Option<Sign> {
130        let (nx, ny, d) = crossing::<T>(self.first, self.second);
131        // orient(a, b, N/D) * D, which has no division: substitute
132        // p = N/D into orient and multiply through by D.
133        let (ax, ay) = (v::<T>(self.a.x), v::<T>(self.a.y));
134        let bax = v::<T>(self.b.x).sub(&ax);
135        let bay = v::<T>(self.b.y).sub(&ay);
136        let scaled = bax
137            .mul(&ny.sub(&ay.mul(&d)))
138            .sub(&bay.mul(&nx.sub(&ax.mul(&d))));
139        // orient = scaled / D, so its sign is sign(scaled) * sign(D).
140        Some(sign_product(scaled.sign()?, d.sign()?))
141    }
142}
143
144/// Which side of the directed line `a -> b` the crossing of `first` and
145/// `second` lies on: [`Sign::Positive`] for left.
146///
147/// Returns `Ok(None)` when the lines are parallel (no single crossing).
148/// The crossing point is never rounded, so a crossing that lies exactly on
149/// `a -> b` reports [`Sign::Zero`].
150pub fn crossing_orientation(
151    first: Line,
152    second: Line,
153    a: Point2,
154    b: Point2,
155) -> Result<Option<Sign>, ExactError> {
156    require_finite(&[a.x, a.y, b.x, b.y])?;
157    if certify(&CrossingDenominator(first, second))? == Sign::Zero {
158        return Ok(None);
159    }
160    certify(&CrossingOrientation {
161        first,
162        second,
163        a,
164        b,
165    })
166    .map(Some)
167}
168
169/// The interval filter alone for [`crossing_orientation`], so the
170/// escalation rate can be observed and tested. [`Certified::Uncertain`]
171/// means the exact tier would run; parallel lines also report uncertain.
172pub fn crossing_orientation_filter(first: Line, second: Line, a: Point2, b: Point2) -> Certified {
173    if !a.x.is_finite() || !a.y.is_finite() || !b.x.is_finite() || !b.y.is_finite() {
174        return filter(&NeverDecides);
175    }
176    filter(&CrossingOrientation {
177        first,
178        second,
179        a,
180        b,
181    })
182}
183
184/// An expression no arithmetic decides: the filter's "uncertain" for input
185/// the exact tier would refuse.
186struct NeverDecides;
187
188impl SignExpr for NeverDecides {
189    fn sign_in<T: Arith>(&self) -> Option<Sign> {
190        None
191    }
192}
193
194// ------------------------------------------------------- line meets circle
195
196/// Which of the two solutions: `Minus` has the smaller parameter.
197#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
198pub enum Branch {
199    /// `t = (-B - sqrt(disc)) / A`: the first hit along the line.
200    Minus,
201    /// `t = (-B + sqrt(disc)) / A`: the second hit along the line.
202    Plus,
203}
204
205/// A point where a line meets a circle, held exactly.
206#[derive(Debug, Clone, Copy, PartialEq)]
207pub struct LineHit {
208    line: Line,
209    circle: Circle,
210    branch: Branch,
211}
212
213/// How a line meets a circle, hits in increasing parameter order.
214#[derive(Debug, Clone, Copy, PartialEq)]
215pub enum HitCount {
216    /// The line misses the circle.
217    Missed,
218    /// The line touches the circle at exactly one point.
219    Tangent(LineHit),
220    /// The line crosses the circle at two points.
221    Secant(LineHit, LineHit),
222}
223
224/// `A t^2 + 2 B t + C = 0` for the line's parameter, with
225/// `disc = B^2 - A*C`, so `t = (-B +- sqrt(disc)) / A`.
226fn quadratic<T: Arith>(line: Line, circle: Circle) -> (T, T, T) {
227    let dx = v::<T>(line.to.x).sub(&v(line.from.x));
228    let dy = v::<T>(line.to.y).sub(&v(line.from.y));
229    let fx = v::<T>(line.from.x).sub(&v(circle.centre.x));
230    let fy = v::<T>(line.from.y).sub(&v(circle.centre.y));
231    let a = dx.square().add(&dy.square());
232    let b = dx.mul(&fx).add(&dy.mul(&fy));
233    let c = fx
234        .square()
235        .add(&fy.square())
236        .sub(&v::<T>(circle.radius).square());
237    let disc = b.square().sub(&a.mul(&c));
238    (a, b, disc)
239}
240
241struct Discriminant(Line, Circle);
242
243impl SignExpr for Discriminant {
244    fn sign_in<T: Arith>(&self) -> Option<Sign> {
245        quadratic::<T>(self.0, self.1).2.sign()
246    }
247}
248
249/// Where `line` meets `circle`, exactly: missed, tangent, or two hits.
250///
251/// Tangency is decided exactly, so a line that grazes the circle is never
252/// reported as two nearly equal hits, or as a miss, by rounding.
253pub fn line_circle_hits(line: Line, circle: Circle) -> Result<HitCount, ExactError> {
254    let hit = |branch| LineHit {
255        line,
256        circle,
257        branch,
258    };
259    Ok(match certify(&Discriminant(line, circle))? {
260        Sign::Negative => HitCount::Missed,
261        Sign::Zero => HitCount::Tangent(hit(Branch::Minus)),
262        _ => HitCount::Secant(hit(Branch::Minus), hit(Branch::Plus)),
263    })
264}
265
266impl LineHit {
267    /// The line this hit lies on.
268    #[must_use]
269    pub const fn line(self) -> Line {
270        self.line
271    }
272
273    /// The circle this hit lies on.
274    #[must_use]
275    pub const fn circle(self) -> Circle {
276        self.circle
277    }
278
279    /// Which solution of the quadratic this is.
280    #[must_use]
281    pub const fn branch(self) -> Branch {
282        self.branch
283    }
284
285    fn root_sign<T: Arith>(self) -> T {
286        match self.branch {
287            Branch::Minus => v::<T>(-1.0),
288            Branch::Plus => v::<T>(1.0),
289        }
290    }
291
292    /// The parameter as `(a + b*sqrt(c)) / d` in arithmetic `T`.
293    fn parameter<T: Arith>(self) -> Root2<T> {
294        let (a, b, disc) = quadratic::<T>(self.line, self.circle);
295        Root2 {
296            a: b.neg(),
297            b: self.root_sign(),
298            c: disc,
299            d: a,
300        }
301    }
302
303    /// Sign of `t - value`: whether this hit comes before (`Negative`),
304    /// at, or after the line parameter `value`. With `0.0` and `1.0` this
305    /// answers "is the hit inside the segment `from..to`".
306    pub fn cmp_param(self, value: f64) -> Result<Sign, ExactError> {
307        require_finite(&[value])?;
308        certify(&HitVersusParam { hit: self, value })
309    }
310
311    /// Which side of the directed line `a -> b` this hit lies on.
312    pub fn orientation(self, a: Point2, b: Point2) -> Result<Sign, ExactError> {
313        require_finite(&[a.x, a.y, b.x, b.y])?;
314        certify(&HitOrientation { hit: self, a, b })
315    }
316
317    /// The parameter rounded to `f64`, for output only.
318    ///
319    /// Never use this to make a decision; ask a sign question instead.
320    #[must_use]
321    pub fn approx_param(self) -> f64 {
322        let (line, circle) = (self.line, self.circle);
323        let (dx, dy) = (line.to.x - line.from.x, line.to.y - line.from.y);
324        let (fx, fy) = (line.from.x - circle.centre.x, line.from.y - circle.centre.y);
325        let a = dx * dx + dy * dy;
326        let b = dx * fx + dy * fy;
327        let c = fx * fx + fy * fy - circle.radius * circle.radius;
328        let root = (b * b - a * c).max(0.0).sqrt();
329        match self.branch {
330            Branch::Minus => (-b - root) / a,
331            Branch::Plus => (-b + root) / a,
332        }
333    }
334
335    /// The point rounded to `f64`, for output only.
336    #[must_use]
337    pub fn approx_point(self) -> Point2 {
338        let t = self.approx_param();
339        let (from, to) = (self.line.from, self.line.to);
340        Point2::new(from.x + t * (to.x - from.x), from.y + t * (to.y - from.y))
341    }
342}
343
344struct HitVersusParam {
345    hit: LineHit,
346    value: f64,
347}
348
349impl SignExpr for HitVersusParam {
350    fn sign_in<T: Arith>(&self) -> Option<Sign> {
351        let mut root = self.hit.parameter::<T>();
352        // t - value = (a - value*d + b*sqrt(c)) / d.
353        root.a = root.a.sub(&v::<T>(self.value).mul(&root.d));
354        root.sign()
355    }
356}
357
358struct HitOrientation {
359    hit: LineHit,
360    a: Point2,
361    b: Point2,
362}
363
364impl SignExpr for HitOrientation {
365    fn sign_in<T: Arith>(&self) -> Option<Sign> {
366        // orient(a, b, from + t*dir) is linear in t: k0 + k1*t, with
367        // k0 = orient(a, b, from) and k1 = cross(b - a, dir). Substituting
368        // t = (-B + s*sqrt(disc)) / A and multiplying by A > 0:
369        // (k0*A - k1*B) + s*k1*sqrt(disc).
370        let line = self.hit.line;
371        let (qa, qb, disc) = quadratic::<T>(line, self.hit.circle);
372        let k0 = orient(self.a, self.b, &v::<T>(line.from.x), &v::<T>(line.from.y));
373        let dx = v::<T>(line.to.x).sub(&v(line.from.x));
374        let dy = v::<T>(line.to.y).sub(&v(line.from.y));
375        let bax = v::<T>(self.b.x).sub(&v(self.a.x));
376        let bay = v::<T>(self.b.y).sub(&v(self.a.y));
377        let k1 = bax.mul(&dy).sub(&bay.mul(&dx));
378        let rational = k0.mul(&qa).sub(&k1.mul(&qb));
379        let irrational = self.hit.root_sign::<T>().mul(&k1);
380        sign_root(&rational, &irrational, &disc)
381    }
382}
383
384// ------------------------------------------------------ ordering along a line
385
386struct HitOrder(LineHit, LineHit);
387
388impl SignExpr for HitOrder {
389    fn sign_in<T: Arith>(&self) -> Option<Sign> {
390        self.0.parameter::<T>().cmp_sign(&self.1.parameter::<T>())
391    }
392}
393
394/// Sign of `t(first) - t(second)`: which hit comes first along their line.
395///
396/// Hits on different circles have different radicands; the comparison is
397/// still exact ([`crate::sign_two_roots`]). Both hits must lie on the same
398/// line, since parameters on different lines are not comparable.
399pub fn compare_along(first: LineHit, second: LineHit) -> Result<Sign, ExactError> {
400    if first.line != second.line {
401        return Err(ExactError::DifferentLines);
402    }
403    if first.circle == second.circle {
404        // Same quadratic: t = (-B -+ sqrt(disc)) / A with A = |to - from|^2
405        // > 0, so Minus < Plus whenever both exist (disc > 0). Deciding this
406        // structurally matters: evaluating it would subtract two equal
407        // intervals, which never certify zero, and escalate every time.
408        return Ok(match (first.branch, second.branch) {
409            (Branch::Minus, Branch::Plus) => Sign::Negative,
410            (Branch::Plus, Branch::Minus) => Sign::Positive,
411            _ => Sign::Zero,
412        });
413    }
414    certify(&HitOrder(first, second))
415}
416
417#[cfg(test)]
418mod tests {
419    use super::*;
420
421    /// The same-circle shortcut in `compare_along` must agree with evaluating
422    /// the comparison, which it exists to skip.
423    #[test]
424    fn same_circle_shortcut_agrees_with_evaluation() {
425        let mut state = 0xC0FF_EE00_1234_5678u64;
426        let mut f = || {
427            state ^= state << 13;
428            state ^= state >> 7;
429            state ^= state << 17;
430            ((state >> 11) as f64 / (1u64 << 53) as f64) * 200.0 - 100.0
431        };
432        let mut checked = 0;
433        while checked < 500 {
434            let line = Line::new(Point2::new(f(), f()), Point2::new(f(), f()));
435            let circle = Circle::new(Point2::new(f(), f()), f().abs());
436            let (Ok(line), Ok(circle)) = (line, circle) else {
437                continue;
438            };
439            if let Ok(HitCount::Secant(a, b)) = line_circle_hits(line, circle) {
440                for (x, y) in [(a, b), (b, a), (a, a), (b, b)] {
441                    assert_eq!(compare_along(x, y), certify(&HitOrder(x, y)));
442                }
443                checked += 1;
444            }
445        }
446    }
447}
axiolid_exact/construct.rs

axiolid_exact/
construct.rs
