tiny_skia/path64/
quad64.rs

Help
1// Copyright 2012 Google Inc.
2// Copyright 2020 Yevhenii Reizner
3//
4// Use of this source code is governed by a BSD-style license that can be
5// found in the LICENSE file.
6
7use super::Scalar64;
8
9#[cfg(all(not(feature = "std"), feature = "no-std-float"))]
10use tiny_skia_path::NoStdFloat;
11
12pub fn push_valid_ts(s: &[f64], real_roots: usize, t: &mut [f64]) -> usize {
13    let mut found_roots = 0;
14    'outer: for index in 0..real_roots {
15        let mut t_value = s[index];
16        if t_value.approximately_zero_or_more() && t_value.approximately_one_or_less() {
17            t_value = t_value.bound(0.0, 1.0);
18
19            for idx2 in 0..found_roots {
20                if t[idx2].approximately_equal(t_value) {
21                    continue 'outer;
22                }
23            }
24
25            t[found_roots] = t_value;
26            found_roots += 1;
27        }
28    }
29
30    found_roots
31}
32
33// note: caller expects multiple results to be sorted smaller first
34// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
35//  analysis of the quadratic equation, suggesting why the following looks at
36//  the sign of B -- and further suggesting that the greatest loss of precision
37//  is in b squared less two a c
38pub fn roots_valid_t(a: f64, b: f64, c: f64, t: &mut [f64]) -> usize {
39    let mut s = [0.0; 3];
40    let real_roots = roots_real(a, b, c, &mut s);
41    push_valid_ts(&s, real_roots, t)
42}
43
44// Numeric Solutions (5.6) suggests to solve the quadratic by computing
45//     Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
46// and using the roots
47//     t1 = Q / A
48//     t2 = C / Q
49//
50// this does not discard real roots <= 0 or >= 1
51pub fn roots_real(a: f64, b: f64, c: f64, s: &mut [f64; 3]) -> usize {
52    if a == 0.0 {
53        return handle_zero(b, c, s);
54    }
55
56    let p = b / (2.0 * a);
57    let q = c / a;
58    if a.approximately_zero() && (p.approximately_zero_inverse() || q.approximately_zero_inverse())
59    {
60        return handle_zero(b, c, s);
61    }
62
63    // normal form: x^2 + px + q = 0
64    let p2 = p * p;
65    if !p2.almost_dequal_ulps(q) && p2 < q {
66        return 0;
67    }
68
69    let mut sqrt_d = 0.0;
70    if p2 > q {
71        sqrt_d = (p2 - q).sqrt();
72    }
73
74    s[0] = sqrt_d - p;
75    s[1] = -sqrt_d - p;
76    1 + usize::from(!s[0].almost_dequal_ulps(s[1]))
77}
78
79fn handle_zero(b: f64, c: f64, s: &mut [f64; 3]) -> usize {
80    if b.approximately_zero() {
81        s[0] = 0.0;
82        (c == 0.0) as usize
83    } else {
84        s[0] = -c / b;
85        1
86    }
87}
tiny_skia/path64/quad64.rs

tiny_skia/path64/
quad64.rs
