source: mainline/uspace/lib/softfloat/add.c@ 8c95dff

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since 8c95dff was faa45c17, checked in by Jakub Jermar <jakub@…>, 13 years ago

Fix a wrong index in softfloat1 test and a couple of cstyle softfloat fixes.
  • Property mode set to 100644
File size: 9.9 KB
Line 
1/*
2 * Copyright (c) 2005 Josef Cejka
3 * Copyright (c) 2011 Petr Koupy
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * - Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * - Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * - The name of the author may not be used to endorse or promote products
16 * derived from this software without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
19 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
20 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
22 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
23 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
25 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
26 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
27 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30/** @addtogroup softfloat
31 * @{
32 */
33/** @file Addition functions.
34 */
35
36#include "sftypes.h"
37#include "add.h"
38#include "comparison.h"
39#include "common.h"
40
41/** Add two single-precision floats with the same sign.
42 *
43 * @param a First input operand.
44 * @param b Second input operand.
45 * @return Result of addition.
46 */
47float32 add_float32(float32 a, float32 b)
48{
49 int expdiff;
50 uint32_t exp1, exp2, frac1, frac2;
51
52 expdiff = a.parts.exp - b.parts.exp;
53 if (expdiff < 0) {
54 if (is_float32_nan(b)) {
55 /* TODO: fix SigNaN */
56 if (is_float32_signan(b)) {
57 }
58
59 return b;
60 }
61
62 if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
63 return b;
64 }
65
66 frac1 = b.parts.fraction;
67 exp1 = b.parts.exp;
68 frac2 = a.parts.fraction;
69 exp2 = a.parts.exp;
70 expdiff *= -1;
71 } else {
72 if ((is_float32_nan(a)) || (is_float32_nan(b))) {
73 /* TODO: fix SigNaN */
74 if (is_float32_signan(a) || is_float32_signan(b)) {
75 }
76 return (is_float32_nan(a) ? a : b);
77 }
78
79 if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
80 return a;
81 }
82
83 frac1 = a.parts.fraction;
84 exp1 = a.parts.exp;
85 frac2 = b.parts.fraction;
86 exp2 = b.parts.exp;
87 }
88
89 if (exp1 == 0) {
90 /* both are denormalized */
91 frac1 += frac2;
92 if (frac1 & FLOAT32_HIDDEN_BIT_MASK) {
93 /* result is not denormalized */
94 a.parts.exp = 1;
95 }
96 a.parts.fraction = frac1;
97 return a;
98 }
99
100 frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
101
102 if (exp2 == 0) {
103 /* second operand is denormalized */
104 --expdiff;
105 } else {
106 /* add hidden bit to second operand */
107 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
108 }
109
110 /* create some space for rounding */
111 frac1 <<= 6;
112 frac2 <<= 6;
113
114 if (expdiff < (FLOAT32_FRACTION_SIZE + 2)) {
115 frac2 >>= expdiff;
116 frac1 += frac2;
117 } else {
118 a.parts.exp = exp1;
119 a.parts.fraction = (frac1 >> 6) & (~(FLOAT32_HIDDEN_BIT_MASK));
120 return a;
121 }
122
123 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
124 ++exp1;
125 frac1 >>= 1;
126 }
127
128 /* rounding - if first bit after fraction is set then round up */
129 frac1 += (0x1 << 5);
130
131 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
132 /* rounding overflow */
133 ++exp1;
134 frac1 >>= 1;
135 }
136
137 if ((exp1 == FLOAT32_MAX_EXPONENT) || (exp2 > exp1)) {
138 /* overflow - set infinity as result */
139 a.parts.exp = FLOAT32_MAX_EXPONENT;
140 a.parts.fraction = 0;
141 return a;
142 }
143
144 a.parts.exp = exp1;
145
146 /* Clear hidden bit and shift */
147 a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
148 return a;
149}
150
151/** Add two double-precision floats with the same sign.
152 *
153 * @param a First input operand.
154 * @param b Second input operand.
155 * @return Result of addition.
156 */
157float64 add_float64(float64 a, float64 b)
158{
159 int expdiff;
160 uint32_t exp1, exp2;
161 uint64_t frac1, frac2;
162
163 expdiff = ((int) a.parts.exp) - b.parts.exp;
164 if (expdiff < 0) {
165 if (is_float64_nan(b)) {
166 /* TODO: fix SigNaN */
167 if (is_float64_signan(b)) {
168 }
169
170 return b;
171 }
172
173 /* b is infinity and a not */
174 if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
175 return b;
176 }
177
178 frac1 = b.parts.fraction;
179 exp1 = b.parts.exp;
180 frac2 = a.parts.fraction;
181 exp2 = a.parts.exp;
182 expdiff *= -1;
183 } else {
184 if (is_float64_nan(a)) {
185 /* TODO: fix SigNaN */
186 if (is_float64_signan(a) || is_float64_signan(b)) {
187 }
188 return a;
189 }
190
191 /* a is infinity and b not */
192 if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
193 return a;
194 }
195
196 frac1 = a.parts.fraction;
197 exp1 = a.parts.exp;
198 frac2 = b.parts.fraction;
199 exp2 = b.parts.exp;
200 }
201
202 if (exp1 == 0) {
203 /* both are denormalized */
204 frac1 += frac2;
205 if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
206 /* result is not denormalized */
207 a.parts.exp = 1;
208 }
209 a.parts.fraction = frac1;
210 return a;
211 }
212
213 /* add hidden bit - frac1 is sure not denormalized */
214 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
215
216 /* second operand ... */
217 if (exp2 == 0) {
218 /* ... is denormalized */
219 --expdiff;
220 } else {
221 /* is not denormalized */
222 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
223 }
224
225 /* create some space for rounding */
226 frac1 <<= 6;
227 frac2 <<= 6;
228
229 if (expdiff < (FLOAT64_FRACTION_SIZE + 2)) {
230 frac2 >>= expdiff;
231 frac1 += frac2;
232 } else {
233 a.parts.exp = exp1;
234 a.parts.fraction = (frac1 >> 6) & (~(FLOAT64_HIDDEN_BIT_MASK));
235 return a;
236 }
237
238 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
239 ++exp1;
240 frac1 >>= 1;
241 }
242
243 /* rounding - if first bit after fraction is set then round up */
244 frac1 += (0x1 << 5);
245
246 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
247 /* rounding overflow */
248 ++exp1;
249 frac1 >>= 1;
250 }
251
252 if ((exp1 == FLOAT64_MAX_EXPONENT) || (exp2 > exp1)) {
253 /* overflow - set infinity as result */
254 a.parts.exp = FLOAT64_MAX_EXPONENT;
255 a.parts.fraction = 0;
256 return a;
257 }
258
259 a.parts.exp = exp1;
260 /* Clear hidden bit and shift */
261 a.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK));
262 return a;
263}
264
265/** Add two quadruple-precision floats with the same sign.
266 *
267 * @param a First input operand.
268 * @param b Second input operand.
269 * @return Result of addition.
270 */
271float128 add_float128(float128 a, float128 b)
272{
273 int expdiff;
274 uint32_t exp1, exp2;
275 uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
276
277 expdiff = ((int) a.parts.exp) - b.parts.exp;
278 if (expdiff < 0) {
279 if (is_float128_nan(b)) {
280 /* TODO: fix SigNaN */
281 if (is_float128_signan(b)) {
282 }
283
284 return b;
285 }
286
287 /* b is infinity and a not */
288 if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
289 return b;
290 }
291
292 frac1_hi = b.parts.frac_hi;
293 frac1_lo = b.parts.frac_lo;
294 exp1 = b.parts.exp;
295 frac2_hi = a.parts.frac_hi;
296 frac2_lo = a.parts.frac_lo;
297 exp2 = a.parts.exp;
298 expdiff *= -1;
299 } else {
300 if (is_float128_nan(a)) {
301 /* TODO: fix SigNaN */
302 if (is_float128_signan(a) || is_float128_signan(b)) {
303 }
304 return a;
305 }
306
307 /* a is infinity and b not */
308 if (a.parts.exp == FLOAT128_MAX_EXPONENT) {
309 return a;
310 }
311
312 frac1_hi = a.parts.frac_hi;
313 frac1_lo = a.parts.frac_lo;
314 exp1 = a.parts.exp;
315 frac2_hi = b.parts.frac_hi;
316 frac2_lo = b.parts.frac_lo;
317 exp2 = b.parts.exp;
318 }
319
320 if (exp1 == 0) {
321 /* both are denormalized */
322 add128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &frac1_hi, &frac1_lo);
323
324 and128(frac1_hi, frac1_lo,
325 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
326 &tmp_hi, &tmp_lo);
327 if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
328 /* result is not denormalized */
329 a.parts.exp = 1;
330 }
331
332 a.parts.frac_hi = frac1_hi;
333 a.parts.frac_lo = frac1_lo;
334 return a;
335 }
336
337 /* add hidden bit - frac1 is sure not denormalized */
338 or128(frac1_hi, frac1_lo,
339 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
340 &frac1_hi, &frac1_lo);
341
342 /* second operand ... */
343 if (exp2 == 0) {
344 /* ... is denormalized */
345 --expdiff;
346 } else {
347 /* is not denormalized */
348 or128(frac2_hi, frac2_lo,
349 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
350 &frac2_hi, &frac2_lo);
351 }
352
353 /* create some space for rounding */
354 lshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
355 lshift128(frac2_hi, frac2_lo, 6, &frac2_hi, &frac2_lo);
356
357 if (expdiff < (FLOAT128_FRACTION_SIZE + 2)) {
358 rshift128(frac2_hi, frac2_lo, expdiff, &frac2_hi, &frac2_lo);
359 add128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &frac1_hi, &frac1_lo);
360 } else {
361 a.parts.exp = exp1;
362
363 rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
364 not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
365 &tmp_hi, &tmp_lo);
366 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
367
368 a.parts.frac_hi = tmp_hi;
369 a.parts.frac_lo = tmp_lo;
370 return a;
371 }
372
373 lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7,
374 &tmp_hi, &tmp_lo);
375 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
376 if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
377 ++exp1;
378 rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
379 }
380
381 /* rounding - if first bit after fraction is set then round up */
382 add128(frac1_hi, frac1_lo, 0x0ll, 0x1ll << 5, &frac1_hi, &frac1_lo);
383
384 lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7,
385 &tmp_hi, &tmp_lo);
386 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
387 if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
388 /* rounding overflow */
389 ++exp1;
390 rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
391 }
392
393 if ((exp1 == FLOAT128_MAX_EXPONENT ) || (exp2 > exp1)) {
394 /* overflow - set infinity as result */
395 a.parts.exp = FLOAT64_MAX_EXPONENT;
396 a.parts.frac_hi = 0;
397 a.parts.frac_lo = 0;
398 return a;
399 }
400
401 a.parts.exp = exp1;
402
403 /* Clear hidden bit and shift */
404 rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
405 not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
406 &tmp_hi, &tmp_lo);
407 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
408
409 a.parts.frac_hi = tmp_hi;
410 a.parts.frac_lo = tmp_lo;
411
412 return a;
413}
414
415/** @}
416 */
Note: See TracBrowser for help on using the repository browser.