| 1 | use crate::float::Float; |
| 2 | use crate::int::{CastInto, Int, MinInt}; |
| 3 | |
| 4 | /// Returns `a + b` |
| 5 | fn add<F: Float>(a: F, b: F) -> F |
| 6 | where |
| 7 | u32: CastInto<F::Int>, |
| 8 | F::Int: CastInto<u32>, |
| 9 | i32: CastInto<F::Int>, |
| 10 | F::Int: CastInto<i32>, |
| 11 | { |
| 12 | let one = F::Int::ONE; |
| 13 | let zero = F::Int::ZERO; |
| 14 | |
| 15 | let bits = F::BITS.cast(); |
| 16 | let significand_bits = F::SIG_BITS; |
| 17 | let max_exponent = F::EXP_SAT; |
| 18 | |
| 19 | let implicit_bit = F::IMPLICIT_BIT; |
| 20 | let significand_mask = F::SIG_MASK; |
| 21 | let sign_bit = F::SIGN_MASK as F::Int; |
| 22 | let abs_mask = sign_bit - one; |
| 23 | let exponent_mask = F::EXP_MASK; |
| 24 | let inf_rep = exponent_mask; |
| 25 | let quiet_bit = implicit_bit >> 1; |
| 26 | let qnan_rep = exponent_mask | quiet_bit; |
| 27 | |
| 28 | let mut a_rep = a.to_bits(); |
| 29 | let mut b_rep = b.to_bits(); |
| 30 | let a_abs = a_rep & abs_mask; |
| 31 | let b_abs = b_rep & abs_mask; |
| 32 | |
| 33 | // Detect if a or b is zero, infinity, or NaN. |
| 34 | if a_abs.wrapping_sub(one) >= inf_rep - one || b_abs.wrapping_sub(one) >= inf_rep - one { |
| 35 | // NaN + anything = qNaN |
| 36 | if a_abs > inf_rep { |
| 37 | return F::from_bits(a_abs | quiet_bit); |
| 38 | } |
| 39 | // anything + NaN = qNaN |
| 40 | if b_abs > inf_rep { |
| 41 | return F::from_bits(b_abs | quiet_bit); |
| 42 | } |
| 43 | |
| 44 | if a_abs == inf_rep { |
| 45 | // +/-infinity + -/+infinity = qNaN |
| 46 | if (a.to_bits() ^ b.to_bits()) == sign_bit { |
| 47 | return F::from_bits(qnan_rep); |
| 48 | } else { |
| 49 | // +/-infinity + anything remaining = +/- infinity |
| 50 | return a; |
| 51 | } |
| 52 | } |
| 53 | |
| 54 | // anything remaining + +/-infinity = +/-infinity |
| 55 | if b_abs == inf_rep { |
| 56 | return b; |
| 57 | } |
| 58 | |
| 59 | // zero + anything = anything |
| 60 | if a_abs == MinInt::ZERO { |
| 61 | // but we need to get the sign right for zero + zero |
| 62 | if b_abs == MinInt::ZERO { |
| 63 | return F::from_bits(a.to_bits() & b.to_bits()); |
| 64 | } else { |
| 65 | return b; |
| 66 | } |
| 67 | } |
| 68 | |
| 69 | // anything + zero = anything |
| 70 | if b_abs == MinInt::ZERO { |
| 71 | return a; |
| 72 | } |
| 73 | } |
| 74 | |
| 75 | // Swap a and b if necessary so that a has the larger absolute value. |
| 76 | if b_abs > a_abs { |
| 77 | // Don't use mem::swap because it may generate references to memcpy in unoptimized code. |
| 78 | let tmp = a_rep; |
| 79 | a_rep = b_rep; |
| 80 | b_rep = tmp; |
| 81 | } |
| 82 | |
| 83 | // Extract the exponent and significand from the (possibly swapped) a and b. |
| 84 | let mut a_exponent: i32 = ((a_rep & exponent_mask) >> significand_bits).cast(); |
| 85 | let mut b_exponent: i32 = ((b_rep & exponent_mask) >> significand_bits).cast(); |
| 86 | let mut a_significand = a_rep & significand_mask; |
| 87 | let mut b_significand = b_rep & significand_mask; |
| 88 | |
| 89 | // normalize any denormals, and adjust the exponent accordingly. |
| 90 | if a_exponent == 0 { |
| 91 | let (exponent, significand) = F::normalize(a_significand); |
| 92 | a_exponent = exponent; |
| 93 | a_significand = significand; |
| 94 | } |
| 95 | if b_exponent == 0 { |
| 96 | let (exponent, significand) = F::normalize(b_significand); |
| 97 | b_exponent = exponent; |
| 98 | b_significand = significand; |
| 99 | } |
| 100 | |
| 101 | // The sign of the result is the sign of the larger operand, a. If they |
| 102 | // have opposite signs, we are performing a subtraction; otherwise addition. |
| 103 | let result_sign = a_rep & sign_bit; |
| 104 | let subtraction = ((a_rep ^ b_rep) & sign_bit) != zero; |
| 105 | |
| 106 | // Shift the significands to give us round, guard and sticky, and or in the |
| 107 | // implicit significand bit. (If we fell through from the denormal path it |
| 108 | // was already set by normalize(), but setting it twice won't hurt |
| 109 | // anything.) |
| 110 | a_significand = (a_significand | implicit_bit) << 3; |
| 111 | b_significand = (b_significand | implicit_bit) << 3; |
| 112 | |
| 113 | // Shift the significand of b by the difference in exponents, with a sticky |
| 114 | // bottom bit to get rounding correct. |
| 115 | let align = a_exponent.wrapping_sub(b_exponent).cast(); |
| 116 | if align != MinInt::ZERO { |
| 117 | if align < bits { |
| 118 | let sticky = |
| 119 | F::Int::from_bool(b_significand << bits.wrapping_sub(align).cast() != MinInt::ZERO); |
| 120 | b_significand = (b_significand >> align.cast()) | sticky; |
| 121 | } else { |
| 122 | b_significand = one; // sticky; b is known to be non-zero. |
| 123 | } |
| 124 | } |
| 125 | if subtraction { |
| 126 | a_significand = a_significand.wrapping_sub(b_significand); |
| 127 | // If a == -b, return +zero. |
| 128 | if a_significand == MinInt::ZERO { |
| 129 | return F::from_bits(MinInt::ZERO); |
| 130 | } |
| 131 | |
| 132 | // If partial cancellation occured, we need to left-shift the result |
| 133 | // and adjust the exponent: |
| 134 | if a_significand < implicit_bit << 3 { |
| 135 | let shift = |
| 136 | a_significand.leading_zeros() as i32 - (implicit_bit << 3).leading_zeros() as i32; |
| 137 | a_significand <<= shift; |
| 138 | a_exponent -= shift; |
| 139 | } |
| 140 | } else { |
| 141 | // addition |
| 142 | a_significand += b_significand; |
| 143 | |
| 144 | // If the addition carried up, we need to right-shift the result and |
| 145 | // adjust the exponent: |
| 146 | if a_significand & (implicit_bit << 4) != MinInt::ZERO { |
| 147 | let sticky = F::Int::from_bool(a_significand & one != MinInt::ZERO); |
| 148 | a_significand = (a_significand >> 1) | sticky; |
| 149 | a_exponent += 1; |
| 150 | } |
| 151 | } |
| 152 | |
| 153 | // If we have overflowed the type, return +/- infinity: |
| 154 | if a_exponent >= max_exponent as i32 { |
| 155 | return F::from_bits(inf_rep | result_sign); |
| 156 | } |
| 157 | |
| 158 | if a_exponent <= 0 { |
| 159 | // Result is denormal before rounding; the exponent is zero and we |
| 160 | // need to shift the significand. |
| 161 | let shift = (1 - a_exponent).cast(); |
| 162 | let sticky = |
| 163 | F::Int::from_bool((a_significand << bits.wrapping_sub(shift).cast()) != MinInt::ZERO); |
| 164 | a_significand = (a_significand >> shift.cast()) | sticky; |
| 165 | a_exponent = 0; |
| 166 | } |
| 167 | |
| 168 | // Low three bits are round, guard, and sticky. |
| 169 | let a_significand_i32: i32 = a_significand.cast(); |
| 170 | let round_guard_sticky: i32 = a_significand_i32 & 0x7; |
| 171 | |
| 172 | // Shift the significand into place, and mask off the implicit bit. |
| 173 | let mut result = (a_significand >> 3) & significand_mask; |
| 174 | |
| 175 | // Insert the exponent and sign. |
| 176 | result |= a_exponent.cast() << significand_bits; |
| 177 | result |= result_sign; |
| 178 | |
| 179 | // Final rounding. The result may overflow to infinity, but that is the |
| 180 | // correct result in that case. |
| 181 | if round_guard_sticky > 0x4 { |
| 182 | result += one; |
| 183 | } |
| 184 | if round_guard_sticky == 0x4 { |
| 185 | result += result & one; |
| 186 | } |
| 187 | |
| 188 | F::from_bits(result) |
| 189 | } |
| 190 | |
| 191 | intrinsics! { |
| 192 | #[avr_skip] |
| 193 | #[aapcs_on_arm] |
| 194 | #[arm_aeabi_alias = __aeabi_fadd] |
| 195 | pub extern "C" fn __addsf3(a: f32, b: f32) -> f32 { |
| 196 | add(a, b) |
| 197 | } |
| 198 | |
| 199 | #[avr_skip] |
| 200 | #[aapcs_on_arm] |
| 201 | #[arm_aeabi_alias = __aeabi_dadd] |
| 202 | pub extern "C" fn __adddf3(a: f64, b: f64) -> f64 { |
| 203 | add(a, b) |
| 204 | } |
| 205 | |
| 206 | #[ppc_alias = __addkf3] |
| 207 | #[cfg (f128_enabled)] |
| 208 | pub extern "C" fn __addtf3(a: f128, b: f128) -> f128 { |
| 209 | add(a, b) |
| 210 | } |
| 211 | } |
| 212 | |