source: mainline/uspace/lib/softfloat/add.c@ ab6edb6

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since ab6edb6 was 1433ecda, checked in by Jiri Svoboda <jiri@…>, 7 years ago

Fix cstyle: make ccheck-fix and commit only files where all the changes are good.
  • Property mode set to 100644
File size: 12.2 KB
RevLine 
[12c6f2d]1/*
[df4ed85]2 * Copyright (c) 2005 Josef Cejka
[c67aff2]3 * Copyright (c) 2011 Petr Koupy
[12c6f2d]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
[750636a]30/** @addtogroup softfloat
[846848a6]31 * @{
32 */
[c67aff2]33/** @file Addition functions.
[846848a6]34 */
35
[2416085]36#include "add.h"
37#include "comparison.h"
38#include "common.h"
[c0c38c7c]39#include "sub.h"
[12c6f2d]40
[88d5c1e]41/** Add two single-precision floats with the same sign.
[c67aff2]42 *
43 * @param a First input operand.
44 * @param b Second input operand.
45 * @return Result of addition.
[12c6f2d]46 */
[88d5c1e]47float32 add_float32(float32 a, float32 b)
[12c6f2d]48{
49 int expdiff;
[c67aff2]50 uint32_t exp1, exp2, frac1, frac2;
[a35b458]51
[4a5abddd]52 expdiff = a.parts.exp - b.parts.exp;
53 if (expdiff < 0) {
[88d5c1e]54 if (is_float32_nan(b)) {
[1266543]55 /* TODO: fix SigNaN */
[88d5c1e]56 if (is_float32_signan(b)) {
[c67aff2]57 }
[12c6f2d]58
59 return b;
[c67aff2]60 }
[a35b458]61
[1b20da0]62 if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
[12c6f2d]63 return b;
64 }
[a35b458]65
[1266543]66 frac1 = b.parts.fraction;
[4a5abddd]67 exp1 = b.parts.exp;
[1266543]68 frac2 = a.parts.fraction;
[4a5abddd]69 exp2 = a.parts.exp;
70 expdiff *= -1;
[12c6f2d]71 } else {
[88d5c1e]72 if ((is_float32_nan(a)) || (is_float32_nan(b))) {
[1266543]73 /* TODO: fix SigNaN */
[88d5c1e]74 if (is_float32_signan(a) || is_float32_signan(b)) {
[c67aff2]75 }
[88d5c1e]76 return (is_float32_nan(a) ? a : b);
[c67aff2]77 }
[a35b458]78
[1b20da0]79 if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
[12c6f2d]80 return a;
81 }
[a35b458]82
[1266543]83 frac1 = a.parts.fraction;
[4a5abddd]84 exp1 = a.parts.exp;
[1266543]85 frac2 = b.parts.fraction;
[4a5abddd]86 exp2 = b.parts.exp;
[c67aff2]87 }
[a35b458]88
[12c6f2d]89 if (exp1 == 0) {
90 /* both are denormalized */
[1266543]91 frac1 += frac2;
[faa45c17]92 if (frac1 & FLOAT32_HIDDEN_BIT_MASK) {
[12c6f2d]93 /* result is not denormalized */
94 a.parts.exp = 1;
[c67aff2]95 }
[1266543]96 a.parts.fraction = frac1;
[12c6f2d]97 return a;
[c67aff2]98 }
[a35b458]99
[1266543]100 frac1 |= FLOAT32_HIDDEN_BIT_MASK; /* add hidden bit */
[12c6f2d]101
102 if (exp2 == 0) {
103 /* second operand is denormalized */
[e6a40ac]104 --expdiff;
[12c6f2d]105 } else {
106 /* add hidden bit to second operand */
[1b20da0]107 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
[c67aff2]108 }
[a35b458]109
[12c6f2d]110 /* create some space for rounding */
[1266543]111 frac1 <<= 6;
112 frac2 <<= 6;
[a35b458]113
[faa45c17]114 if (expdiff < (FLOAT32_FRACTION_SIZE + 2)) {
[1266543]115 frac2 >>= expdiff;
116 frac1 += frac2;
[d3ca210]117 } else {
118 a.parts.exp = exp1;
119 a.parts.fraction = (frac1 >> 6) & (~(FLOAT32_HIDDEN_BIT_MASK));
120 return a;
121 }
[a35b458]122
[faa45c17]123 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
[12c6f2d]124 ++exp1;
[1266543]125 frac1 >>= 1;
[c67aff2]126 }
[a35b458]127
[1266543]128 /* rounding - if first bit after fraction is set then round up */
129 frac1 += (0x1 << 5);
[a35b458]130
[1b20da0]131 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
[1266543]132 /* rounding overflow */
[12c6f2d]133 ++exp1;
[1266543]134 frac1 >>= 1;
[c67aff2]135 }
[a35b458]136
[faa45c17]137 if ((exp1 == FLOAT32_MAX_EXPONENT) || (exp2 > exp1)) {
[c67aff2]138 /* overflow - set infinity as result */
139 a.parts.exp = FLOAT32_MAX_EXPONENT;
140 a.parts.fraction = 0;
141 return a;
142 }
[a35b458]143
[12c6f2d]144 a.parts.exp = exp1;
[a35b458]145
[750636a]146 /* Clear hidden bit and shift */
[1b20da0]147 a.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
[12c6f2d]148 return a;
149}
150
[88d5c1e]151/** Add two double-precision floats with the same sign.
[c67aff2]152 *
153 * @param a First input operand.
154 * @param b Second input operand.
155 * @return Result of addition.
[12c6f2d]156 */
[88d5c1e]157float64 add_float64(float64 a, float64 b)
[12c6f2d]158{
159 int expdiff;
[aa59fa0]160 uint32_t exp1, exp2;
161 uint64_t frac1, frac2;
[a35b458]162
[c67aff2]163 expdiff = ((int) a.parts.exp) - b.parts.exp;
[4a5abddd]164 if (expdiff < 0) {
[88d5c1e]165 if (is_float64_nan(b)) {
[1266543]166 /* TODO: fix SigNaN */
[88d5c1e]167 if (is_float64_signan(b)) {
[c67aff2]168 }
[12c6f2d]169
170 return b;
[c67aff2]171 }
[a35b458]172
[88d5c1e]173 /* b is infinity and a not */
174 if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
[12c6f2d]175 return b;
176 }
[a35b458]177
[1266543]178 frac1 = b.parts.fraction;
[4a5abddd]179 exp1 = b.parts.exp;
[1266543]180 frac2 = a.parts.fraction;
[4a5abddd]181 exp2 = a.parts.exp;
182 expdiff *= -1;
[12c6f2d]183 } else {
[88d5c1e]184 if (is_float64_nan(a)) {
[1266543]185 /* TODO: fix SigNaN */
[88d5c1e]186 if (is_float64_signan(a) || is_float64_signan(b)) {
[c67aff2]187 }
[12c6f2d]188 return a;
[c67aff2]189 }
[a35b458]190
[12c6f2d]191 /* a is infinity and b not */
[1b20da0]192 if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
[12c6f2d]193 return a;
194 }
[a35b458]195
[1266543]196 frac1 = a.parts.fraction;
[12c6f2d]197 exp1 = a.parts.exp;
[1266543]198 frac2 = b.parts.fraction;
[12c6f2d]199 exp2 = b.parts.exp;
[c67aff2]200 }
[a35b458]201
[12c6f2d]202 if (exp1 == 0) {
203 /* both are denormalized */
[1266543]204 frac1 += frac2;
[1b20da0]205 if (frac1 & FLOAT64_HIDDEN_BIT_MASK) {
[12c6f2d]206 /* result is not denormalized */
207 a.parts.exp = 1;
[c67aff2]208 }
[1266543]209 a.parts.fraction = frac1;
[12c6f2d]210 return a;
[c67aff2]211 }
[a35b458]212
[1266543]213 /* add hidden bit - frac1 is sure not denormalized */
214 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
[12c6f2d]215
216 /* second operand ... */
217 if (exp2 == 0) {
218 /* ... is denormalized */
[1b20da0]219 --expdiff;
[12c6f2d]220 } else {
221 /* is not denormalized */
[1266543]222 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
[c67aff2]223 }
[a35b458]224
[12c6f2d]225 /* create some space for rounding */
[1266543]226 frac1 <<= 6;
227 frac2 <<= 6;
[a35b458]228
[c67aff2]229 if (expdiff < (FLOAT64_FRACTION_SIZE + 2)) {
[1266543]230 frac2 >>= expdiff;
231 frac1 += frac2;
[d3ca210]232 } else {
233 a.parts.exp = exp1;
234 a.parts.fraction = (frac1 >> 6) & (~(FLOAT64_HIDDEN_BIT_MASK));
235 return a;
236 }
[a35b458]237
[c67aff2]238 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
[12c6f2d]239 ++exp1;
[1266543]240 frac1 >>= 1;
[c67aff2]241 }
[a35b458]242
[1266543]243 /* rounding - if first bit after fraction is set then round up */
[1b20da0]244 frac1 += (0x1 << 5);
[a35b458]245
[1b20da0]246 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
[1266543]247 /* rounding overflow */
[12c6f2d]248 ++exp1;
[1266543]249 frac1 >>= 1;
[c67aff2]250 }
[a35b458]251
[faa45c17]252 if ((exp1 == FLOAT64_MAX_EXPONENT) || (exp2 > exp1)) {
[c67aff2]253 /* overflow - set infinity as result */
254 a.parts.exp = FLOAT64_MAX_EXPONENT;
255 a.parts.fraction = 0;
256 return a;
257 }
[a35b458]258
[12c6f2d]259 a.parts.exp = exp1;
[750636a]260 /* Clear hidden bit and shift */
[faa45c17]261 a.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK));
[c67aff2]262 return a;
263}
264
[88d5c1e]265/** Add two quadruple-precision floats with the same sign.
[c67aff2]266 *
267 * @param a First input operand.
268 * @param b Second input operand.
269 * @return Result of addition.
270 */
[88d5c1e]271float128 add_float128(float128 a, float128 b)
[c67aff2]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) {
[88d5c1e]279 if (is_float128_nan(b)) {
[c67aff2]280 /* TODO: fix SigNaN */
[88d5c1e]281 if (is_float128_signan(b)) {
[c67aff2]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 {
[88d5c1e]300 if (is_float128_nan(a)) {
[c67aff2]301 /* TODO: fix SigNaN */
[88d5c1e]302 if (is_float128_signan(a) || is_float128_signan(b)) {
[c67aff2]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,
[1433ecda]385 &tmp_hi, &tmp_lo);
[c67aff2]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
[1433ecda]393 if ((exp1 == FLOAT128_MAX_EXPONENT) || (exp2 > exp1)) {
[c67aff2]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;
[a35b458]402
[c67aff2]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
[12c6f2d]412 return a;
413}
414
[c0c38c7c]415#ifdef float32_t
416
417float32_t __addsf3(float32_t a, float32_t b)
418{
419 float32_u ua;
420 ua.val = a;
[a35b458]421
[c0c38c7c]422 float32_u ub;
423 ub.val = b;
[a35b458]424
[c0c38c7c]425 float32_u res;
[a35b458]426
[c0c38c7c]427 if (ua.data.parts.sign != ub.data.parts.sign) {
428 if (ua.data.parts.sign) {
429 ua.data.parts.sign = 0;
430 res.data = sub_float32(ub.data, ua.data);
431 } else {
432 ub.data.parts.sign = 0;
433 res.data = sub_float32(ua.data, ub.data);
434 }
435 } else
436 res.data = add_float32(ua.data, ub.data);
[a35b458]437
[c0c38c7c]438 return res.val;
439}
440
441float32_t __aeabi_fadd(float32_t a, float32_t b)
442{
443 float32_u ua;
444 ua.val = a;
[a35b458]445
[c0c38c7c]446 float32_u ub;
447 ub.val = b;
[a35b458]448
[c0c38c7c]449 float32_u res;
[a35b458]450
[c0c38c7c]451 if (ua.data.parts.sign != ub.data.parts.sign) {
452 if (ua.data.parts.sign) {
453 ua.data.parts.sign = 0;
454 res.data = sub_float32(ub.data, ua.data);
455 } else {
456 ub.data.parts.sign = 0;
457 res.data = sub_float32(ua.data, ub.data);
458 }
459 } else
460 res.data = add_float32(ua.data, ub.data);
[a35b458]461
[c0c38c7c]462 return res.val;
463}
464
465#endif
466
467#ifdef float64_t
468
469float64_t __adddf3(float64_t a, float64_t b)
470{
471 float64_u ua;
472 ua.val = a;
[a35b458]473
[c0c38c7c]474 float64_u ub;
475 ub.val = b;
[a35b458]476
[c0c38c7c]477 float64_u res;
[a35b458]478
[c0c38c7c]479 if (ua.data.parts.sign != ub.data.parts.sign) {
480 if (ua.data.parts.sign) {
481 ua.data.parts.sign = 0;
482 res.data = sub_float64(ub.data, ua.data);
483 } else {
484 ub.data.parts.sign = 0;
485 res.data = sub_float64(ua.data, ub.data);
486 }
487 } else
488 res.data = add_float64(ua.data, ub.data);
[a35b458]489
[c0c38c7c]490 return res.val;
491}
492
493float64_t __aeabi_dadd(float64_t a, float64_t b)
494{
495 float64_u ua;
496 ua.val = a;
[a35b458]497
[c0c38c7c]498 float64_u ub;
499 ub.val = b;
[a35b458]500
[c0c38c7c]501 float64_u res;
[a35b458]502
[c0c38c7c]503 if (ua.data.parts.sign != ub.data.parts.sign) {
504 if (ua.data.parts.sign) {
505 ua.data.parts.sign = 0;
506 res.data = sub_float64(ub.data, ua.data);
507 } else {
508 ub.data.parts.sign = 0;
509 res.data = sub_float64(ua.data, ub.data);
510 }
511 } else
512 res.data = add_float64(ua.data, ub.data);
[a35b458]513
[c0c38c7c]514 return res.val;
515}
516
517#endif
518
519#ifdef float128_t
520
521float128_t __addtf3(float128_t a, float128_t b)
522{
523 float128_u ua;
524 ua.val = a;
[a35b458]525
[c0c38c7c]526 float128_u ub;
527 ub.val = b;
[a35b458]528
[c0c38c7c]529 float128_u res;
[a35b458]530
[c0c38c7c]531 if (ua.data.parts.sign != ub.data.parts.sign) {
532 if (ua.data.parts.sign) {
533 ua.data.parts.sign = 0;
534 res.data = sub_float128(ub.data, ua.data);
535 } else {
536 ub.data.parts.sign = 0;
537 res.data = sub_float128(ua.data, ub.data);
538 }
539 } else
540 res.data = add_float128(ua.data, ub.data);
[a35b458]541
[c0c38c7c]542 return res.val;
543}
544
545void _Qp_add(float128_t *c, float128_t *a, float128_t *b)
546{
547 *c = __addtf3(*a, *b);
548}
549
550#endif
551
[231a60a]552/** @}
[846848a6]553 */
Note: See TracBrowser for help on using the repository browser.