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

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

softfloat: move sources to a single directory
implement more ARM EABI bindings
  • Property mode set to 100644
File size: 9.9 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 "sftypes.h"
37#include "add.h"
38#include "comparison.h"
39#include "common.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;
[12c6f2d]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 }
[12c6f2d]61
[4a5abddd]62 if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
[12c6f2d]63 return b;
64 }
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 }
[12c6f2d]78
[4a5abddd]79 if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
[12c6f2d]80 return a;
81 }
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 }
[12c6f2d]88
89 if (exp1 == 0) {
90 /* both are denormalized */
[1266543]91 frac1 += frac2;
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 }
[12c6f2d]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 */
[1266543]107 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
[c67aff2]108 }
[12c6f2d]109
110 /* create some space for rounding */
[1266543]111 frac1 <<= 6;
112 frac2 <<= 6;
[12c6f2d]113
[1266543]114 if (expdiff < (FLOAT32_FRACTION_SIZE + 2) ) {
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 }
[12c6f2d]122
[1266543]123 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7) ) {
[12c6f2d]124 ++exp1;
[1266543]125 frac1 >>= 1;
[c67aff2]126 }
[12c6f2d]127
[1266543]128 /* rounding - if first bit after fraction is set then round up */
129 frac1 += (0x1 << 5);
[12c6f2d]130
[1266543]131 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
132 /* rounding overflow */
[12c6f2d]133 ++exp1;
[1266543]134 frac1 >>= 1;
[c67aff2]135 }
[e6a40ac]136
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 }
[1266543]143
[12c6f2d]144 a.parts.exp = exp1;
145
[750636a]146 /* Clear hidden bit and shift */
[c67aff2]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;
[12c6f2d]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 }
[12c6f2d]172
[88d5c1e]173 /* b is infinity and a not */
174 if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
[12c6f2d]175 return b;
176 }
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 }
[12c6f2d]190
191 /* a is infinity and b not */
[c67aff2]192 if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
[12c6f2d]193 return a;
194 }
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 }
[12c6f2d]201
202 if (exp1 == 0) {
203 /* both are denormalized */
[1266543]204 frac1 += frac2;
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 }
[12c6f2d]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 */
219 --expdiff;
220 } else {
221 /* is not denormalized */
[1266543]222 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
[c67aff2]223 }
[12c6f2d]224
225 /* create some space for rounding */
[1266543]226 frac1 <<= 6;
227 frac2 <<= 6;
[12c6f2d]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 }
[12c6f2d]237
[c67aff2]238 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
[12c6f2d]239 ++exp1;
[1266543]240 frac1 >>= 1;
[c67aff2]241 }
[12c6f2d]242
[1266543]243 /* rounding - if first bit after fraction is set then round up */
244 frac1 += (0x1 << 5);
[12c6f2d]245
[1266543]246 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
247 /* rounding overflow */
[12c6f2d]248 ++exp1;
[1266543]249 frac1 >>= 1;
[c67aff2]250 }
[12c6f2d]251
[d3ca210]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 }
[1266543]258
[12c6f2d]259 a.parts.exp = exp1;
[750636a]260 /* Clear hidden bit and shift */
[c67aff2]261 a.parts.fraction = ((frac1 >> 6 ) & (~FLOAT64_HIDDEN_BIT_MASK));
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,
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;
[d3ca210]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
[231a60a]415/** @}
[846848a6]416 */
Note: See TracBrowser for help on using the repository browser.