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: usize = 0;
14 'outer: for index: usize in 0..real_roots {
15 let mut t_value: f64 = s[index];
16 if t_value.approximately_zero_or_more() && t_value.approximately_one_or_less() {
17 t_value = t_value.bound(min:0.0, max:1.0);
18
19 for idx2: usize 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: [f64; 3] = [0.0; 3];
40 let real_roots: usize = 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}
88