brep_kernel/geometry/
poly.rs

1//! Polynomial roots shared by geometric intersections.
2const DEGENERATE: f64 = 1e-9;
3
4/// Real quadratic roots, in minus-square-root then plus-square-root order.
5/// Linear/constant degeneracy thresholds remain the caller's responsibility.
6/// A repeated root is returned twice to preserve two-branch callers.
7pub fn solve_quadratic(c0: f64, c1: f64, c2: f64) -> Vec<f64> {
8    if c2 == 0.0 {
9        return if c1 != 0.0 { vec![-c0 / c1] } else { Vec::new() };
10    }
11    let discriminant = c1 * c1 - 4.0 * c2 * c0;
12    if !(discriminant >= 0.0) {
13        return Vec::new();
14    }
15    let q = -0.5 * (c1 + c1.signum() * discriminant.sqrt());
16    if q == 0.0 {
17        return vec![-c1 / (2.0 * c2); 2];
18    }
19    if c1.is_sign_negative() {
20        vec![c0 / q, q / c2]
21    } else {
22        vec![q / c2, c0 / q]
23    }
24}
25
26/// Real roots up to degree four; decline near-multiple roots and ill-scaled leading terms.
27pub fn real_roots(coefficients: &[f64]) -> Option<Vec<f64>> {
28    let mut coefficients: Vec<f64> = coefficients.to_vec();
29    if coefficients.iter().any(|c| !c.is_finite()) {
30        return None;
31    }
32    let largest = coefficients.iter().fold(0.0_f64, |m, c| m.max(c.abs()));
33    if !(largest > 0.0) {
34        return None;
35    }
36    let leading = *coefficients.last()?;
37    if leading.abs() <= DEGENERATE * largest {
38        return None;
39    }
40    while coefficients.len() > 1 && coefficients.last() == Some(&0.0) {
41        coefficients.pop();
42    }
43    let degree = coefficients.len() - 1;
44    if degree == 0 {
45        return Some(Vec::new());
46    }
47    let evaluate = |s: f64| coefficients.iter().rev().fold(0.0, |acc, c| acc * s + c);
48    let terms = |s: f64| coefficients.iter().enumerate().fold(0.0, |acc, (i, c)| acc + (c * s.powi(i as i32)).abs());
49    if degree == 1 {
50        return Some(vec![-coefficients[0] / coefficients[1]]);
51    }
52    // Cauchy's bound on the root magnitude.
53    let bound = 1.0 + coefficients[..degree].iter().map(|c| (c / leading).abs()).fold(0.0, f64::max);
54    let derivative: Vec<f64> = coefficients.iter().enumerate().skip(1).map(|(i, c)| c * i as f64).collect();
55    let mut critical = real_roots(&derivative)?;
56    critical.retain(|c| c.abs() <= bound);
57    critical.sort_by(f64::total_cmp);
58    let mut nodes = vec![-bound];
59    nodes.extend(critical.iter().copied());
60    nodes.push(bound);
61    // A critical value at zero is a double root: decline.
62    for &c in &critical {
63        if evaluate(c).abs() <= DEGENERATE * terms(c).max(largest * DEGENERATE) {
64            return None;
65        }
66    }
67    let mut roots = Vec::new();
68    for pair in nodes.windows(2) {
69        let (mut lo, mut hi) = (pair[0], pair[1]);
70        let (mut f_lo, f_hi) = (evaluate(lo), evaluate(hi));
71        if f_lo == 0.0 {
72            roots.push(lo);
73            continue;
74        }
75        if f_lo.signum() == f_hi.signum() {
76            continue;
77        }
78        for _ in 0..200 {
79            let mid = 0.5 * (lo + hi);
80            if mid <= lo || mid >= hi {
81                break;
82            }
83            let f_mid = evaluate(mid);
84            if f_mid == 0.0 {
85                lo = mid;
86                hi = mid;
87                break;
88            }
89            if f_mid.signum() == f_lo.signum() {
90                lo = mid;
91                f_lo = f_mid;
92            } else {
93                hi = mid;
94            }
95        }
96        let mut root = 0.5 * (lo + hi);
97        for _ in 0..3 {
98            let slope = derivative.iter().rev().fold(0.0, |acc, c| acc * root + c);
99            if slope == 0.0 {
100                break;
101            }
102            let next = root - evaluate(root) / slope;
103            if !next.is_finite() || next < pair[0] || next > pair[1] {
104                break;
105            }
106            root = next;
107        }
108        roots.push(root);
109    }
110    if let Some(&last) = nodes.last() {
111        if evaluate(last) == 0.0 {
112            roots.push(last);
113        }
114    }
115    Some(roots)
116}
117
brep_kernel/geometry/poly.rs

brep_kernel/geometry/
poly.rs
