1 | use crate::float::Float; |
2 | use crate::int::{CastInto, Int}; |
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::SIGNIFICAND_BITS; |
17 | let max_exponent = F::EXPONENT_MAX; |
18 | |
19 | let implicit_bit = F::IMPLICIT_BIT; |
20 | let significand_mask = F::SIGNIFICAND_MASK; |
21 | let sign_bit = F::SIGN_MASK as F::Int; |
22 | let abs_mask = sign_bit - one; |
23 | let exponent_mask = F::EXPONENT_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.repr(); |
29 | let mut b_rep = b.repr(); |
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_repr(a_abs | quiet_bit); |
38 | } |
39 | // anything + NaN = qNaN |
40 | if b_abs > inf_rep { |
41 | return F::from_repr(b_abs | quiet_bit); |
42 | } |
43 | |
44 | if a_abs == inf_rep { |
45 | // +/-infinity + -/+infinity = qNaN |
46 | if (a.repr() ^ b.repr()) == sign_bit { |
47 | return F::from_repr(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 == Int::ZERO { |
61 | // but we need to get the sign right for zero + zero |
62 | if b_abs == Int::ZERO { |
63 | return F::from_repr(a.repr() & b.repr()); |
64 | } else { |
65 | return b; |
66 | } |
67 | } |
68 | |
69 | // anything + zero = anything |
70 | if b_abs == Int::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 != Int::ZERO { |
117 | if align < bits { |
118 | let sticky = |
119 | F::Int::from_bool(b_significand << bits.wrapping_sub(align).cast() != Int::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 == Int::ZERO { |
129 | return F::from_repr(Int::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 != Int::ZERO { |
147 | let sticky = F::Int::from_bool(a_significand & one != Int::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_repr(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()) != Int::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_repr(result) |
189 | } |
190 | |
191 | intrinsics! { |
192 | #[aapcs_on_arm] |
193 | #[arm_aeabi_alias = __aeabi_fadd] |
194 | pub extern "C" fn __addsf3(a: f32, b: f32) -> f32 { |
195 | add(a, b) |
196 | } |
197 | |
198 | #[aapcs_on_arm] |
199 | #[arm_aeabi_alias = __aeabi_dadd] |
200 | pub extern "C" fn __adddf3(a: f64, b: f64) -> f64 { |
201 | add(a, b) |
202 | } |
203 | |
204 | #[cfg (target_arch = "arm" )] |
205 | pub extern "C" fn __addsf3vfp(a: f32, b: f32) -> f32 { |
206 | a + b |
207 | } |
208 | |
209 | #[cfg (target_arch = "arm" )] |
210 | pub extern "C" fn __adddf3vfp(a: f64, b: f64) -> f64 { |
211 | a + b |
212 | } |
213 | } |
214 | |