1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151

// Copyright 2012 Google Inc.
// Copyright 2020 Yevhenii Reizner
//
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

use tiny_skia_path::{Scalar, SCALAR_MAX};

#[cfg(all(not(feature = "std"), feature = "no-std-float"))]
use tiny_skia_path::NoStdFloat;

// Must be first, because of macro scope rules.
#[macro_use]
pub mod point64;

pub mod cubic64;
pub mod line_cubic_intersections;
mod quad64;

// The code below is from SkPathOpsTypes.

const DBL_EPSILON_ERR: f64 = f64::EPSILON * 4.0;
const FLT_EPSILON_HALF: f64 = (f32::EPSILON / 2.0) as f64;
const FLT_EPSILON_CUBED: f64 = (f32::EPSILON * f32::EPSILON * f32::EPSILON) as f64;
const FLT_EPSILON_INVERSE: f64 = 1.0 / f32::EPSILON as f64;

pub trait Scalar64 {
    fn bound(self, min: Self, max: Self) -> Self;
    fn between(self, a: f64, b: f64) -> bool;
    fn precisely_zero(self) -> bool;
    fn approximately_zero_or_more(self) -> bool;
    fn approximately_one_or_less(self) -> bool;
    fn approximately_zero(self) -> bool;
    fn approximately_zero_inverse(self) -> bool;
    fn approximately_zero_cubed(self) -> bool;
    fn approximately_zero_half(self) -> bool;
    fn approximately_zero_when_compared_to(self, other: Self) -> bool;
    fn approximately_equal(self, other: Self) -> bool;
    fn approximately_equal_half(self, other: Self) -> bool;
    fn almost_dequal_ulps(self, other: Self) -> bool;
}

impl Scalar64 for f64 {
    // Works just like SkTPin, returning `max` for NaN/inf
    fn bound(self, min: Self, max: Self) -> Self {
        max.min(self).max(min)
    }

    /// Returns true if (a <= self <= b) || (a >= self >= b).
    fn between(self, a: f64, b: f64) -> bool {
        debug_assert!(
            ((a <= self && self <= b) || (a >= self && self >= b))
                == ((a - self) * (b - self) <= 0.0)
                || (a.precisely_zero() && self.precisely_zero() && b.precisely_zero())
        );

        (a - self) * (b - self) <= 0.0
    }

    fn precisely_zero(self) -> bool {
        self.abs() < DBL_EPSILON_ERR
    }

    fn approximately_zero_or_more(self) -> bool {
        self > -f64::EPSILON
    }

    fn approximately_one_or_less(self) -> bool {
        self < 1.0 + f64::EPSILON
    }

    fn approximately_zero(self) -> bool {
        self.abs() < f64::EPSILON
    }

    fn approximately_zero_inverse(self) -> bool {
        self.abs() > FLT_EPSILON_INVERSE
    }

    fn approximately_zero_cubed(self) -> bool {
        self.abs() < FLT_EPSILON_CUBED
    }

    fn approximately_zero_half(self) -> bool {
        self < FLT_EPSILON_HALF
    }

    fn approximately_zero_when_compared_to(self, other: Self) -> bool {
        self == 0.0 || self.abs() < (other * (f32::EPSILON as f64)).abs()
    }

    // Use this for comparing Ts in the range of 0 to 1. For general numbers (larger and smaller) use
    // AlmostEqualUlps instead.
    fn approximately_equal(self, other: Self) -> bool {
        (self - other).approximately_zero()
    }

    fn approximately_equal_half(self, other: Self) -> bool {
        (self - other).approximately_zero_half()
    }

    fn almost_dequal_ulps(self, other: Self) -> bool {
        if self.abs() < SCALAR_MAX as f64 && other.abs() < SCALAR_MAX as f64 {
            (self as f32).almost_dequal_ulps(other as f32)
        } else {
            (self - other).abs() / self.abs().max(other.abs()) < (f32::EPSILON * 16.0) as f64
        }
    }
}

pub fn cube_root(x: f64) -> f64 {
    if x.approximately_zero_cubed() {
        return 0.0;
    }

    let result = halley_cbrt3d(x.abs());
    if x < 0.0 {
        -result
    } else {
        result
    }
}

// cube root approximation using 3 iterations of Halley's method (double)
fn halley_cbrt3d(d: f64) -> f64 {
    let mut a = cbrt_5d(d);
    a = cbrta_halleyd(a, d);
    a = cbrta_halleyd(a, d);
    cbrta_halleyd(a, d)
}

// cube root approximation using bit hack for 64-bit float
// adapted from Kahan's cbrt
fn cbrt_5d(d: f64) -> f64 {
    let b1 = 715094163;
    let mut t: f64 = 0.0;
    let pt: &mut [u32; 2] = bytemuck::cast_mut(&mut t);
    let px: [u32; 2] = bytemuck::cast(d);
    pt[1] = px[1] / 3 + b1;
    t
}

// iterative cube root approximation using Halley's method (double)
fn cbrta_halleyd(a: f64, r: f64) -> f64 {
    let a3 = a * a * a;
    a * (a3 + r + r) / (a3 + a3 + r)
}

fn interp(a: f64, b: f64, t: f64) -> f64 {
    a + (b - a) * t
}