axiolid_exact/
root.rs

1//! Numbers of the form `(a + b*sqrt(c)) / d`.
2//!
3//! A line meets a circle at such a parameter; comparing two such hits, or
4//! asking which side of a line one lies on, is a sign question about them.
5//! Signs are decided by squaring with case analysis (the approach of CGAL's
6//! `Root_of_2`), never by evaluating a root, so in exact arithmetic the
7//! answer is exact, including an exact zero.
8//!
9//! Degree note: `sign_root` squares once (degree 2 in `a`, `b`, `c`);
10//! `sign_two_roots` squares twice. Mantissa length, and so exact-tier cost,
11//! grows with that degree. Both run the interval filter first.
12
13use axiolid_guarantees::Sign;
14
15use crate::arith::{sign_product, Arith};
16
17/// The value `(a + b*sqrt(c)) / d`, with `c >= 0` and `d != 0`.
18#[derive(Debug, Clone, PartialEq)]
19pub struct Root2<T> {
20    /// Rational part of the numerator.
21    pub a: T,
22    /// Coefficient of the root.
23    pub b: T,
24    /// Radicand; must not be negative.
25    pub c: T,
26    /// Denominator; must not be zero.
27    pub d: T,
28}
29
30impl<T: Arith> Root2<T> {
31    /// Sign of the value, or `None` when `T` cannot decide (or the
32    /// preconditions fail: negative radicand, zero denominator).
33    #[must_use]
34    pub fn sign(&self) -> Option<Sign> {
35        let denominator = self.d.sign()?;
36        if denominator == Sign::Zero {
37            return None;
38        }
39        Some(sign_product(
40            sign_root(&self.a, &self.b, &self.c)?,
41            denominator,
42        ))
43    }
44
45    /// Sign of `self - other`: which of the two values is larger.
46    ///
47    /// Works for different radicands. With `x = (a1 + b1*sqrt(c1)) / d1`
48    /// and `y` likewise, `x - y` has the sign of `d1 * d2` times
49    /// `(a1*d2 - a2*d1) + b1*d2*sqrt(c1) - b2*d1*sqrt(c2)`.
50    #[must_use]
51    pub fn cmp_sign(&self, other: &Self) -> Option<Sign> {
52        let d1 = self.d.sign()?;
53        let d2 = other.d.sign()?;
54        if d1 == Sign::Zero || d2 == Sign::Zero {
55            return None;
56        }
57        let p = self.a.mul(&other.d).sub(&other.a.mul(&self.d));
58        let q = self.b.mul(&other.d);
59        let r = other.b.mul(&self.d).neg();
60        let inner = sign_two_roots(&p, &q, &self.c, &r, &other.c)?;
61        Some(sign_product(sign_product(inner, d1), d2))
62    }
63}
64
65/// Sign of `a + b*sqrt(c)` for `c >= 0`.
66///
67/// If `a` and `b*sqrt(c)` agree in sign (or one is zero), that is the
68/// answer. If they disagree, the larger magnitude wins, and comparing
69/// magnitudes is comparing squares: the result is `sign(a) *
70/// sign(a^2 - b^2*c)`. An exact zero comes out when they cancel exactly.
71///
72/// Returns `None` when `T` cannot decide, or when `c` is negative (the
73/// value is not real).
74#[must_use]
75pub fn sign_root<T: Arith>(a: &T, b: &T, c: &T) -> Option<Sign> {
76    let radicand = c.sign()?;
77    if radicand == Sign::Negative {
78        return None;
79    }
80    let sa = a.sign()?;
81    // b*sqrt(c) has the sign of b, unless c is zero.
82    let sb = if radicand == Sign::Zero {
83        Sign::Zero
84    } else {
85        b.sign()?
86    };
87    if sb == Sign::Zero {
88        return Some(sa);
89    }
90    if sa == Sign::Zero || sa == sb {
91        return Some(sb);
92    }
93    let dominance = a.square().sub(&b.square().mul(c)).sign()?;
94    Some(sign_product(sa, dominance))
95}
96
97/// Sign of `p + q*sqrt(c) + r*sqrt(e)` for `c, e >= 0`.
98///
99/// Split as `u + v` with `u = p + q*sqrt(c)` and `v = r*sqrt(e)`. If they
100/// agree in sign (or one is zero) that is the answer. Otherwise the result
101/// is `sign(u) * sign(u^2 - v^2)`, and `u^2 - v^2 = (p^2 + q^2*c - r^2*e) +
102/// 2*p*q*sqrt(c)` is again one root, decided by [`sign_root`].
103#[must_use]
104pub fn sign_two_roots<T: Arith>(p: &T, q: &T, c: &T, r: &T, e: &T) -> Option<Sign> {
105    let su = sign_root(p, q, c)?;
106    let radicand = e.sign()?;
107    if radicand == Sign::Negative {
108        return None;
109    }
110    let sv = if radicand == Sign::Zero {
111        Sign::Zero
112    } else {
113        r.sign()?
114    };
115    if sv == Sign::Zero {
116        return Some(su);
117    }
118    if su == Sign::Zero || su == sv {
119        return Some(sv);
120    }
121    let rational = p.square().add(&q.square().mul(c)).sub(&r.square().mul(e));
122    let two = T::from_f64(2.0);
123    let irrational = two.mul(p).mul(q);
124    let dominance = sign_root(&rational, &irrational, c)?;
125    Some(sign_product(su, dominance))
126}
axiolid_exact/root.rs

axiolid_exact/
root.rs
