pxfm/sin_cosf/argument_reduction.rs
1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1. Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2. Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3. Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::bits::get_exponent_f32;
30use crate::common::f_fmla;
31use crate::round::RoundFinite;
32
33#[derive(Debug, Copy, Clone)]
34pub(crate) struct ArgumentReducer {
35 pub(crate) x: f64,
36 pub(crate) x_abs: u32,
37}
38
39impl ArgumentReducer {
40 // Return k and y, where
41 // k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
42 #[cfg(any(
43 all(
44 any(target_arch = "x86", target_arch = "x86_64"),
45 target_feature = "fma"
46 ),
47 all(target_arch = "aarch64", target_feature = "neon")
48 ))]
49 #[inline]
50 pub(crate) fn reduce_small(self) -> (f64, i64) {
51 /*
52 Generated in SageMath:
53 z = RealField(300)(32) / RealField(300).pi()
54 n = 53
55 x_hi = RealField(n)(z) # convert to f64
56 x_mid = RealField(n)(z - RealField(300)(x_hi))
57 x_lo = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid))
58 print(double_to_hex(x_hi), ",")
59 print(double_to_hex(x_mid), ",")
60 print(double_to_hex(x_lo), ",")
61 */
62 const THIRTYTWO_OVER_PI: [u64; 3] =
63 [0x40245f306dc9c883, 0xbcc6b01ec5417056, 0xb966447e493ad4ce];
64
65 let kd = (self.x * f64::from_bits(THIRTYTWO_OVER_PI[0])).round_finite();
66 let mut y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[0]), -kd);
67 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[1]), y);
68 (y, unsafe {
69 kd.to_int_unchecked::<i64>() // indeterminate values is always filtered out before this call, as well only lowest bits are used
70 })
71 }
72
73 // Return k and y, where
74 // k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
75 #[cfg(not(any(
76 all(
77 any(target_arch = "x86", target_arch = "x86_64"),
78 target_feature = "fma"
79 ),
80 all(target_arch = "aarch64", target_feature = "neon")
81 )))]
82 #[inline]
83 pub(crate) fn reduce_small(self) -> (f64, i64) {
84 /*
85 Generated in SageMath:
86 z = RealField(300)(32) / RealField(300).pi()
87 n = 28
88 x_hi = RealField(n)(z) # convert to f64
89 x_mid = RealField(n)(z - RealField(300)(x_hi))
90 x_lo = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid))
91 print(double_to_hex(x_hi), ",")
92 print(double_to_hex(x_mid), ",")
93 print(double_to_hex(x_lo), ",")
94 */
95 const THIRTYTWO_OVER_PI: [u64; 3] =
96 [0x40245f306e000000, 0xbe3b1bbeae000000, 0x3c63f84eb0000000];
97 let prod = self.x * f64::from_bits(THIRTYTWO_OVER_PI[0]);
98 let kd = prod.round_finite();
99 let mut y = prod - kd;
100 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[1]), y);
101 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[2]), y);
102 (y, unsafe {
103 kd.to_int_unchecked::<i64>() // indeterminate values is always filtered out before this call, as well only lowest bits are used
104 })
105 }
106
107 // Return k and y, where
108 // k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
109 #[cfg(any(
110 all(
111 any(target_arch = "x86", target_arch = "x86_64"),
112 target_feature = "fma"
113 ),
114 all(target_arch = "aarch64", target_feature = "neon")
115 ))]
116 #[inline]
117 pub(crate) fn reduce_large(self, exponent: i32) -> (f64, i64) {
118 /*
119 Generated in SageMath:
120 z = RealField(300)(32) / RealField(300).pi()
121 n = 53
122 x_hi = RealField(n)(z) # convert to f64
123 x_mid = RealField(n)(z - RealField(300)(x_hi))
124 x_lo = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid))
125 x_lo_2 = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid) - RealField(300)(x_lo))
126 x_lo_3 = z - RealField(300)(x_hi) - RealField(300)(x_mid) - RealField(300)(x_lo) - RealField(300)(x_lo_2)
127 print(double_to_hex(x_hi), ",")
128 print(double_to_hex(x_mid), ",")
129 print(double_to_hex(x_lo), ",")
130 print(double_to_hex(x_lo_2), ",")
131 print(double_to_hex(x_lo_3), ",")
132 */
133 const THIRTYTWO_OVER_PI: [u64; 5] = [
134 0x40245f306dc9c883,
135 0xbcc6b01ec5417056,
136 0xb966447e493ad4ce,
137 0x360e21c820ff28b2,
138 0xb29508510ea79237,
139 ];
140
141 // 2^45 <= |x| < 2^99
142 if exponent < 99 {
143 // - When x < 2^99, the full exact product of x * THIRTYTWO_OVER_PI[0]
144 // contains at least one integral bit <= 2^5.
145 // - When 2^45 <= |x| < 2^55, the lowest 6 unit bits are contained
146 // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[0]).
147 // - When |x| >= 2^55, the LSB of double(x * THIRTYTWO_OVER_PI[0]) is at
148 // least 2^6.
149 let mut prod_hi = self.x * f64::from_bits(THIRTYTWO_OVER_PI[0]);
150 prod_hi = f64::from_bits(
151 prod_hi.to_bits()
152 & (if exponent < 55 {
153 0xfffffffffffff000
154 } else {
155 0xffffffffffffffff
156 }),
157 ); // |x| < 2^55
158 let k_hi = prod_hi.round_finite();
159 let truncated_prod = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[0]), -k_hi);
160 let prod_lo = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[1]), truncated_prod);
161 let k_lo = prod_lo.round_finite();
162 let mut y = f_fmla(
163 self.x,
164 f64::from_bits(THIRTYTWO_OVER_PI[1]),
165 truncated_prod - k_lo,
166 );
167 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[2]), y);
168 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[3]), y);
169
170 return (y, k_lo as i64);
171 }
172
173 // - When x >= 2^110, the full exact product of x * THIRTYTWO_OVER_PI[0] does
174 // not contain any of the lowest 6 unit bits, so we can ignore it completely.
175 // - When 2^99 <= |x| < 2^110, the lowest 6 unit bits are contained
176 // in the last 12 bits of double(x * THIRTYTWO_OVER_PI[1]).
177 // - When |x| >= 2^110, the LSB of double(x * THIRTYTWO_OVER_PI[1]) is at
178 // least 64.
179 let mut prod_hi = self.x * f64::from_bits(THIRTYTWO_OVER_PI[1]);
180 prod_hi = f64::from_bits(
181 prod_hi.to_bits()
182 & (if exponent < 110 {
183 0xfffffffffffff000
184 } else {
185 0xffffffffffffffff
186 }),
187 ); // |x| < 2^110
188 let k_hi = prod_hi.round_finite();
189 let truncated_prod = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[1]), -k_hi);
190 let prod_lo = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[2]), truncated_prod);
191 let k_lo = prod_lo.round_finite();
192 let mut y = f_fmla(
193 self.x,
194 f64::from_bits(THIRTYTWO_OVER_PI[2]),
195 truncated_prod - k_lo,
196 );
197 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[3]), y);
198 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[4]), y);
199
200 (y, k_lo as i64)
201 }
202
203 // Return k and y, where
204 // k = round(x * 32 / pi) and y = (x * 32 / pi) - k.
205 #[cfg(not(any(
206 all(
207 any(target_arch = "x86", target_arch = "x86_64"),
208 target_feature = "fma"
209 ),
210 all(target_arch = "aarch64", target_feature = "neon")
211 )))]
212 #[inline]
213 pub(crate) fn reduce_large(self, exponent: i32) -> (f64, i64) {
214 // For large range, there are at most 2 parts of THIRTYTWO_OVER_PI_28
215 // contributing to the lowest 6 binary digits (k & 63). If the least
216 // significant bit of x * the least significant bit of THIRTYTWO_OVER_PI_28[i]
217 // >= 64, we can completely ignore THIRTYTWO_OVER_PI_28[i].
218
219 // Generated by SageMath:
220 // z = RealField(300)(32) / RealField(300).pi()
221 // n = 28
222 // x_hi = RealField(n)(z) # convert to f64
223 // x_mid = RealField(n)(z - RealField(300)(x_hi))
224 // x_lo = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid))
225 // x_lo_2 = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid) - RealField(300)(x_lo))
226 // x_lo_3 = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid) - RealField(300)(x_lo) - RealField(300)(x_lo_2))
227 // x_lo_4 = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid) - RealField(300)(x_lo) - RealField(300)(x_lo_2) - RealField(300)(x_lo_3))
228 // x_lo_5 = RealField(n)(z - RealField(300)(x_hi) - RealField(300)(x_mid) - RealField(300)(x_lo) - RealField(300)(x_lo_2) - RealField(300)(x_lo_3) - RealField(300)(x_lo_4))
229 // x_lo_6 = (z - RealField(300)(x_hi) - RealField(300)(x_mid) - RealField(300)(x_lo) - RealField(300)(x_lo_2) - RealField(300)(x_lo_3) - RealField(300)(x_lo_4) - RealField(300)(x_lo_5))
230 //
231 //
232 // print(double_to_hex(x_hi), ",")
233 // print(double_to_hex(x_mid), ",")
234 // print(double_to_hex(x_lo), ",")
235 // print(double_to_hex(x_lo_2), ",")
236 // print(double_to_hex(x_lo_3), ",")
237 // print(double_to_hex(x_lo_4), ",")
238 // print(double_to_hex(x_lo_5), ",")
239 // print(double_to_hex(x_lo_6), ",")
240 static THIRTYTWO_OVER_PI: [u64; 8] = [
241 0x40245f306e000000,
242 0xbe3b1bbeae000000,
243 0x3c63f84eb0000000,
244 0xba87056592000000,
245 0x38bc0db62a000000,
246 0xb6e4cd8778000000,
247 0xb51bef806c000000,
248 0x33363abdebbc561b,
249 ];
250
251 let mut idx = 0;
252 const FRACTION_LEN: i32 = 23;
253 let x_lsb_exp_m4 = exponent - FRACTION_LEN;
254
255 // Exponents of the least significant bits of the corresponding entries in
256 // THIRTYTWO_OVER_PI_28.
257 static THIRTYTWO_OVER_PI_28_LSB_EXP: [i32; 8] =
258 [-24, -55, -81, -114, -143, -170, -200, -230];
259
260 // Skipping the first parts of 32/pi such that:
261 // LSB of x * LSB of THIRTYTWO_OVER_PI_28[i] >= 32.
262 while x_lsb_exp_m4 + THIRTYTWO_OVER_PI_28_LSB_EXP[idx] > 5 {
263 idx += 1;
264 }
265
266 let prod_hi = self.x * f64::from_bits(THIRTYTWO_OVER_PI[idx]);
267 // Get the integral part of x * THIRTYTWO_OVER_PI_28[idx]
268 let k_hi = prod_hi.round_finite();
269 // Get the fractional part of x * THIRTYTWO_OVER_PI_28[idx]
270 let frac = prod_hi - k_hi;
271 let prod_lo = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[idx + 1]), frac);
272 let k_lo = prod_lo.round_finite();
273
274 // Now y is the fractional parts.
275 let mut y = prod_lo - k_lo;
276 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[idx + 2]), y);
277 y = f_fmla(self.x, f64::from_bits(THIRTYTWO_OVER_PI[idx + 3]), y);
278
279 (y, (k_hi as i64).wrapping_add(k_lo as i64))
280 }
281
282 #[inline]
283 pub(crate) fn reduce(self) -> (f64, i64) {
284 #[cfg(any(
285 all(
286 any(target_arch = "x86", target_arch = "x86_64"),
287 target_feature = "fma"
288 ),
289 all(target_arch = "aarch64", target_feature = "neon")
290 ))]
291 const SMALL_PASS_BOUND: u32 = 0x5600_0000u32;
292 #[cfg(not(any(
293 all(
294 any(target_arch = "x86", target_arch = "x86_64"),
295 target_feature = "fma"
296 ),
297 all(target_arch = "aarch64", target_feature = "neon")
298 )))]
299 const SMALL_PASS_BOUND: u32 = 0x4a80_0000u32;
300 if self.x_abs < SMALL_PASS_BOUND {
301 // 2^45
302 self.reduce_small()
303 } else {
304 self.reduce_large(get_exponent_f32(f32::from_bits(self.x_abs)))
305 }
306 }
307}