source: mainline/uspace/lib/softfloat/add.c@ 22fb7ab

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since 22fb7ab was c0c38c7c, checked in by Martin Decky <martin@…>, 10 years ago

software floating point overhaul
use proper type mapping
fix cosine calculation
  • Property mode set to 100644
File size: 12.2 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 "add.h"
37#include "comparison.h"
38#include "common.h"
39#include "sub.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#ifdef float32_t
416
417float32_t __addsf3(float32_t a, float32_t b)
418{
419 float32_u ua;
420 ua.val = a;
421
422 float32_u ub;
423 ub.val = b;
424
425 float32_u res;
426
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);
437
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;
445
446 float32_u ub;
447 ub.val = b;
448
449 float32_u res;
450
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);
461
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;
473
474 float64_u ub;
475 ub.val = b;
476
477 float64_u res;
478
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);
489
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;
497
498 float64_u ub;
499 ub.val = b;
500
501 float64_u res;
502
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);
513
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;
525
526 float128_u ub;
527 ub.val = b;
528
529 float128_u res;
530
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);
541
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
552/** @}
553 */
Note: See TracBrowser for help on using the repository browser.