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 tiny_skia_path::{Scalar, SCALAR_MAX};
8
9#[cfg(all(not(feature = "std"), feature = "no-std-float"))]
10use tiny_skia_path::NoStdFloat;
11
12// Must be first, because of macro scope rules.
13#[macro_use]
14pub mod point64;
15
16pub mod cubic64;
17pub mod line_cubic_intersections;
18mod quad64;
19
20// The code below is from SkPathOpsTypes.
21
22const DBL_EPSILON_ERR: f64 = f64::EPSILON * 4.0;
23const FLT_EPSILON_HALF: f64 = (f32::EPSILON / 2.0) as f64;
24const FLT_EPSILON_CUBED: f64 = (f32::EPSILON * f32::EPSILON * f32::EPSILON) as f64;
25const FLT_EPSILON_INVERSE: f64 = 1.0 / f32::EPSILON as f64;
26
27pub trait Scalar64 {
28 fn bound(self, min: Self, max: Self) -> Self;
29 fn between(self, a: f64, b: f64) -> bool;
30 fn precisely_zero(self) -> bool;
31 fn approximately_zero_or_more(self) -> bool;
32 fn approximately_one_or_less(self) -> bool;
33 fn approximately_zero(self) -> bool;
34 fn approximately_zero_inverse(self) -> bool;
35 fn approximately_zero_cubed(self) -> bool;
36 fn approximately_zero_half(self) -> bool;
37 fn approximately_zero_when_compared_to(self, other: Self) -> bool;
38 fn approximately_equal(self, other: Self) -> bool;
39 fn approximately_equal_half(self, other: Self) -> bool;
40 fn almost_dequal_ulps(self, other: Self) -> bool;
41}
42
43impl Scalar64 for f64 {
44 // Works just like SkTPin, returning `max` for NaN/inf
45 fn bound(self, min: Self, max: Self) -> Self {
46 max.min(self).max(min)
47 }
48
49 /// Returns true if (a <= self <= b) || (a >= self >= b).
50 fn between(self, a: f64, b: f64) -> bool {
51 debug_assert!(
52 ((a <= self && self <= b) || (a >= self && self >= b))
53 == ((a - self) * (b - self) <= 0.0)
54 || (a.precisely_zero() && self.precisely_zero() && b.precisely_zero())
55 );
56
57 (a - self) * (b - self) <= 0.0
58 }
59
60 fn precisely_zero(self) -> bool {
61 self.abs() < DBL_EPSILON_ERR
62 }
63
64 fn approximately_zero_or_more(self) -> bool {
65 self > -f64::EPSILON
66 }
67
68 fn approximately_one_or_less(self) -> bool {
69 self < 1.0 + f64::EPSILON
70 }
71
72 fn approximately_zero(self) -> bool {
73 self.abs() < f64::EPSILON
74 }
75
76 fn approximately_zero_inverse(self) -> bool {
77 self.abs() > FLT_EPSILON_INVERSE
78 }
79
80 fn approximately_zero_cubed(self) -> bool {
81 self.abs() < FLT_EPSILON_CUBED
82 }
83
84 fn approximately_zero_half(self) -> bool {
85 self < FLT_EPSILON_HALF
86 }
87
88 fn approximately_zero_when_compared_to(self, other: Self) -> bool {
89 self == 0.0 || self.abs() < (other * (f32::EPSILON as f64)).abs()
90 }
91
92 // Use this for comparing Ts in the range of 0 to 1. For general numbers (larger and smaller) use
93 // AlmostEqualUlps instead.
94 fn approximately_equal(self, other: Self) -> bool {
95 (self - other).approximately_zero()
96 }
97
98 fn approximately_equal_half(self, other: Self) -> bool {
99 (self - other).approximately_zero_half()
100 }
101
102 fn almost_dequal_ulps(self, other: Self) -> bool {
103 if self.abs() < SCALAR_MAX as f64 && other.abs() < SCALAR_MAX as f64 {
104 (self as f32).almost_dequal_ulps(other as f32)
105 } else {
106 (self - other).abs() / self.abs().max(other.abs()) < (f32::EPSILON * 16.0) as f64
107 }
108 }
109}
110
111pub fn cube_root(x: f64) -> f64 {
112 if x.approximately_zero_cubed() {
113 return 0.0;
114 }
115
116 let result: f64 = halley_cbrt3d(x.abs());
117 if x < 0.0 {
118 -result
119 } else {
120 result
121 }
122}
123
124// cube root approximation using 3 iterations of Halley's method (double)
125fn halley_cbrt3d(d: f64) -> f64 {
126 let mut a: f64 = cbrt_5d(d);
127 a = cbrta_halleyd(a, r:d);
128 a = cbrta_halleyd(a, r:d);
129 cbrta_halleyd(a, r:d)
130}
131
132// cube root approximation using bit hack for 64-bit float
133// adapted from Kahan's cbrt
134fn cbrt_5d(d: f64) -> f64 {
135 let b1: u32 = 715094163;
136 let mut t: f64 = 0.0;
137 let pt: &mut [u32; 2] = bytemuck::cast_mut(&mut t);
138 let px: [u32; 2] = bytemuck::cast(d);
139 pt[1] = px[1] / 3 + b1;
140 t
141}
142
143// iterative cube root approximation using Halley's method (double)
144fn cbrta_halleyd(a: f64, r: f64) -> f64 {
145 let a3: f64 = a * a * a;
146 a * (a3 + r + r) / (a3 + a3 + r)
147}
148
149fn interp(a: f64, b: f64, t: f64) -> f64 {
150 a + (b - a) * t
151}
152