source: mainline/uspace/lib/softfloat/generic/sub.c@ 88d5c1e

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

softfloat redesign

  • avoid hardwired type sizes, actual sizes are determined at compile-time
  • add basic support for x87 extended-precision data types (stored as 96bit long double)
  • a lot of coding style changes (removal of CamelCase, etc.)
  • Property mode set to 100644
File size: 10.5 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 Substraction functions.
[846848a6]34 */
35
[750636a]36#include <sftypes.h>
37#include <sub.h>
38#include <comparison.h>
[c67aff2]39#include <common.h>
[12c6f2d]40
[88d5c1e]41/** Subtract two single-precision floats with the same sign.
[c67aff2]42 *
43 * @param a First input operand.
44 * @param b Second input operand.
[88d5c1e]45 *
[c67aff2]46 * @return Result of substraction.
[88d5c1e]47 *
[12c6f2d]48 */
[88d5c1e]49float32 sub_float32(float32 a, float32 b)
[12c6f2d]50{
51 int expdiff;
[aa59fa0]52 uint32_t exp1, exp2, frac1, frac2;
[12c6f2d]53 float32 result;
[88d5c1e]54
55 result.bin = 0;
[12c6f2d]56
[a96c570]57 expdiff = a.parts.exp - b.parts.exp;
[88d5c1e]58 if ((expdiff < 0 ) || ((expdiff == 0) &&
59 (a.parts.fraction < b.parts.fraction))) {
60 if (is_float32_nan(b)) {
61 if (is_float32_signan(b)) {
62 // TODO: fix SigNaN
[c67aff2]63 }
[88d5c1e]64
[12c6f2d]65 return b;
[c67aff2]66 }
[12c6f2d]67
[88d5c1e]68 if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
69 /* num -(+-inf) = -+inf */
70 b.parts.sign = !b.parts.sign;
[12c6f2d]71 return b;
72 }
73
[88d5c1e]74 result.parts.sign = !a.parts.sign;
[12c6f2d]75
[1266543]76 frac1 = b.parts.fraction;
[a96c570]77 exp1 = b.parts.exp;
[1266543]78 frac2 = a.parts.fraction;
[a96c570]79 exp2 = a.parts.exp;
80 expdiff *= -1;
[12c6f2d]81 } else {
[88d5c1e]82 if (is_float32_nan(a)) {
83 if ((is_float32_signan(a)) || (is_float32_signan(b))) {
84 // TODO: fix SigNaN
[c67aff2]85 }
[88d5c1e]86
[12c6f2d]87 return a;
[c67aff2]88 }
[12c6f2d]89
[88d5c1e]90 if (a.parts.exp == FLOAT32_MAX_EXPONENT) {
[a96c570]91 if (b.parts.exp == FLOAT32_MAX_EXPONENT) {
[12c6f2d]92 /* inf - inf => nan */
[88d5c1e]93 // TODO: fix exception
94 result.bin = FLOAT32_NAN;
[12c6f2d]95 return result;
[c67aff2]96 }
[88d5c1e]97
[12c6f2d]98 return a;
99 }
100
101 result.parts.sign = a.parts.sign;
102
[1266543]103 frac1 = a.parts.fraction;
[a96c570]104 exp1 = a.parts.exp;
[1266543]105 frac2 = b.parts.fraction;
[88d5c1e]106 exp2 = b.parts.exp;
[c67aff2]107 }
[12c6f2d]108
[a96c570]109 if (exp1 == 0) {
[1266543]110 /* both are denormalized */
[c67aff2]111 result.parts.fraction = frac1 - frac2;
[1266543]112 if (result.parts.fraction > frac1) {
[88d5c1e]113 // TODO: underflow exception
[12c6f2d]114 return result;
[c67aff2]115 }
[88d5c1e]116
[a96c570]117 result.parts.exp = 0;
[12c6f2d]118 return result;
[c67aff2]119 }
[88d5c1e]120
[a96c570]121 /* add hidden bit */
[88d5c1e]122 frac1 |= FLOAT32_HIDDEN_BIT_MASK;
[a96c570]123
124 if (exp2 == 0) {
125 /* denormalized */
[88d5c1e]126 --expdiff;
[a96c570]127 } else {
128 /* normalized */
[1266543]129 frac2 |= FLOAT32_HIDDEN_BIT_MASK;
[c67aff2]130 }
[a96c570]131
132 /* create some space for rounding */
[1266543]133 frac1 <<= 6;
134 frac2 <<= 6;
[a96c570]135
[88d5c1e]136 if (expdiff > FLOAT32_FRACTION_SIZE + 1)
[c67aff2]137 goto done;
[12c6f2d]138
[1266543]139 frac1 = frac1 - (frac2 >> expdiff);
[88d5c1e]140
[a96c570]141done:
[1266543]142 /* TODO: find first nonzero digit and shift result and detect possibly underflow */
143 while ((exp1 > 0) && (!(frac1 & (FLOAT32_HIDDEN_BIT_MASK << 6 )))) {
[a96c570]144 --exp1;
[1266543]145 frac1 <<= 1;
[c67aff2]146 /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
147 }
[12c6f2d]148
[1266543]149 /* rounding - if first bit after fraction is set then round up */
150 frac1 += 0x20;
[88d5c1e]151
[1266543]152 if (frac1 & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
[a96c570]153 ++exp1;
[1266543]154 frac1 >>= 1;
[c67aff2]155 }
[a96c570]156
[c67aff2]157 /* Clear hidden bit and shift */
[88d5c1e]158 result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
[a96c570]159 result.parts.exp = exp1;
160
161 return result;
162}
163
[88d5c1e]164/** Subtract two double-precision floats with the same sign.
[c67aff2]165 *
166 * @param a First input operand.
167 * @param b Second input operand.
[88d5c1e]168 *
[c67aff2]169 * @return Result of substraction.
[88d5c1e]170 *
[a96c570]171 */
[88d5c1e]172float64 sub_float64(float64 a, float64 b)
[a96c570]173{
174 int expdiff;
[aa59fa0]175 uint32_t exp1, exp2;
176 uint64_t frac1, frac2;
[a96c570]177 float64 result;
[88d5c1e]178
179 result.bin = 0;
[a96c570]180
181 expdiff = a.parts.exp - b.parts.exp;
[88d5c1e]182 if ((expdiff < 0 ) ||
183 ((expdiff == 0) && (a.parts.fraction < b.parts.fraction))) {
184 if (is_float64_nan(b)) {
185 if (is_float64_signan(b)) {
186 // TODO: fix SigNaN
[c67aff2]187 }
[88d5c1e]188
[a96c570]189 return b;
[c67aff2]190 }
[a96c570]191
[88d5c1e]192 if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
193 /* num -(+-inf) = -+inf */
194 b.parts.sign = !b.parts.sign;
[a96c570]195 return b;
196 }
197
[88d5c1e]198 result.parts.sign = !a.parts.sign;
[a96c570]199
[1266543]200 frac1 = b.parts.fraction;
[a96c570]201 exp1 = b.parts.exp;
[1266543]202 frac2 = a.parts.fraction;
[a96c570]203 exp2 = a.parts.exp;
204 expdiff *= -1;
205 } else {
[88d5c1e]206 if (is_float64_nan(a)) {
207 if (is_float64_signan(a) || is_float64_signan(b)) {
208 // TODO: fix SigNaN
[c67aff2]209 }
[88d5c1e]210
[a96c570]211 return a;
[c67aff2]212 }
[a96c570]213
[88d5c1e]214 if (a.parts.exp == FLOAT64_MAX_EXPONENT) {
[a96c570]215 if (b.parts.exp == FLOAT64_MAX_EXPONENT) {
216 /* inf - inf => nan */
[88d5c1e]217 // TODO: fix exception
218 result.bin = FLOAT64_NAN;
[a96c570]219 return result;
[c67aff2]220 }
[88d5c1e]221
[a96c570]222 return a;
223 }
224
225 result.parts.sign = a.parts.sign;
226
[1266543]227 frac1 = a.parts.fraction;
[a96c570]228 exp1 = a.parts.exp;
[1266543]229 frac2 = b.parts.fraction;
[88d5c1e]230 exp2 = b.parts.exp;
[c67aff2]231 }
[12c6f2d]232
[a96c570]233 if (exp1 == 0) {
[1266543]234 /* both are denormalized */
235 result.parts.fraction = frac1 - frac2;
236 if (result.parts.fraction > frac1) {
[88d5c1e]237 // TODO: underflow exception
[a96c570]238 return result;
[c67aff2]239 }
[88d5c1e]240
[a96c570]241 result.parts.exp = 0;
242 return result;
[c67aff2]243 }
[88d5c1e]244
[a96c570]245 /* add hidden bit */
[88d5c1e]246 frac1 |= FLOAT64_HIDDEN_BIT_MASK;
[12c6f2d]247
[a96c570]248 if (exp2 == 0) {
249 /* denormalized */
[88d5c1e]250 --expdiff;
[12c6f2d]251 } else {
[a96c570]252 /* normalized */
[1266543]253 frac2 |= FLOAT64_HIDDEN_BIT_MASK;
[c67aff2]254 }
[12c6f2d]255
[a96c570]256 /* create some space for rounding */
[1266543]257 frac1 <<= 6;
258 frac2 <<= 6;
[a96c570]259
[88d5c1e]260 if (expdiff > FLOAT64_FRACTION_SIZE + 1)
[c67aff2]261 goto done;
[12c6f2d]262
[1266543]263 frac1 = frac1 - (frac2 >> expdiff);
[88d5c1e]264
[12c6f2d]265done:
[1266543]266 /* TODO: find first nonzero digit and shift result and detect possibly underflow */
267 while ((exp1 > 0) && (!(frac1 & (FLOAT64_HIDDEN_BIT_MASK << 6 )))) {
[a96c570]268 --exp1;
[1266543]269 frac1 <<= 1;
[c67aff2]270 /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
271 }
[12c6f2d]272
[1266543]273 /* rounding - if first bit after fraction is set then round up */
274 frac1 += 0x20;
[88d5c1e]275
[1266543]276 if (frac1 & (FLOAT64_HIDDEN_BIT_MASK << 7)) {
[12c6f2d]277 ++exp1;
[1266543]278 frac1 >>= 1;
[c67aff2]279 }
[12c6f2d]280
[c67aff2]281 /* Clear hidden bit and shift */
[88d5c1e]282 result.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK));
[12c6f2d]283 result.parts.exp = exp1;
284
285 return result;
[a96c570]286}
[12c6f2d]287
[88d5c1e]288/** Subtract two quadruple-precision floats with the same sign.
[c67aff2]289 *
290 * @param a First input operand.
291 * @param b Second input operand.
[88d5c1e]292 *
[c67aff2]293 * @return Result of substraction.
[88d5c1e]294 *
[c67aff2]295 */
[88d5c1e]296float128 sub_float128(float128 a, float128 b)
[c67aff2]297{
298 int expdiff;
299 uint32_t exp1, exp2;
300 uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo;
301 float128 result;
[88d5c1e]302
303 result.bin.hi = 0;
304 result.bin.lo = 0;
305
[c67aff2]306 expdiff = a.parts.exp - b.parts.exp;
307 if ((expdiff < 0 ) || ((expdiff == 0) &&
308 lt128(a.parts.frac_hi, a.parts.frac_lo, b.parts.frac_hi, b.parts.frac_lo))) {
[88d5c1e]309 if (is_float128_nan(b)) {
310 if (is_float128_signan(b)) {
311 // TODO: fix SigNaN
[c67aff2]312 }
[88d5c1e]313
[c67aff2]314 return b;
315 }
[88d5c1e]316
[c67aff2]317 if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
[88d5c1e]318 /* num -(+-inf) = -+inf */
319 b.parts.sign = !b.parts.sign;
[c67aff2]320 return b;
321 }
[88d5c1e]322
[c67aff2]323 result.parts.sign = !a.parts.sign;
[88d5c1e]324
[c67aff2]325 frac1_hi = b.parts.frac_hi;
326 frac1_lo = b.parts.frac_lo;
327 exp1 = b.parts.exp;
328 frac2_hi = a.parts.frac_hi;
329 frac2_lo = a.parts.frac_lo;
330 exp2 = a.parts.exp;
331 expdiff *= -1;
332 } else {
[88d5c1e]333 if (is_float128_nan(a)) {
334 if (is_float128_signan(a) || is_float128_signan(b)) {
335 // TODO: fix SigNaN
[c67aff2]336 }
[88d5c1e]337
[c67aff2]338 return a;
339 }
[88d5c1e]340
[c67aff2]341 if (a.parts.exp == FLOAT128_MAX_EXPONENT) {
342 if (b.parts.exp == FLOAT128_MAX_EXPONENT) {
343 /* inf - inf => nan */
[88d5c1e]344 // TODO: fix exception
345 result.bin.hi = FLOAT128_NAN_HI;
346 result.bin.lo = FLOAT128_NAN_LO;
[c67aff2]347 return result;
348 }
349 return a;
350 }
[88d5c1e]351
[c67aff2]352 result.parts.sign = a.parts.sign;
[88d5c1e]353
[c67aff2]354 frac1_hi = a.parts.frac_hi;
355 frac1_lo = a.parts.frac_lo;
356 exp1 = a.parts.exp;
357 frac2_hi = b.parts.frac_hi;
358 frac2_lo = b.parts.frac_lo;
359 exp2 = b.parts.exp;
360 }
[88d5c1e]361
[c67aff2]362 if (exp1 == 0) {
363 /* both are denormalized */
364 sub128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &tmp_hi, &tmp_lo);
365 result.parts.frac_hi = tmp_hi;
366 result.parts.frac_lo = tmp_lo;
367 if (lt128(frac1_hi, frac1_lo, result.parts.frac_hi, result.parts.frac_lo)) {
[88d5c1e]368 // TODO: underflow exception
[c67aff2]369 return result;
370 }
[88d5c1e]371
[c67aff2]372 result.parts.exp = 0;
373 return result;
374 }
[88d5c1e]375
[c67aff2]376 /* add hidden bit */
377 or128(frac1_hi, frac1_lo,
378 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
379 &frac1_hi, &frac1_lo);
[88d5c1e]380
[c67aff2]381 if (exp2 == 0) {
382 /* denormalized */
383 --expdiff;
384 } else {
385 /* normalized */
386 or128(frac2_hi, frac2_lo,
387 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
388 &frac2_hi, &frac2_lo);
389 }
[88d5c1e]390
[c67aff2]391 /* create some space for rounding */
392 lshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
393 lshift128(frac2_hi, frac2_lo, 6, &frac2_hi, &frac2_lo);
[88d5c1e]394
395 if (expdiff > FLOAT128_FRACTION_SIZE + 1)
[c67aff2]396 goto done;
[88d5c1e]397
[c67aff2]398 rshift128(frac2_hi, frac2_lo, expdiff, &tmp_hi, &tmp_lo);
399 sub128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo);
[88d5c1e]400
[c67aff2]401done:
402 /* TODO: find first nonzero digit and shift result and detect possibly underflow */
403 lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6,
404 &tmp_hi, &tmp_lo);
405 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
406 while ((exp1 > 0) && (!lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo))) {
407 --exp1;
408 lshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
409 /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */
[88d5c1e]410
[c67aff2]411 lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6,
412 &tmp_hi, &tmp_lo);
413 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
414 }
[88d5c1e]415
[c67aff2]416 /* rounding - if first bit after fraction is set then round up */
417 add128(frac1_hi, frac1_lo, 0x0ll, 0x20ll, &frac1_hi, &frac1_lo);
[88d5c1e]418
[c67aff2]419 lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7,
420 &tmp_hi, &tmp_lo);
421 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
422 if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) {
423 ++exp1;
424 rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo);
425 }
[88d5c1e]426
[c67aff2]427 /* Clear hidden bit and shift */
428 rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo);
429 not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
430 &tmp_hi, &tmp_lo);
431 and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo);
432 result.parts.frac_hi = tmp_hi;
433 result.parts.frac_lo = tmp_lo;
[88d5c1e]434
[c67aff2]435 result.parts.exp = exp1;
[88d5c1e]436
[c67aff2]437 return result;
438}
439
[750636a]440/** @}
[846848a6]441 */
Note: See TracBrowser for help on using the repository browser.