1//! Floating point division routines.
2//!
3//! This module documentation gives an overview of the method used. More documentation is inline.
4//!
5//! # Relevant notation
6//!
7//! - `m_a`: the mantissa of `a`, in base 2
8//! - `p_a`: the exponent of `a`, in base 2. I.e. `a = m_a * 2^p_a`
9//! - `uqN` (e.g. `uq1`): this refers to Q notation for fixed-point numbers. UQ1.31 is an unsigned
10//! fixed-point number with 1 integral bit, and 31 decimal bits. A `uqN` variable of type `uM`
11//! will have N bits of integer and M-N bits of fraction.
12//! - `hw`: half width, i.e. for `f64` this will be a `u32`.
13//! - `x` is the best estimate of `1/m_b`
14//!
15//! # Method Overview
16//!
17//! Division routines must solve for `a / b`, which is `res = m_a*2^p_a / m_b*2^p_b`. The basic
18//! process is as follows:
19//!
20//! - Rearange the exponent and significand to simplify the operations:
21//! `res = (m_a / m_b) * 2^{p_a - p_b}`.
22//! - Check for early exits (infinity, zero, etc).
23//! - If `a` or `b` are subnormal, normalize by shifting the mantissa and adjusting the exponent.
24//! - Set the implicit bit so math is correct.
25//! - Shift mantissa significant digits (with implicit bit) fully left such that fixed-point UQ1
26//! or UQ0 numbers can be used for mantissa math. These will have greater precision than the
27//! actual mantissa, which is important for correct rounding.
28//! - Calculate the reciprocal of `m_b`, `x`.
29//! - Use the reciprocal to multiply rather than divide: `res = m_a * x_b * 2^{p_a - p_b}`.
30//! - Reapply rounding.
31//!
32//! # Reciprocal calculation
33//!
34//! Calculating the reciprocal is the most complicated part of this process. It uses the
35//! [Newton-Raphson method], which picks an initial estimation (of the reciprocal) and performs
36//! a number of iterations to increase its precision.
37//!
38//! In general, Newton's method takes the following form:
39//!
40//! ```text
41//! `x_n` is a guess or the result of a previous iteration. Increasing `n` converges to the
42//! desired result.
43//!
44//! The result approaches a zero of `f(x)` by applying a correction to the previous gues.
45//!
46//! x_{n+1} = x_n - f(x_n) / f'(x_n)
47//! ```
48//!
49//! Applying this to find the reciprocal:
50//!
51//! ```text
52//! 1 / x = b
53//!
54//! Rearrange so we can solve by finding a zero
55//! 0 = (1 / x) - b = f(x)
56//!
57//! f'(x) = -x^{-2}
58//!
59//! x_{n+1} = 2*x_n - b*x_n^2
60//! ```
61//!
62//! This is a process that can be repeated to calculate the reciprocal with enough precision to
63//! achieve a correctly rounded result for the overall division operation. The maximum required
64//! number of iterations is known since precision doubles with each iteration.
65//!
66//! # Half-width operations
67//!
68//! Calculating the reciprocal requires widening multiplication and performing arithmetic on the
69//! results, meaning that emulated integer arithmetic on `u128` (for `f64`) and `u256` (for `f128`)
70//! gets used instead of native math.
71//!
72//! To make this more efficient, all but the final operation can be computed using half-width
73//! integers. For example, rather than computing four iterations using 128-bit integers for `f64`,
74//! we can instead perform three iterations using native 64-bit integers and only one final
75//! iteration using the full 128 bits.
76//!
77//! This works because of precision doubling. Some leeway is allowed here because the fixed-point
78//! number has more bits than the final mantissa will.
79//!
80//! [Newton-Raphson method]: https://en.wikipedia.org/wiki/Newton%27s_method
81
82use super::HalfRep;
83use crate::float::Float;
84use crate::int::{CastFrom, CastInto, DInt, HInt, Int, MinInt};
85use core::mem::size_of;
86use core::ops;
87
88fn div<F: Float>(a: F, b: F) -> F
89where
90 F::Int: CastInto<i32>,
91 F::Int: From<HalfRep<F>>,
92 F::Int: From<u8>,
93 F::Int: HInt + DInt,
94 <F::Int as HInt>::D: ops::Shr<u32, Output = <F::Int as HInt>::D>,
95 F::Int: From<u32>,
96 u16: CastInto<F::Int>,
97 i32: CastInto<F::Int>,
98 u32: CastInto<F::Int>,
99 u128: CastInto<HalfRep<F>>,
100{
101 let one = F::Int::ONE;
102 let zero = F::Int::ZERO;
103 let one_hw = HalfRep::<F>::ONE;
104 let zero_hw = HalfRep::<F>::ZERO;
105 let hw = F::BITS / 2;
106 let lo_mask = F::Int::MAX >> hw;
107
108 let significand_bits = F::SIG_BITS;
109 // Saturated exponent, representing infinity
110 let exponent_sat: F::Int = F::EXP_SAT.cast();
111
112 let exponent_bias = F::EXP_BIAS;
113 let implicit_bit = F::IMPLICIT_BIT;
114 let significand_mask = F::SIG_MASK;
115 let sign_bit = F::SIGN_MASK;
116 let abs_mask = sign_bit - one;
117 let exponent_mask = F::EXP_MASK;
118 let inf_rep = exponent_mask;
119 let quiet_bit = implicit_bit >> 1;
120 let qnan_rep = exponent_mask | quiet_bit;
121 let (mut half_iterations, full_iterations) = get_iterations::<F>();
122 let recip_precision = reciprocal_precision::<F>();
123
124 if F::BITS == 128 {
125 // FIXME(tgross35): f128 seems to require one more half iteration than expected
126 half_iterations += 1;
127 }
128
129 let a_rep = a.to_bits();
130 let b_rep = b.to_bits();
131
132 // Exponent numeric representationm not accounting for bias
133 let a_exponent = (a_rep >> significand_bits) & exponent_sat;
134 let b_exponent = (b_rep >> significand_bits) & exponent_sat;
135 let quotient_sign = (a_rep ^ b_rep) & sign_bit;
136
137 let mut a_significand = a_rep & significand_mask;
138 let mut b_significand = b_rep & significand_mask;
139
140 // The exponent of our final result in its encoded form
141 let mut res_exponent: i32 =
142 i32::cast_from(a_exponent) - i32::cast_from(b_exponent) + (exponent_bias as i32);
143
144 // Detect if a or b is zero, denormal, infinity, or NaN.
145 if a_exponent.wrapping_sub(one) >= (exponent_sat - one)
146 || b_exponent.wrapping_sub(one) >= (exponent_sat - one)
147 {
148 let a_abs = a_rep & abs_mask;
149 let b_abs = b_rep & abs_mask;
150
151 // NaN / anything = qNaN
152 if a_abs > inf_rep {
153 return F::from_bits(a_rep | quiet_bit);
154 }
155
156 // anything / NaN = qNaN
157 if b_abs > inf_rep {
158 return F::from_bits(b_rep | quiet_bit);
159 }
160
161 if a_abs == inf_rep {
162 if b_abs == inf_rep {
163 // infinity / infinity = NaN
164 return F::from_bits(qnan_rep);
165 } else {
166 // infinity / anything else = +/- infinity
167 return F::from_bits(a_abs | quotient_sign);
168 }
169 }
170
171 // anything else / infinity = +/- 0
172 if b_abs == inf_rep {
173 return F::from_bits(quotient_sign);
174 }
175
176 if a_abs == zero {
177 if b_abs == zero {
178 // zero / zero = NaN
179 return F::from_bits(qnan_rep);
180 } else {
181 // zero / anything else = +/- zero
182 return F::from_bits(quotient_sign);
183 }
184 }
185
186 // anything else / zero = +/- infinity
187 if b_abs == zero {
188 return F::from_bits(inf_rep | quotient_sign);
189 }
190
191 // a is denormal. Renormalize it and set the scale to include the necessary exponent
192 // adjustment.
193 if a_abs < implicit_bit {
194 let (exponent, significand) = F::normalize(a_significand);
195 res_exponent += exponent;
196 a_significand = significand;
197 }
198
199 // b is denormal. Renormalize it and set the scale to include the necessary exponent
200 // adjustment.
201 if b_abs < implicit_bit {
202 let (exponent, significand) = F::normalize(b_significand);
203 res_exponent -= exponent;
204 b_significand = significand;
205 }
206 }
207
208 // Set the implicit significand bit. If we fell through from the
209 // denormal path it was already set by normalize( ), but setting it twice
210 // won't hurt anything.
211 a_significand |= implicit_bit;
212 b_significand |= implicit_bit;
213
214 // Transform to a fixed-point representation by shifting the significand to the high bits. We
215 // know this is in the range [1.0, 2.0] since the implicit bit is set to 1 above.
216 let b_uq1 = b_significand << (F::BITS - significand_bits - 1);
217
218 // Align the significand of b as a UQ1.(n-1) fixed-point number in the range
219 // [1.0, 2.0) and get a UQ0.n approximate reciprocal using a small minimax
220 // polynomial approximation: x0 = 3/4 + 1/sqrt(2) - b/2.
221 // The max error for this approximation is achieved at endpoints, so
222 // abs(x0(b) - 1/b) <= abs(x0(1) - 1/1) = 3/4 - 1/sqrt(2) = 0.04289...,
223 // which is about 4.5 bits.
224 // The initial approximation is between x0(1.0) = 0.9571... and x0(2.0) = 0.4571...
225 //
226 // Then, refine the reciprocal estimate using a quadratically converging
227 // Newton-Raphson iteration:
228 // x_{n+1} = x_n * (2 - x_n * b)
229 //
230 // Let b be the original divisor considered "in infinite precision" and
231 // obtained from IEEE754 representation of function argument (with the
232 // implicit bit set). Corresponds to rep_t-sized b_UQ1 represented in
233 // UQ1.(W-1).
234 //
235 // Let b_hw be an infinitely precise number obtained from the highest (HW-1)
236 // bits of divisor significand (with the implicit bit set). Corresponds to
237 // half_rep_t-sized b_UQ1_hw represented in UQ1.(HW-1) that is a **truncated**
238 // version of b_UQ1.
239 //
240 // Let e_n := x_n - 1/b_hw
241 // E_n := x_n - 1/b
242 // abs(E_n) <= abs(e_n) + (1/b_hw - 1/b)
243 // = abs(e_n) + (b - b_hw) / (b*b_hw)
244 // <= abs(e_n) + 2 * 2^-HW
245 //
246 // rep_t-sized iterations may be slower than the corresponding half-width
247 // variant depending on the handware and whether single/double/quad precision
248 // is selected.
249 //
250 // NB: Using half-width iterations increases computation errors due to
251 // rounding, so error estimations have to be computed taking the selected
252 // mode into account!
253 let mut x_uq0 = if half_iterations > 0 {
254 // Starting with (n-1) half-width iterations
255 let b_uq1_hw: HalfRep<F> = b_uq1.hi();
256
257 // C is (3/4 + 1/sqrt(2)) - 1 truncated to W0 fractional bits as UQ0.HW
258 // with W0 being either 16 or 32 and W0 <= HW.
259 // That is, C is the aforementioned 3/4 + 1/sqrt(2) constant (from which
260 // b/2 is subtracted to obtain x0) wrapped to [0, 1) range.
261 let c_hw = c_hw::<F>();
262
263 // Check that the top bit is set, i.e. value is within `[1, 2)`.
264 debug_assert!(b_uq1_hw & (one_hw << (HalfRep::<F>::BITS - 1)) > zero_hw);
265
266 // b >= 1, thus an upper bound for 3/4 + 1/sqrt(2) - b/2 is about 0.9572,
267 // so x0 fits to UQ0.HW without wrapping.
268 let mut x_uq0_hw: HalfRep<F> =
269 c_hw.wrapping_sub(b_uq1_hw /* exact b_hw/2 as UQ0.HW */);
270
271 // An e_0 error is comprised of errors due to
272 // * x0 being an inherently imprecise first approximation of 1/b_hw
273 // * C_hw being some (irrational) number **truncated** to W0 bits
274 // Please note that e_0 is calculated against the infinitely precise
275 // reciprocal of b_hw (that is, **truncated** version of b).
276 //
277 // e_0 <= 3/4 - 1/sqrt(2) + 2^-W0
278 //
279 // By construction, 1 <= b < 2
280 // f(x) = x * (2 - b*x) = 2*x - b*x^2
281 // f'(x) = 2 * (1 - b*x)
282 //
283 // On the [0, 1] interval, f(0) = 0,
284 // then it increses until f(1/b) = 1 / b, maximum on (0, 1),
285 // then it decreses to f(1) = 2 - b
286 //
287 // Let g(x) = x - f(x) = b*x^2 - x.
288 // On (0, 1/b), g(x) < 0 <=> f(x) > x
289 // On (1/b, 1], g(x) > 0 <=> f(x) < x
290 //
291 // For half-width iterations, b_hw is used instead of b.
292 for _ in 0..half_iterations {
293 // corr_UQ1_hw can be **larger** than 2 - b_hw*x by at most 1*Ulp
294 // of corr_UQ1_hw.
295 // "0.0 - (...)" is equivalent to "2.0 - (...)" in UQ1.(HW-1).
296 // On the other hand, corr_UQ1_hw should not overflow from 2.0 to 0.0 provided
297 // no overflow occurred earlier: ((rep_t)x_UQ0_hw * b_UQ1_hw >> HW) is
298 // expected to be strictly positive because b_UQ1_hw has its highest bit set
299 // and x_UQ0_hw should be rather large (it converges to 1/2 < 1/b_hw <= 1).
300 //
301 // Now, we should multiply UQ0.HW and UQ1.(HW-1) numbers, naturally
302 // obtaining an UQ1.(HW-1) number and proving its highest bit could be
303 // considered to be 0 to be able to represent it in UQ0.HW.
304 // From the above analysis of f(x), if corr_UQ1_hw would be represented
305 // without any intermediate loss of precision (that is, in twice_rep_t)
306 // x_UQ0_hw could be at most [1.]000... if b_hw is exactly 1.0 and strictly
307 // less otherwise. On the other hand, to obtain [1.]000..., one have to pass
308 // 1/b_hw == 1.0 to f(x), so this cannot occur at all without overflow (due
309 // to 1.0 being not representable as UQ0.HW).
310 // The fact corr_UQ1_hw was virtually round up (due to result of
311 // multiplication being **first** truncated, then negated - to improve
312 // error estimations) can increase x_UQ0_hw by up to 2*Ulp of x_UQ0_hw.
313 //
314 // Now, either no overflow occurred or x_UQ0_hw is 0 or 1 in its half_rep_t
315 // representation. In the latter case, x_UQ0_hw will be either 0 or 1 after
316 // any number of iterations, so just subtract 2 from the reciprocal
317 // approximation after last iteration.
318 //
319 // In infinite precision, with 0 <= eps1, eps2 <= U = 2^-HW:
320 // corr_UQ1_hw = 2 - (1/b_hw + e_n) * b_hw + 2*eps1
321 // = 1 - e_n * b_hw + 2*eps1
322 // x_UQ0_hw = (1/b_hw + e_n) * (1 - e_n*b_hw + 2*eps1) - eps2
323 // = 1/b_hw - e_n + 2*eps1/b_hw + e_n - e_n^2*b_hw + 2*e_n*eps1 - eps2
324 // = 1/b_hw + 2*eps1/b_hw - e_n^2*b_hw + 2*e_n*eps1 - eps2
325 // e_{n+1} = -e_n^2*b_hw + 2*eps1/b_hw + 2*e_n*eps1 - eps2
326 // = 2*e_n*eps1 - (e_n^2*b_hw + eps2) + 2*eps1/b_hw
327 // \------ >0 -------/ \-- >0 ---/
328 // abs(e_{n+1}) <= 2*abs(e_n)*U + max(2*e_n^2 + U, 2 * U)
329 x_uq0_hw = next_guess(x_uq0_hw, b_uq1_hw);
330 }
331
332 // For initial half-width iterations, U = 2^-HW
333 // Let abs(e_n) <= u_n * U,
334 // then abs(e_{n+1}) <= 2 * u_n * U^2 + max(2 * u_n^2 * U^2 + U, 2 * U)
335 // u_{n+1} <= 2 * u_n * U + max(2 * u_n^2 * U + 1, 2)
336 //
337 // Account for possible overflow (see above). For an overflow to occur for the
338 // first time, for "ideal" corr_UQ1_hw (that is, without intermediate
339 // truncation), the result of x_UQ0_hw * corr_UQ1_hw should be either maximum
340 // value representable in UQ0.HW or less by 1. This means that 1/b_hw have to
341 // be not below that value (see g(x) above), so it is safe to decrement just
342 // once after the final iteration. On the other hand, an effective value of
343 // divisor changes after this point (from b_hw to b), so adjust here.
344 x_uq0_hw = x_uq0_hw.wrapping_sub(one_hw);
345
346 // Error estimations for full-precision iterations are calculated just
347 // as above, but with U := 2^-W and taking extra decrementing into account.
348 // We need at least one such iteration.
349 //
350 // Simulating operations on a twice_rep_t to perform a single final full-width
351 // iteration. Using ad-hoc multiplication implementations to take advantage
352 // of particular structure of operands.
353 let blo: F::Int = b_uq1 & lo_mask;
354
355 // x_UQ0 = x_UQ0_hw * 2^HW - 1
356 // x_UQ0 * b_UQ1 = (x_UQ0_hw * 2^HW) * (b_UQ1_hw * 2^HW + blo) - b_UQ1
357 //
358 // <--- higher half ---><--- lower half --->
359 // [x_UQ0_hw * b_UQ1_hw]
360 // + [ x_UQ0_hw * blo ]
361 // - [ b_UQ1 ]
362 // = [ result ][.... discarded ...]
363 let corr_uq1: F::Int = (F::Int::from(x_uq0_hw) * F::Int::from(b_uq1_hw)
364 + ((F::Int::from(x_uq0_hw) * blo) >> hw))
365 .wrapping_sub(one)
366 .wrapping_neg(); // account for *possible* carry
367
368 let lo_corr: F::Int = corr_uq1 & lo_mask;
369 let hi_corr: F::Int = corr_uq1 >> hw;
370
371 // x_UQ0 * corr_UQ1 = (x_UQ0_hw * 2^HW) * (hi_corr * 2^HW + lo_corr) - corr_UQ1
372 let mut x_uq0: F::Int = ((F::Int::from(x_uq0_hw) * hi_corr) << 1)
373 .wrapping_add((F::Int::from(x_uq0_hw) * lo_corr) >> (hw - 1))
374 // 1 to account for the highest bit of corr_UQ1 can be 1
375 // 1 to account for possible carry
376 // Just like the case of half-width iterations but with possibility
377 // of overflowing by one extra Ulp of x_UQ0.
378 .wrapping_sub(F::Int::from(2u8));
379
380 x_uq0 -= one;
381 // ... and then traditional fixup by 2 should work
382
383 // On error estimation:
384 // abs(E_{N-1}) <= (u_{N-1} + 2 /* due to conversion e_n -> E_n */) * 2^-HW
385 // + (2^-HW + 2^-W))
386 // abs(E_{N-1}) <= (u_{N-1} + 3.01) * 2^-HW
387 //
388 // Then like for the half-width iterations:
389 // With 0 <= eps1, eps2 < 2^-W
390 // E_N = 4 * E_{N-1} * eps1 - (E_{N-1}^2 * b + 4 * eps2) + 4 * eps1 / b
391 // abs(E_N) <= 2^-W * [ 4 * abs(E_{N-1}) + max(2 * abs(E_{N-1})^2 * 2^W + 4, 8)) ]
392 // abs(E_N) <= 2^-W * [ 4 * (u_{N-1} + 3.01) * 2^-HW + max(4 + 2 * (u_{N-1} + 3.01)^2, 8) ]
393 x_uq0
394 } else {
395 // C is (3/4 + 1/sqrt(2)) - 1 truncated to 64 fractional bits as UQ0.n
396 let c: F::Int = F::Int::from(0x7504F333u32) << (F::BITS - 32);
397 let mut x_uq0: F::Int = c.wrapping_sub(b_uq1);
398
399 // E_0 <= 3/4 - 1/sqrt(2) + 2 * 2^-64
400 // x_uq0
401 for _ in 0..full_iterations {
402 x_uq0 = next_guess(x_uq0, b_uq1);
403 }
404
405 x_uq0
406 };
407
408 // Finally, account for possible overflow, as explained above.
409 x_uq0 = x_uq0.wrapping_sub(2.cast());
410
411 // Suppose 1/b - P * 2^-W < x < 1/b + P * 2^-W
412 x_uq0 -= recip_precision.cast();
413
414 // Now 1/b - (2*P) * 2^-W < x < 1/b
415 // FIXME Is x_UQ0 still >= 0.5?
416
417 let mut quotient_uq1: F::Int = x_uq0.widen_mul(a_significand << 1).hi();
418 // Now, a/b - 4*P * 2^-W < q < a/b for q=<quotient_UQ1:dummy> in UQ1.(SB+1+W).
419
420 // quotient_UQ1 is in [0.5, 2.0) as UQ1.(SB+1),
421 // adjust it to be in [1.0, 2.0) as UQ1.SB.
422 let mut residual_lo = if quotient_uq1 < (implicit_bit << 1) {
423 // Highest bit is 0, so just reinterpret quotient_UQ1 as UQ1.SB,
424 // effectively doubling its value as well as its error estimation.
425 let residual_lo = (a_significand << (significand_bits + 1))
426 .wrapping_sub(quotient_uq1.wrapping_mul(b_significand));
427 res_exponent -= 1;
428 a_significand <<= 1;
429 residual_lo
430 } else {
431 // Highest bit is 1 (the UQ1.(SB+1) value is in [1, 2)), convert it
432 // to UQ1.SB by right shifting by 1. Least significant bit is omitted.
433 quotient_uq1 >>= 1;
434 (a_significand << significand_bits).wrapping_sub(quotient_uq1.wrapping_mul(b_significand))
435 };
436
437 // drop mutability
438 let quotient = quotient_uq1;
439
440 // NB: residualLo is calculated above for the normal result case.
441 // It is re-computed on denormal path that is expected to be not so
442 // performance-sensitive.
443 //
444 // Now, q cannot be greater than a/b and can differ by at most 8*P * 2^-W + 2^-SB
445 // Each NextAfter() increments the floating point value by at least 2^-SB
446 // (more, if exponent was incremented).
447 // Different cases (<---> is of 2^-SB length, * = a/b that is shown as a midpoint):
448 // q
449 // | | * | | | | |
450 // <---> 2^t
451 // | | | | | * | |
452 // q
453 // To require at most one NextAfter(), an error should be less than 1.5 * 2^-SB.
454 // (8*P) * 2^-W + 2^-SB < 1.5 * 2^-SB
455 // (8*P) * 2^-W < 0.5 * 2^-SB
456 // P < 2^(W-4-SB)
457 // Generally, for at most R NextAfter() to be enough,
458 // P < (2*R - 1) * 2^(W-4-SB)
459 // For f32 (0+3): 10 < 32 (OK)
460 // For f32 (2+1): 32 < 74 < 32 * 3, so two NextAfter() are required
461 // For f64: 220 < 256 (OK)
462 // For f128: 4096 * 3 < 13922 < 4096 * 5 (three NextAfter() are required)
463 //
464 // If we have overflowed the exponent, return infinity
465 if res_exponent >= i32::cast_from(exponent_sat) {
466 return F::from_bits(inf_rep | quotient_sign);
467 }
468
469 // Now, quotient <= the correctly-rounded result
470 // and may need taking NextAfter() up to 3 times (see error estimates above)
471 // r = a - b * q
472 let mut abs_result = if res_exponent > 0 {
473 let mut ret = quotient & significand_mask;
474 ret |= F::Int::from(res_exponent as u32) << significand_bits;
475 residual_lo <<= 1;
476 ret
477 } else {
478 if ((significand_bits as i32) + res_exponent) < 0 {
479 return F::from_bits(quotient_sign);
480 }
481
482 let ret = quotient.wrapping_shr(u32::cast_from(res_exponent.wrapping_neg()) + 1);
483 residual_lo = a_significand
484 .wrapping_shl(significand_bits.wrapping_add(CastInto::<u32>::cast(res_exponent)))
485 .wrapping_sub(ret.wrapping_mul(b_significand) << 1);
486 ret
487 };
488
489 residual_lo += abs_result & one; // tie to even
490 // conditionally turns the below LT comparison into LTE
491 abs_result += u8::from(residual_lo > b_significand).into();
492
493 if F::BITS == 128 || (F::BITS == 32 && half_iterations > 0) {
494 // Do not round Infinity to NaN
495 abs_result +=
496 u8::from(abs_result < inf_rep && residual_lo > (2 + 1).cast() * b_significand).into();
497 }
498
499 if F::BITS == 128 {
500 abs_result +=
501 u8::from(abs_result < inf_rep && residual_lo > (4 + 1).cast() * b_significand).into();
502 }
503
504 F::from_bits(abs_result | quotient_sign)
505}
506
507/// Calculate the number of iterations required for a float type's precision.
508///
509/// This returns `(h, f)` where `h` is the number of iterations to be done using integers at half
510/// the float's bit width, and `f` is the number of iterations done using integers of the float's
511/// full width. This is further explained in the module documentation.
512///
513/// # Requirements
514///
515/// The initial estimate should have at least 8 bits of precision. If this is not true, results
516/// will be inaccurate.
517const fn get_iterations<F: Float>() -> (usize, usize) {
518 // Precision doubles with each iteration. Assume we start with 8 bits of precision.
519 let total_iterations: usize = F::BITS.ilog2() as usize - 2;
520
521 if 2 * size_of::<F>() <= size_of::<*const ()>() {
522 // If widening multiplication will be efficient (uses word-sized integers), there is no
523 // reason to use half-sized iterations.
524 (0, total_iterations)
525 } else {
526 // Otherwise, do as many iterations as possible at half width.
527 (total_iterations - 1, 1)
528 }
529}
530
531/// `u_n` for different precisions (with N-1 half-width iterations).
532///
533/// W0 is the precision of C
534/// u_0 = (3/4 - 1/sqrt(2) + 2^-W0) * 2^HW
535///
536/// Estimated with bc:
537///
538/// ```text
539/// define half1(un) { return 2.0 * (un + un^2) / 2.0^hw + 1.0; }
540/// define half2(un) { return 2.0 * un / 2.0^hw + 2.0; }
541/// define full1(un) { return 4.0 * (un + 3.01) / 2.0^hw + 2.0 * (un + 3.01)^2 + 4.0; }
542/// define full2(un) { return 4.0 * (un + 3.01) / 2.0^hw + 8.0; }
543///
544/// | f32 (0 + 3) | f32 (2 + 1) | f64 (3 + 1) | f128 (4 + 1)
545/// u_0 | < 184224974 | < 2812.1 | < 184224974 | < 791240234244348797
546/// u_1 | < 15804007 | < 242.7 | < 15804007 | < 67877681371350440
547/// u_2 | < 116308 | < 2.81 | < 116308 | < 499533100252317
548/// u_3 | < 7.31 | | < 7.31 | < 27054456580
549/// u_4 | | | | < 80.4
550/// Final (U_N) | same as u_3 | < 72 | < 218 | < 13920
551/// ````
552///
553/// Add 2 to `U_N` due to final decrement.
554const fn reciprocal_precision<F: Float>() -> u16 {
555 let (half_iterations: usize, full_iterations: usize) = get_iterations::<F>();
556
557 if full_iterations < 1 {
558 panic!("Must have at least one full iteration");
559 }
560
561 // FIXME(tgross35): calculate this programmatically
562 if F::BITS == 32 && half_iterations == 2 && full_iterations == 1 {
563 74u16
564 } else if F::BITS == 32 && half_iterations == 0 && full_iterations == 3 {
565 10
566 } else if F::BITS == 64 && half_iterations == 3 && full_iterations == 1 {
567 220
568 } else if F::BITS == 128 && half_iterations == 4 && full_iterations == 1 {
569 13922
570 } else {
571 panic!("Invalid number of iterations")
572 }
573}
574
575/// The value of `C` adjusted to half width.
576///
577/// C is (3/4 + 1/sqrt(2)) - 1 truncated to W0 fractional bits as UQ0.HW with W0 being either
578/// 16 or 32 and W0 <= HW. That is, C is the aforementioned 3/4 + 1/sqrt(2) constant (from
579/// which b/2 is subtracted to obtain x0) wrapped to [0, 1) range.
580fn c_hw<F: Float>() -> HalfRep<F>
581where
582 F::Int: DInt,
583 u128: CastInto<HalfRep<F>>,
584{
585 const C_U128: u128 = 0x7504f333f9de6108b2fb1366eaa6a542;
586 const { C_U128 >> (u128::BITS - <HalfRep<F>>::BITS) }.cast()
587}
588
589/// Perform one iteration at any width to approach `1/b`, given previous guess `x`. Returns
590/// the next `x` as a UQ0 number.
591///
592/// This is the `x_{n+1} = 2*x_n - b*x_n^2` algorithm, implemented as `x_n * (2 - b*x_n)`. It
593/// uses widening multiplication to calculate the result with necessary precision.
594fn next_guess<I>(x_uq0: I, b_uq1: I) -> I
595where
596 I: Int + HInt,
597 <I as HInt>::D: ops::Shr<u32, Output = <I as HInt>::D>,
598{
599 // `corr = 2 - b*x_n`
600 //
601 // This looks like `0 - b*x_n`. However, this works - in `UQ1`, `0.0 - x = 2.0 - x`.
602 let corr_uq1: I = I::ZERO.wrapping_sub(x_uq0.widen_mul(b_uq1).hi());
603
604 // `x_n * corr = x_n * (2 - b*x_n)`
605 (x_uq0.widen_mul(corr_uq1) >> (I::BITS - 1)).lo()
606}
607
608intrinsics! {
609 #[avr_skip]
610 #[arm_aeabi_alias = __aeabi_fdiv]
611 pub extern "C" fn __divsf3(a: f32, b: f32) -> f32 {
612 div(a, b)
613 }
614
615 #[avr_skip]
616 #[arm_aeabi_alias = __aeabi_ddiv]
617 pub extern "C" fn __divdf3(a: f64, b: f64) -> f64 {
618 div(a, b)
619 }
620
621 #[avr_skip]
622 #[ppc_alias = __divkf3]
623 #[cfg(f128_enabled)]
624 pub extern "C" fn __divtf3(a: f128, b: f128) -> f128 {
625 div(a, b)
626 }
627
628 #[cfg(target_arch = "arm")]
629 pub extern "C" fn __divsf3vfp(a: f32, b: f32) -> f32 {
630 a / b
631 }
632
633 #[cfg(target_arch = "arm")]
634 pub extern "C" fn __divdf3vfp(a: f64, b: f64) -> f64 {
635 a / b
636 }
637}
638