Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit 93e3091

Browse files
authoredNov 14, 2017
Merge pull request #206 from aeleos/div
Implement divsf3 and divdf3
2 parents bb2c81b + 8bb3002 commit 93e3091

File tree

6 files changed

+659
-2
lines changed

6 files changed

+659
-2
lines changed
 

‎README.md

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -133,11 +133,11 @@ features = ["c"]
133133
- [ ] arm/unordsf2vfp.S
134134
- [x] ashldi3.c
135135
- [x] ashrdi3.c
136-
- [ ] divdf3.c
136+
- [x] divdf3.c
137137
- [x] divdi3.c
138138
- [x] divmoddi4.c
139139
- [x] divmodsi4.c
140-
- [ ] divsf3.c
140+
- [x] divsf3.c
141141
- [x] divsi3.c
142142
- [ ] extendhfsf2.c
143143
- [ ] extendsfdf2.c

‎build.rs

Lines changed: 185 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -125,6 +125,10 @@ mod tests {
125125
Mulsf3,
126126
Muldf3,
127127

128+
// float/div.rs
129+
Divsf3,
130+
Divdf3,
131+
128132
// int/mul.rs
129133
Muldi3,
130134
Mulodi4,
@@ -3399,6 +3403,187 @@ fn muldf3() {
33993403
}
34003404
}
34013405

