doubled/
f32.rs

1use super::*;
2
3impl Upper for f32 {
4    #[inline]
5    fn upper(self) -> Self {
6        f32::from_bits(self.to_bits() & 0x_ffff_f000)
7    }
8}
9
10impl FromMask for Doubled<f32> {
11    type Mask = u32;
12    fn from_mask(u0: Self::Mask, u1: Self::Mask) -> Self {
13        Self::new(f32::from_bits(u0), f32::from_bits(u1))
14    }
15}
16
17impl Check for f32 {
18    fn check(self) -> bool {
19        self.is_infinite() || self.is_nan()
20    }
21}
22
23#[inline]
24fn fabsfk(x: f32) -> f32 {
25    f32::from_bits(0x_7fff_ffff & x.to_bits())
26}
27
28impl core::convert::From<f32> for Doubled<f32> {
29    #[inline]
30    fn from(f: f32) -> Self {
31        Self::new(f, 0.)
32    }
33}
34
35impl core::convert::From<Doubled<f32>> for f32 {
36    #[inline]
37    fn from(f: Doubled<f32>) -> Self {
38        f.0 + f.1
39    }
40}
41
42impl core::convert::From<f64> for Doubled<f32> {
43    #[inline]
44    fn from(f: f64) -> Self {
45        Self::from_f64(f)
46    }
47}
48
49impl Doubled<f32> {
50    #[inline]
51    pub fn abs(self) -> Self {
52        if self.0 < 0. {
53            Self::new(-self.0, -self.1)
54        } else {
55            self
56        }
57    }
58
59    #[inline]
60    pub fn square(self) -> Self {
61        let xh = self.0.upper();
62        let xl = self.0 - xh;
63        let r0 = self.0 * self.0;
64        Self::new(
65            r0,
66            xh * xh - r0 + (xh + xh) * xl + xl * xl + self.0 * (self.1 + self.1),
67        )
68    }
69
70    #[inline]
71    pub fn square_as_f(self) -> f32 {
72        let xh = self.0.upper();
73        let xl = self.0 - xh;
74        xh * self.1 + xh * self.1 + xl * xl + (xh * xl + xh * xl) + xh * xh
75    }
76
77    #[inline]
78    pub fn mul_as_f(self, other: Self) -> f32 {
79        let xh = self.0.upper();
80        let xl = self.0 - xh;
81        let yh = other.0.upper();
82        let yl = other.0 - yh;
83        self.1 * yh + xh * other.1 + xl * yl + xh * yl + xl * yh + xh * yh
84    }
85
86    #[inline]
87    pub fn from_f64(f: f64) -> Self {
88        let x = f as f32;
89        Self::new(x, (f - (x as f64)) as f32)
90    }
91
92    pub fn recip(self) -> Self {
93        let t = 1. / self.0;
94        let dh = self.0.upper();
95        let dl = self.0 - dh;
96        let th = t.upper();
97        let tl = t - th;
98        Doubled::new(
99            t,
100            t * (1. - dh * th - dh * tl - dl * th - dl * tl - self.1 * t),
101        )
102    }
103}
104
105impl CheckOrder for Doubled<f32> {
106    #[inline]
107    fn check_order(self, other: Self) {
108        self.0.check_order(other.0)
109    }
110}
111
112impl CheckOrder<f32> for Doubled<f32> {
113    #[inline]
114    fn check_order(self, other: f32) {
115        self.0.check_order(other)
116    }
117}
118
119impl CheckOrder<Doubled<f32>> for f32 {
120    #[inline]
121    fn check_order(self, other: Doubled<f32>) {
122        self.check_order(other.0)
123    }
124}
125
126impl CheckOrder for f32 {
127    #[inline]
128    fn check_order(self, other: Self) {
129        debug_assert!(self.check() || other.check() || fabsfk(self) >= fabsfk(other));
130    }
131}
132
133impl core::ops::Add<Doubled<f32>> for f32 {
134    type Output = Doubled<f32>;
135    #[inline]
136    fn add(self, other: Doubled<f32>) -> Self::Output {
137        let r0 = self + other.0;
138        let v = r0 - self; // == other.0
139        Doubled::new(r0, self - (r0 - v) + (other.0 - v) + other.1) // [other.0+self, other.1]
140    }
141}
142
143impl core::ops::Mul for Doubled<f32> {
144    type Output = Self;
145    #[inline]
146    fn mul(self, other: Self) -> Self {
147        let xh = self.0.upper();
148        let xl = self.0 - xh;
149        let yh = other.0.upper();
150        let yl = other.0 - yh;
151        let r0 = self.0 * other.0;
152        Self::new(
153            r0,
154            xh * yh - r0 + xl * yh + xh * yl + xl * yl + self.0 * other.1 + self.1 * other.0,
155        )
156    }
157}
158
159impl core::ops::Mul<f32> for Doubled<f32> {
160    type Output = Self;
161    #[inline]
162    fn mul(self, other: f32) -> Self {
163        let xh = self.0.upper();
164        let xl = self.0 - xh;
165        let yh = other.upper();
166        let yl = other - yh;
167        let r0 = self.0 * other;
168        Self::new(
169            r0,
170            xh * yh - r0 + xl * yh + xh * yl + xl * yl + self.1 * other,
171        )
172    }
173}
174
175impl core::ops::Div for Doubled<f32> {
176    type Output = Self;
177    #[inline]
178    fn div(self, other: Self) -> Self {
179        let t = 1. / other.0;
180        let dh = other.0.upper();
181        let dl = other.0 - dh;
182        let th = t.upper();
183        let tl = t - th;
184        let nhh = self.0.upper();
185        let nhl = self.0 - nhh;
186
187        let q0 = self.0 * t;
188
189        let u = -q0
190            + nhh * th
191            + nhh * tl
192            + nhl * th
193            + nhl * tl
194            + q0 * (1. - dh * th - dh * tl - dl * th - dl * tl);
195
196        Self::new(q0, t * (self.1 - q0 * other.1) + u)
197    }
198}
199
200impl AsDoubled for f32 {
201    #[inline]
202    fn as_doubled(self) -> Doubled<Self> {
203        Doubled::new(self, 0.)
204    }
205}
206
207impl MulAsDoubled for f32 {
208    #[inline]
209    fn mul_as_doubled(self, other: Self) -> Doubled<Self> {
210        let xh = self.upper();
211        let xl = self - xh;
212        let yh = other.upper();
213        let yl = other - yh;
214        let r0 = self * other;
215        Doubled::new(r0, xh * yh - r0 + xl * yh + xh * yl + xl * yl)
216    }
217}
218
219impl RecipAsDoubled for f32 {
220    fn recip_as_doubled(self) -> Doubled<f32> {
221        let t = 1. / self;
222        let dh = self.upper();
223        let dl = self - dh;
224        let th = t.upper();
225        let tl = t - th;
226        Doubled::new(t, t * (1. - dh * th - dh * tl - dl * th - dl * tl))
227    }
228}