2016-08-18 04:51:37 +08:00
|
|
|
use float::Float;
|
2019-05-15 05:33:08 +08:00
|
|
|
use int::{CastInto, Int};
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-06-24 02:05:25 +08:00
|
|
|
/// Returns `a + b`
|
2019-05-15 05:33:08 +08:00
|
|
|
fn add<F: Float>(a: F, b: F) -> F
|
|
|
|
where
|
2017-09-14 23:33:44 +08:00
|
|
|
u32: CastInto<F::Int>,
|
|
|
|
F::Int: CastInto<u32>,
|
|
|
|
i32: CastInto<F::Int>,
|
|
|
|
F::Int: CastInto<i32>,
|
|
|
|
{
|
|
|
|
let one = F::Int::ONE;
|
|
|
|
let zero = F::Int::ZERO;
|
|
|
|
|
2019-05-15 05:33:08 +08:00
|
|
|
let bits = F::BITS.cast();
|
2017-09-14 23:33:44 +08:00
|
|
|
let significand_bits = F::SIGNIFICAND_BITS;
|
2019-05-15 05:33:08 +08:00
|
|
|
let max_exponent = F::EXPONENT_MAX;
|
2017-09-14 23:33:44 +08:00
|
|
|
|
2019-05-15 05:33:08 +08:00
|
|
|
let implicit_bit = F::IMPLICIT_BIT;
|
2017-09-14 23:33:44 +08:00
|
|
|
let significand_mask = F::SIGNIFICAND_MASK;
|
2019-05-15 05:33:08 +08:00
|
|
|
let sign_bit = F::SIGN_MASK as F::Int;
|
|
|
|
let abs_mask = sign_bit - one;
|
|
|
|
let exponent_mask = F::EXPONENT_MASK;
|
|
|
|
let inf_rep = exponent_mask;
|
|
|
|
let quiet_bit = implicit_bit >> 1;
|
|
|
|
let qnan_rep = exponent_mask | quiet_bit;
|
2017-09-14 23:33:44 +08:00
|
|
|
|
|
|
|
let mut a_rep = a.repr();
|
|
|
|
let mut b_rep = b.repr();
|
|
|
|
let a_abs = a_rep & abs_mask;
|
|
|
|
let b_abs = b_rep & abs_mask;
|
|
|
|
|
|
|
|
// Detect if a or b is zero, infinity, or NaN.
|
2019-05-15 05:33:08 +08:00
|
|
|
if a_abs.wrapping_sub(one) >= inf_rep - one || b_abs.wrapping_sub(one) >= inf_rep - one {
|
2017-09-14 23:33:44 +08:00
|
|
|
// NaN + anything = qNaN
|
|
|
|
if a_abs > inf_rep {
|
|
|
|
return F::from_repr(a_abs | quiet_bit);
|
|
|
|
}
|
|
|
|
// anything + NaN = qNaN
|
|
|
|
if b_abs > inf_rep {
|
|
|
|
return F::from_repr(b_abs | quiet_bit);
|
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
if a_abs == inf_rep {
|
|
|
|
// +/-infinity + -/+infinity = qNaN
|
|
|
|
if (a.repr() ^ b.repr()) == sign_bit {
|
|
|
|
return F::from_repr(qnan_rep);
|
|
|
|
} else {
|
|
|
|
// +/-infinity + anything remaining = +/- infinity
|
|
|
|
return a;
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2017-09-14 23:33:44 +08:00
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// anything remaining + +/-infinity = +/-infinity
|
|
|
|
if b_abs == inf_rep {
|
|
|
|
return b;
|
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// zero + anything = anything
|
|
|
|
if a_abs == Int::ZERO {
|
|
|
|
// but we need to get the sign right for zero + zero
|
|
|
|
if b_abs == Int::ZERO {
|
|
|
|
return F::from_repr(a.repr() & b.repr());
|
|
|
|
} else {
|
|
|
|
return b;
|
2016-08-18 04:51:37 +08:00
|
|
|
}
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// anything + zero = anything
|
|
|
|
if b_abs == Int::ZERO {
|
2019-05-15 05:33:08 +08:00
|
|
|
return a;
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2017-09-14 23:33:44 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
// Swap a and b if necessary so that a has the larger absolute value.
|
|
|
|
if b_abs > a_abs {
|
|
|
|
// Don't use mem::swap because it may generate references to memcpy in unoptimized code.
|
|
|
|
let tmp = a_rep;
|
|
|
|
a_rep = b_rep;
|
|
|
|
b_rep = tmp;
|
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// Extract the exponent and significand from the (possibly swapped) a and b.
|
|
|
|
let mut a_exponent: i32 = ((a_rep & exponent_mask) >> significand_bits).cast();
|
|
|
|
let mut b_exponent: i32 = ((b_rep & exponent_mask) >> significand_bits).cast();
|
|
|
|
let mut a_significand = a_rep & significand_mask;
|
|
|
|
let mut b_significand = b_rep & significand_mask;
|
|
|
|
|
|
|
|
// normalize any denormals, and adjust the exponent accordingly.
|
|
|
|
if a_exponent == 0 {
|
|
|
|
let (exponent, significand) = F::normalize(a_significand);
|
|
|
|
a_exponent = exponent;
|
|
|
|
a_significand = significand;
|
|
|
|
}
|
|
|
|
if b_exponent == 0 {
|
|
|
|
let (exponent, significand) = F::normalize(b_significand);
|
|
|
|
b_exponent = exponent;
|
|
|
|
b_significand = significand;
|
|
|
|
}
|
|
|
|
|
|
|
|
// The sign of the result is the sign of the larger operand, a. If they
|
|
|
|
// have opposite signs, we are performing a subtraction; otherwise addition.
|
|
|
|
let result_sign = a_rep & sign_bit;
|
|
|
|
let subtraction = ((a_rep ^ b_rep) & sign_bit) != zero;
|
|
|
|
|
|
|
|
// Shift the significands to give us round, guard and sticky, and or in the
|
|
|
|
// implicit significand bit. (If we fell through from the denormal path it
|
|
|
|
// was already set by normalize(), but setting it twice won't hurt
|
|
|
|
// anything.)
|
|
|
|
a_significand = (a_significand | implicit_bit) << 3;
|
|
|
|
b_significand = (b_significand | implicit_bit) << 3;
|
|
|
|
|
|
|
|
// Shift the significand of b by the difference in exponents, with a sticky
|
|
|
|
// bottom bit to get rounding correct.
|
|
|
|
let align = a_exponent.wrapping_sub(b_exponent).cast();
|
|
|
|
if align != Int::ZERO {
|
|
|
|
if align < bits {
|
2019-05-15 05:33:08 +08:00
|
|
|
let sticky =
|
|
|
|
F::Int::from_bool(b_significand << bits.wrapping_sub(align).cast() != Int::ZERO);
|
2017-09-14 23:33:44 +08:00
|
|
|
b_significand = (b_significand >> align.cast()) | sticky;
|
|
|
|
} else {
|
|
|
|
b_significand = one; // sticky; b is known to be non-zero.
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2017-09-14 23:33:44 +08:00
|
|
|
}
|
|
|
|
if subtraction {
|
|
|
|
a_significand = a_significand.wrapping_sub(b_significand);
|
|
|
|
// If a == -b, return +zero.
|
|
|
|
if a_significand == Int::ZERO {
|
|
|
|
return F::from_repr(Int::ZERO);
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// If partial cancellation occured, we need to left-shift the result
|
|
|
|
// and adjust the exponent:
|
|
|
|
if a_significand < implicit_bit << 3 {
|
2019-05-15 05:33:08 +08:00
|
|
|
let shift =
|
|
|
|
a_significand.leading_zeros() as i32 - (implicit_bit << 3).leading_zeros() as i32;
|
2017-09-14 23:33:44 +08:00
|
|
|
a_significand <<= shift;
|
|
|
|
a_exponent -= shift;
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2019-05-15 05:33:08 +08:00
|
|
|
} else
|
|
|
|
/* addition */
|
|
|
|
{
|
2017-09-14 23:33:44 +08:00
|
|
|
a_significand += b_significand;
|
|
|
|
|
|
|
|
// If the addition carried up, we need to right-shift the result and
|
|
|
|
// adjust the exponent:
|
|
|
|
if a_significand & implicit_bit << 4 != Int::ZERO {
|
|
|
|
let sticky = F::Int::from_bool(a_significand & one != Int::ZERO);
|
|
|
|
a_significand = a_significand >> 1 | sticky;
|
|
|
|
a_exponent += 1;
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2017-09-14 23:33:44 +08:00
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// If we have overflowed the type, return +/- infinity:
|
|
|
|
if a_exponent >= max_exponent as i32 {
|
|
|
|
return F::from_repr(inf_rep | result_sign);
|
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
if a_exponent <= 0 {
|
|
|
|
// Result is denormal before rounding; the exponent is zero and we
|
|
|
|
// need to shift the significand.
|
|
|
|
let shift = (1 - a_exponent).cast();
|
2019-05-15 05:33:08 +08:00
|
|
|
let sticky =
|
|
|
|
F::Int::from_bool((a_significand << bits.wrapping_sub(shift).cast()) != Int::ZERO);
|
2017-09-14 23:33:44 +08:00
|
|
|
a_significand = a_significand >> shift.cast() | sticky;
|
|
|
|
a_exponent = 0;
|
|
|
|
}
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// Low three bits are round, guard, and sticky.
|
|
|
|
let a_significand_i32: i32 = a_significand.cast();
|
|
|
|
let round_guard_sticky: i32 = a_significand_i32 & 0x7;
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// Shift the significand into place, and mask off the implicit bit.
|
|
|
|
let mut result = a_significand >> 3 & significand_mask;
|
2016-10-16 19:53:31 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// Insert the exponent and sign.
|
|
|
|
result |= a_exponent.cast() << significand_bits;
|
|
|
|
result |= result_sign;
|
2016-08-18 04:51:37 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
// Final rounding. The result may overflow to infinity, but that is the
|
|
|
|
// correct result in that case.
|
2019-05-15 05:33:08 +08:00
|
|
|
if round_guard_sticky > 0x4 {
|
|
|
|
result += one;
|
|
|
|
}
|
|
|
|
if round_guard_sticky == 0x4 {
|
|
|
|
result += result & one;
|
|
|
|
}
|
2017-02-08 01:51:18 +08:00
|
|
|
|
2017-09-14 23:33:44 +08:00
|
|
|
F::from_repr(result)
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2017-02-08 01:51:18 +08:00
|
|
|
|
2017-06-24 02:05:25 +08:00
|
|
|
intrinsics! {
|
|
|
|
#[aapcs_on_arm]
|
2017-06-24 06:16:07 +08:00
|
|
|
#[arm_aeabi_alias = __aeabi_fadd]
|
2017-06-24 02:05:25 +08:00
|
|
|
pub extern "C" fn __addsf3(a: f32, b: f32) -> f32 {
|
2017-09-14 23:33:44 +08:00
|
|
|
add(a, b)
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2017-02-08 01:51:18 +08:00
|
|
|
|
2017-06-24 02:05:25 +08:00
|
|
|
#[aapcs_on_arm]
|
2017-06-24 06:16:07 +08:00
|
|
|
#[arm_aeabi_alias = __aeabi_dadd]
|
2017-06-24 02:05:25 +08:00
|
|
|
pub extern "C" fn __adddf3(a: f64, b: f64) -> f64 {
|
2017-09-14 23:33:44 +08:00
|
|
|
add(a, b)
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|
2018-02-05 01:34:31 +08:00
|
|
|
|
|
|
|
#[cfg(target_arch = "arm")]
|
|
|
|
pub extern "C" fn __addsf3vfp(a: f32, b: f32) -> f32 {
|
|
|
|
a + b
|
|
|
|
}
|
|
|
|
|
|
|
|
#[cfg(target_arch = "arm")]
|
|
|
|
pub extern "C" fn __adddf3vfp(a: f64, b: f64) -> f64 {
|
|
|
|
a + b
|
|
|
|
}
|
2017-06-24 02:05:25 +08:00
|
|
|
}
|