1 // The functions are complex with many branches, and explicit
2 // `return`s makes it clear where function exit points are
3 #![allow(clippy::needless_return)]
4
5 use float::Float;
6 use int::{CastInto, DInt, HInt, Int};
7
div32<F: Float>(a: F, b: F) -> F where u32: CastInto<F::Int>, F::Int: CastInto<u32>, i32: CastInto<F::Int>, F::Int: CastInto<i32>, F::Int: HInt,8 fn div32<F: Float>(a: F, b: F) -> F
9 where
10 u32: CastInto<F::Int>,
11 F::Int: CastInto<u32>,
12 i32: CastInto<F::Int>,
13 F::Int: CastInto<i32>,
14 F::Int: HInt,
15 {
16 let one = F::Int::ONE;
17 let zero = F::Int::ZERO;
18
19 // let bits = F::BITS;
20 let significand_bits = F::SIGNIFICAND_BITS;
21 let max_exponent = F::EXPONENT_MAX;
22
23 let exponent_bias = F::EXPONENT_BIAS;
24
25 let implicit_bit = F::IMPLICIT_BIT;
26 let significand_mask = F::SIGNIFICAND_MASK;
27 let sign_bit = F::SIGN_MASK as F::Int;
28 let abs_mask = sign_bit - one;
29 let exponent_mask = F::EXPONENT_MASK;
30 let inf_rep = exponent_mask;
31 let quiet_bit = implicit_bit >> 1;
32 let qnan_rep = exponent_mask | quiet_bit;
33
34 #[inline(always)]
35 fn negate_u32(a: u32) -> u32 {
36 (<i32>::wrapping_neg(a as i32)) as u32
37 }
38
39 let a_rep = a.repr();
40 let b_rep = b.repr();
41
42 let a_exponent = (a_rep >> significand_bits) & max_exponent.cast();
43 let b_exponent = (b_rep >> significand_bits) & max_exponent.cast();
44 let quotient_sign = (a_rep ^ b_rep) & sign_bit;
45
46 let mut a_significand = a_rep & significand_mask;
47 let mut b_significand = b_rep & significand_mask;
48 let mut scale = 0;
49
50 // Detect if a or b is zero, denormal, infinity, or NaN.
51 if a_exponent.wrapping_sub(one) >= (max_exponent - 1).cast()
52 || b_exponent.wrapping_sub(one) >= (max_exponent - 1).cast()
53 {
54 let a_abs = a_rep & abs_mask;
55 let b_abs = b_rep & abs_mask;
56
57 // NaN / anything = qNaN
58 if a_abs > inf_rep {
59 return F::from_repr(a_rep | quiet_bit);
60 }
61 // anything / NaN = qNaN
62 if b_abs > inf_rep {
63 return F::from_repr(b_rep | quiet_bit);
64 }
65
66 if a_abs == inf_rep {
67 if b_abs == inf_rep {
68 // infinity / infinity = NaN
69 return F::from_repr(qnan_rep);
70 } else {
71 // infinity / anything else = +/- infinity
72 return F::from_repr(a_abs | quotient_sign);
73 }
74 }
75
76 // anything else / infinity = +/- 0
77 if b_abs == inf_rep {
78 return F::from_repr(quotient_sign);
79 }
80
81 if a_abs == zero {
82 if b_abs == zero {
83 // zero / zero = NaN
84 return F::from_repr(qnan_rep);
85 } else {
86 // zero / anything else = +/- zero
87 return F::from_repr(quotient_sign);
88 }
89 }
90
91 // anything else / zero = +/- infinity
92 if b_abs == zero {
93 return F::from_repr(inf_rep | quotient_sign);
94 }
95
96 // one or both of a or b is denormal, the other (if applicable) is a
97 // normal number. Renormalize one or both of a and b, and set scale to
98 // include the necessary exponent adjustment.
99 if a_abs < implicit_bit {
100 let (exponent, significand) = F::normalize(a_significand);
101 scale += exponent;
102 a_significand = significand;
103 }
104
105 if b_abs < implicit_bit {
106 let (exponent, significand) = F::normalize(b_significand);
107 scale -= exponent;
108 b_significand = significand;
109 }
110 }
111
112 // Or in the implicit significand bit. (If we fell through from the
113 // denormal path it was already set by normalize( ), but setting it twice
114 // won't hurt anything.)
115 a_significand |= implicit_bit;
116 b_significand |= implicit_bit;
117 let mut quotient_exponent: i32 = CastInto::<i32>::cast(a_exponent)
118 .wrapping_sub(CastInto::<i32>::cast(b_exponent))
119 .wrapping_add(scale);
120
121 // Align the significand of b as a Q31 fixed-point number in the range
122 // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax
123 // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
124 // is accurate to about 3.5 binary digits.
125 let q31b = CastInto::<u32>::cast(b_significand << 8.cast());
126 let mut reciprocal = (0x7504f333u32).wrapping_sub(q31b);
127
128 // Now refine the reciprocal estimate using a Newton-Raphson iteration:
129 //
130 // x1 = x0 * (2 - x0 * b)
131 //
132 // This doubles the number of correct binary digits in the approximation
133 // with each iteration, so after three iterations, we have about 28 binary
134 // digits of accuracy.
135 let mut correction: u32;
136 correction = negate_u32(((reciprocal as u64).wrapping_mul(q31b as u64) >> 32) as u32);
137 reciprocal = ((reciprocal as u64).wrapping_mul(correction as u64) as u64 >> 31) as u32;
138 correction = negate_u32(((reciprocal as u64).wrapping_mul(q31b as u64) >> 32) as u32);
139 reciprocal = ((reciprocal as u64).wrapping_mul(correction as u64) as u64 >> 31) as u32;
140 correction = negate_u32(((reciprocal as u64).wrapping_mul(q31b as u64) >> 32) as u32);
141 reciprocal = ((reciprocal as u64).wrapping_mul(correction as u64) as u64 >> 31) as u32;
142
143 // Exhaustive testing shows that the error in reciprocal after three steps
144 // is in the interval [-0x1.f58108p-31, 0x1.d0e48cp-29], in line with our
145 // expectations. We bump the reciprocal by a tiny value to force the error
146 // to be strictly positive (in the range [0x1.4fdfp-37,0x1.287246p-29], to
147 // be specific). This also causes 1/1 to give a sensible approximation
148 // instead of zero (due to overflow).
149 reciprocal = reciprocal.wrapping_sub(2);
150
151 // The numerical reciprocal is accurate to within 2^-28, lies in the
152 // interval [0x1.000000eep-1, 0x1.fffffffcp-1], and is strictly smaller
153 // than the true reciprocal of b. Multiplying a by this reciprocal thus
154 // gives a numerical q = a/b in Q24 with the following properties:
155 //
156 // 1. q < a/b
157 // 2. q is in the interval [0x1.000000eep-1, 0x1.fffffffcp0)
158 // 3. the error in q is at most 2^-24 + 2^-27 -- the 2^24 term comes
159 // from the fact that we truncate the product, and the 2^27 term
160 // is the error in the reciprocal of b scaled by the maximum
161 // possible value of a. As a consequence of this error bound,
162 // either q or nextafter(q) is the correctly rounded
163 let mut quotient = (a_significand << 1).widen_mul(reciprocal.cast()).hi();
164
165 // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
166 // In either case, we are going to compute a residual of the form
167 //
168 // r = a - q*b
169 //
170 // We know from the construction of q that r satisfies:
171 //
172 // 0 <= r < ulp(q)*b
173 //
174 // if r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
175 // already have the correct result. The exact halfway case cannot occur.
176 // We also take this time to right shift quotient if it falls in the [1,2)
177 // range and adjust the exponent accordingly.
178 let residual = if quotient < (implicit_bit << 1) {
179 quotient_exponent = quotient_exponent.wrapping_sub(1);
180 (a_significand << (significand_bits + 1)).wrapping_sub(quotient.wrapping_mul(b_significand))
181 } else {
182 quotient >>= 1;
183 (a_significand << significand_bits).wrapping_sub(quotient.wrapping_mul(b_significand))
184 };
185
186 let written_exponent = quotient_exponent.wrapping_add(exponent_bias as i32);
187
188 if written_exponent >= max_exponent as i32 {
189 // If we have overflowed the exponent, return infinity.
190 return F::from_repr(inf_rep | quotient_sign);
191 } else if written_exponent < 1 {
192 // Flush denormals to zero. In the future, it would be nice to add
193 // code to round them correctly.
194 return F::from_repr(quotient_sign);
195 } else {
196 let round = ((residual << 1) > b_significand) as u32;
197 // Clear the implicit bits
198 let mut abs_result = quotient & significand_mask;
199 // Insert the exponent
200 abs_result |= written_exponent.cast() << significand_bits;
201 // Round
202 abs_result = abs_result.wrapping_add(round.cast());
203 // Insert the sign and return
204 return F::from_repr(abs_result | quotient_sign);
205 }
206 }
207
div64<F: Float>(a: F, b: F) -> F where u32: CastInto<F::Int>, F::Int: CastInto<u32>, i32: CastInto<F::Int>, F::Int: CastInto<i32>, u64: CastInto<F::Int>, F::Int: CastInto<u64>, i64: CastInto<F::Int>, F::Int: CastInto<i64>, F::Int: HInt,208 fn div64<F: Float>(a: F, b: F) -> F
209 where
210 u32: CastInto<F::Int>,
211 F::Int: CastInto<u32>,
212 i32: CastInto<F::Int>,
213 F::Int: CastInto<i32>,
214 u64: CastInto<F::Int>,
215 F::Int: CastInto<u64>,
216 i64: CastInto<F::Int>,
217 F::Int: CastInto<i64>,
218 F::Int: HInt,
219 {
220 let one = F::Int::ONE;
221 let zero = F::Int::ZERO;
222
223 // let bits = F::BITS;
224 let significand_bits = F::SIGNIFICAND_BITS;
225 let max_exponent = F::EXPONENT_MAX;
226
227 let exponent_bias = F::EXPONENT_BIAS;
228
229 let implicit_bit = F::IMPLICIT_BIT;
230 let significand_mask = F::SIGNIFICAND_MASK;
231 let sign_bit = F::SIGN_MASK as F::Int;
232 let abs_mask = sign_bit - one;
233 let exponent_mask = F::EXPONENT_MASK;
234 let inf_rep = exponent_mask;
235 let quiet_bit = implicit_bit >> 1;
236 let qnan_rep = exponent_mask | quiet_bit;
237 // let exponent_bits = F::EXPONENT_BITS;
238
239 #[inline(always)]
240 fn negate_u32(a: u32) -> u32 {
241 (<i32>::wrapping_neg(a as i32)) as u32
242 }
243
244 #[inline(always)]
245 fn negate_u64(a: u64) -> u64 {
246 (<i64>::wrapping_neg(a as i64)) as u64
247 }
248
249 let a_rep = a.repr();
250 let b_rep = b.repr();
251
252 let a_exponent = (a_rep >> significand_bits) & max_exponent.cast();
253 let b_exponent = (b_rep >> significand_bits) & max_exponent.cast();
254 let quotient_sign = (a_rep ^ b_rep) & sign_bit;
255
256 let mut a_significand = a_rep & significand_mask;
257 let mut b_significand = b_rep & significand_mask;
258 let mut scale = 0;
259
260 // Detect if a or b is zero, denormal, infinity, or NaN.
261 if a_exponent.wrapping_sub(one) >= (max_exponent - 1).cast()
262 || b_exponent.wrapping_sub(one) >= (max_exponent - 1).cast()
263 {
264 let a_abs = a_rep & abs_mask;
265 let b_abs = b_rep & abs_mask;
266
267 // NaN / anything = qNaN
268 if a_abs > inf_rep {
269 return F::from_repr(a_rep | quiet_bit);
270 }
271 // anything / NaN = qNaN
272 if b_abs > inf_rep {
273 return F::from_repr(b_rep | quiet_bit);
274 }
275
276 if a_abs == inf_rep {
277 if b_abs == inf_rep {
278 // infinity / infinity = NaN
279 return F::from_repr(qnan_rep);
280 } else {
281 // infinity / anything else = +/- infinity
282 return F::from_repr(a_abs | quotient_sign);
283 }
284 }
285
286 // anything else / infinity = +/- 0
287 if b_abs == inf_rep {
288 return F::from_repr(quotient_sign);
289 }
290
291 if a_abs == zero {
292 if b_abs == zero {
293 // zero / zero = NaN
294 return F::from_repr(qnan_rep);
295 } else {
296 // zero / anything else = +/- zero
297 return F::from_repr(quotient_sign);
298 }
299 }
300
301 // anything else / zero = +/- infinity
302 if b_abs == zero {
303 return F::from_repr(inf_rep | quotient_sign);
304 }
305
306 // one or both of a or b is denormal, the other (if applicable) is a
307 // normal number. Renormalize one or both of a and b, and set scale to
308 // include the necessary exponent adjustment.
309 if a_abs < implicit_bit {
310 let (exponent, significand) = F::normalize(a_significand);
311 scale += exponent;
312 a_significand = significand;
313 }
314
315 if b_abs < implicit_bit {
316 let (exponent, significand) = F::normalize(b_significand);
317 scale -= exponent;
318 b_significand = significand;
319 }
320 }
321
322 // Or in the implicit significand bit. (If we fell through from the
323 // denormal path it was already set by normalize( ), but setting it twice
324 // won't hurt anything.)
325 a_significand |= implicit_bit;
326 b_significand |= implicit_bit;
327 let mut quotient_exponent: i32 = CastInto::<i32>::cast(a_exponent)
328 .wrapping_sub(CastInto::<i32>::cast(b_exponent))
329 .wrapping_add(scale);
330
331 // Align the significand of b as a Q31 fixed-point number in the range
332 // [1, 2.0) and get a Q32 approximate reciprocal using a small minimax
333 // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2. This
334 // is accurate to about 3.5 binary digits.
335 let q31b = CastInto::<u32>::cast(b_significand >> 21.cast());
336 let mut recip32 = (0x7504f333u32).wrapping_sub(q31b);
337
338 // Now refine the reciprocal estimate using a Newton-Raphson iteration:
339 //
340 // x1 = x0 * (2 - x0 * b)
341 //
342 // This doubles the number of correct binary digits in the approximation
343 // with each iteration, so after three iterations, we have about 28 binary
344 // digits of accuracy.
345 let mut correction32: u32;
346 correction32 = negate_u32(((recip32 as u64).wrapping_mul(q31b as u64) >> 32) as u32);
347 recip32 = ((recip32 as u64).wrapping_mul(correction32 as u64) >> 31) as u32;
348 correction32 = negate_u32(((recip32 as u64).wrapping_mul(q31b as u64) >> 32) as u32);
349 recip32 = ((recip32 as u64).wrapping_mul(correction32 as u64) >> 31) as u32;
350 correction32 = negate_u32(((recip32 as u64).wrapping_mul(q31b as u64) >> 32) as u32);
351 recip32 = ((recip32 as u64).wrapping_mul(correction32 as u64) >> 31) as u32;
352
353 // recip32 might have overflowed to exactly zero in the preceeding
354 // computation if the high word of b is exactly 1.0. This would sabotage
355 // the full-width final stage of the computation that follows, so we adjust
356 // recip32 downward by one bit.
357 recip32 = recip32.wrapping_sub(1);
358
359 // We need to perform one more iteration to get us to 56 binary digits;
360 // The last iteration needs to happen with extra precision.
361 let q63blo = CastInto::<u32>::cast(b_significand << 11.cast());
362 let correction: u64;
363 let mut reciprocal: u64;
364 correction = negate_u64(
365 (recip32 as u64)
366 .wrapping_mul(q31b as u64)
367 .wrapping_add((recip32 as u64).wrapping_mul(q63blo as u64) >> 32),
368 );
369 let c_hi = (correction >> 32) as u32;
370 let c_lo = correction as u32;
371 reciprocal = (recip32 as u64)
372 .wrapping_mul(c_hi as u64)
373 .wrapping_add((recip32 as u64).wrapping_mul(c_lo as u64) >> 32);
374
375 // We already adjusted the 32-bit estimate, now we need to adjust the final
376 // 64-bit reciprocal estimate downward to ensure that it is strictly smaller
377 // than the infinitely precise exact reciprocal. Because the computation
378 // of the Newton-Raphson step is truncating at every step, this adjustment
379 // is small; most of the work is already done.
380 reciprocal = reciprocal.wrapping_sub(2);
381
382 // The numerical reciprocal is accurate to within 2^-56, lies in the
383 // interval [0.5, 1.0), and is strictly smaller than the true reciprocal
384 // of b. Multiplying a by this reciprocal thus gives a numerical q = a/b
385 // in Q53 with the following properties:
386 //
387 // 1. q < a/b
388 // 2. q is in the interval [0.5, 2.0)
389 // 3. the error in q is bounded away from 2^-53 (actually, we have a
390 // couple of bits to spare, but this is all we need).
391
392 // We need a 64 x 64 multiply high to compute q, which isn't a basic
393 // operation in C, so we need to be a little bit fussy.
394 // let mut quotient: F::Int = ((((reciprocal as u64)
395 // .wrapping_mul(CastInto::<u32>::cast(a_significand << 1) as u64))
396 // >> 32) as u32)
397 // .cast();
398
399 // We need a 64 x 64 multiply high to compute q, which isn't a basic
400 // operation in C, so we need to be a little bit fussy.
401 let mut quotient = (a_significand << 2).widen_mul(reciprocal.cast()).hi();
402
403 // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
404 // In either case, we are going to compute a residual of the form
405 //
406 // r = a - q*b
407 //
408 // We know from the construction of q that r satisfies:
409 //
410 // 0 <= r < ulp(q)*b
411 //
412 // if r is greater than 1/2 ulp(q)*b, then q rounds up. Otherwise, we
413 // already have the correct result. The exact halfway case cannot occur.
414 // We also take this time to right shift quotient if it falls in the [1,2)
415 // range and adjust the exponent accordingly.
416 let residual = if quotient < (implicit_bit << 1) {
417 quotient_exponent = quotient_exponent.wrapping_sub(1);
418 (a_significand << (significand_bits + 1)).wrapping_sub(quotient.wrapping_mul(b_significand))
419 } else {
420 quotient >>= 1;
421 (a_significand << significand_bits).wrapping_sub(quotient.wrapping_mul(b_significand))
422 };
423
424 let written_exponent = quotient_exponent.wrapping_add(exponent_bias as i32);
425
426 if written_exponent >= max_exponent as i32 {
427 // If we have overflowed the exponent, return infinity.
428 return F::from_repr(inf_rep | quotient_sign);
429 } else if written_exponent < 1 {
430 // Flush denormals to zero. In the future, it would be nice to add
431 // code to round them correctly.
432 return F::from_repr(quotient_sign);
433 } else {
434 let round = ((residual << 1) > b_significand) as u32;
435 // Clear the implicit bits
436 let mut abs_result = quotient & significand_mask;
437 // Insert the exponent
438 abs_result |= written_exponent.cast() << significand_bits;
439 // Round
440 abs_result = abs_result.wrapping_add(round.cast());
441 // Insert the sign and return
442 return F::from_repr(abs_result | quotient_sign);
443 }
444 }
445
446 intrinsics! {
447 #[arm_aeabi_alias = __aeabi_fdiv]
448 pub extern "C" fn __divsf3(a: f32, b: f32) -> f32 {
449 div32(a, b)
450 }
451
452 #[arm_aeabi_alias = __aeabi_ddiv]
453 pub extern "C" fn __divdf3(a: f64, b: f64) -> f64 {
454 div64(a, b)
455 }
456
457 #[cfg(target_arch = "arm")]
458 pub extern "C" fn __divsf3vfp(a: f32, b: f32) -> f32 {
459 a / b
460 }
461
462 #[cfg(target_arch = "arm")]
463 pub extern "C" fn __divdf3vfp(a: f64, b: f64) -> f64 {
464 a / b
465 }
466 }
467