3406+
#[derive(Eq, Hash, PartialEq)]
3407+
pub struct Divsf3 {
3408+
a: u32, // f32
3409+
b: u32, // f32
3410+
c: u32, // f32
3411+
}
3412+
3413+
impl TestCase for Divsf3 {
3414+
fn name() -> &'static str {
3415+
"divsf3"
3416+
}
3417+
3418+
fn generate<R>(rng: &mut R) -> Option<Self>
3419+
where
3420+
R: Rng,
3421+
Self: Sized,
3422+
{
3423+
let a = gen_large_f32(rng);
3424+
let b = gen_large_f32(rng);
3425+
if b == 0.0 {
3426+
return None;
3427+
}
3428+
let c = a / b;
3429+
// TODO accept NaNs. We don't do that right now because we can't check
3430+
// for NaN-ness on the thumb targets (due to missing intrinsics)
3431+
if a.is_nan() || b.is_nan() || c.is_nan()|| c.abs() <= unsafe { mem::transmute(16777215u32) } {
3432+
return None;
3433+
}
3434+
3435+
Some(
3436+
Divsf3 {
3437+
a: to_u32(a),
3438+
b: to_u32(b),
3439+
c: to_u32(c),
3440+
},
3441+
)
3442+
}
3443+
3444+
fn to_string(&self, buffer: &mut String) {
3445+
writeln!(
3446+
buffer,
3447+
"(({a}, {b}), {c}),",
3448+
a = self.a,
3449+
b = self.b,
3450+
c = self.c
3451+
)
3452+
.unwrap();
3453+
}
3454+
3455+
fn prologue() -> &'static str {
3456+
r#"
3457+
#[cfg(all(target_arch = "arm",
3458+
not(any(target_env = "gnu", target_env = "musl")),
3459+
target_os = "linux",
3460+
test))]
3461+
use core::mem;
3462+
#[cfg(not(all(target_arch = "arm",
3463+
not(any(target_env = "gnu", target_env = "musl")),
3464+
target_os = "linux",
3465+
test)))]
3466+
use std::mem;
3467+
use compiler_builtins::float::div::__divsf3;
3468+
3469+
fn mk_f32(x: u32) -> f32 {
3470+
unsafe { mem::transmute(x) }
3471+
}
3472+
3473+
fn to_u32(x: f32) -> u32 {
3474+
unsafe { mem::transmute(x) }
3475+
}
3476+
3477+
static TEST_CASES: &[((u32, u32), u32)] = &[
3478+
"#
3479+
}
3480+
3481+
fn epilogue() -> &'static str {
3482+
"
3483+
];
3484+
3485+
#[test]
3486+
fn divsf3() {
3487+
for &((a, b), c) in TEST_CASES {
3488+
let c_ = __divsf3(mk_f32(a), mk_f32(b));
3489+
assert_eq!(((a, b), c), ((a, b), to_u32(c_)));
3490+
}
3491+
}
3492+
"
3493+
}
3494+
}
3495+
3496+
#[derive(Eq, Hash, PartialEq)]
3497+
pub struct Divdf3 {
3498+
a: u64, // f64
3499+
b: u64, // f64
3500+
c: u64, // f64
3501+
}
3502+
3503+
impl TestCase for Divdf3 {
3504+
fn name() -> &'static str {
3505+
"divdf3"
3506+
}
3507+
3508+
fn generate<R>(rng: &mut R) -> Option<Self>
3509+
where
3510+
R: Rng,
3511+
Self: Sized,
3512+
{
3513+
let a = gen_large_f64(rng);
3514+
let b = gen_large_f64(rng);
3515+
if b == 0.0 {
3516+
return None;
3517+
}
3518+
let c = a / b;
3519+
// TODO accept NaNs. We don't do that right now because we can't check
3520+
// for NaN-ness on the thumb targets (due to missing intrinsics)
3521+
if a.is_nan() || b.is_nan() || c.is_nan()
3522+
|| c.abs() <= unsafe { mem::transmute(4503599627370495u64) } {
3523+
return None;
3524+
}
3525+
3526+
Some(
3527+
Divdf3 {
3528+
a: to_u64(a),
3529+
b: to_u64(b),
3530+
c: to_u64(c),
3531+
},
3532+
)
3533+
}
3534+
3535+
fn to_string(&self, buffer: &mut String) {
3536+
writeln!(
3537+
buffer,
3538+
"(({a}, {b}), {c}),",
3539+
a = self.a,
3540+
b = self.b,
3541+
c = self.c
3542+
)
3543+
.unwrap();
3544+
}
3545+
3546+
fn prologue() -> &'static str {
3547+
r#"
3548+
#[cfg(all(target_arch = "arm",
3549+
not(any(target_env = "gnu", target_env = "musl")),
3550+
target_os = "linux",
3551+
test))]
3552+
use core::mem;
3553+
#[cfg(not(all(target_arch = "arm",
3554+
not(any(target_env = "gnu", target_env = "musl")),
3555+
target_os = "linux",
3556+
test)))]
3557+
use std::mem;
3558+
use compiler_builtins::float::div::__divdf3;
3559+
3560+
fn mk_f64(x: u64) -> f64 {
3561+
unsafe { mem::transmute(x) }
3562+
}
3563+
3564+
fn to_u64(x: f64) -> u64 {
3565+
unsafe { mem::transmute(x) }
3566+
}
3567+
3568+
static TEST_CASES: &[((u64, u64), u64)] = &[
3569+
"#
3570+
}
3571+
3572+
fn epilogue() -> &'static str {
3573+
"
3574+
];
3575+
3576+
#[test]
3577+
fn divdf3() {
3578+
for &((a, b), c) in TEST_CASES {
3579+
let c_ = __divdf3(mk_f64(a), mk_f64(b));
3580+
assert_eq!(((a, b), c), ((a, b), to_u64(c_)));
3581+
}
3582+
}
3583+
"
3584+
}
3585+
}
3586+
34023587

34033588
#[derive(Eq, Hash, PartialEq)]
34043589
pub struct Udivdi3 {

‎src/float/div.rs

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

‎src/float/mod.rs

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ pub mod add;
88
pub mod pow;
99
pub mod sub;
1010
pub mod mul;
11+
pub mod div;
1112

1213
/// Trait for some basic operations on floats
1314
pub trait Float:

‎tests/divdf3.rs

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
#![feature(compiler_builtins_lib)]
2+
#![feature(i128_type)]
3+
#![cfg_attr(all(target_arch = "arm", not(any(target_env = "gnu", target_env = "musl")),
4+
target_os = "linux", test),
5+
no_std)]
6+
7+
include!(concat!(env!("OUT_DIR"), "/divdf3.rs"));

‎tests/divsf3.rs

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
#![feature(compiler_builtins_lib)]
2+
#![feature(i128_type)]
3+
#![cfg_attr(all(target_arch = "arm", not(any(target_env = "gnu", target_env = "musl")),
4+
target_os = "linux", test),
5+
no_std)]
6+
7+
include!(concat!(env!("OUT_DIR"), "/divsf3.rs"));

0 commit comments

Comments
 (0)
Please sign in to comment.