Changes in uspace/lib/softfloat/generic/mul.c [750636a:c67aff2] in mainline
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
uspace/lib/softfloat/generic/mul.c
r750636a rc67aff2 1 1 /* 2 2 * Copyright (c) 2005 Josef Cejka 3 * Copyright (c) 2011 Petr Koupy 3 4 * All rights reserved. 4 5 * … … 30 31 * @{ 31 32 */ 32 /** @file 33 /** @file Multiplication functions. 33 34 */ 34 35 … … 38 39 #include <common.h> 39 40 40 /** Multiply two 32 bit float numbers 41 * 41 /** 42 * Multiply two single-precision floats. 43 * 44 * @param a First input operand. 45 * @param b Second input operand. 46 * @return Result of multiplication. 42 47 */ 43 48 float32 mulFloat32(float32 a, float32 b) … … 49 54 result.parts.sign = a.parts.sign ^ b.parts.sign; 50 55 51 if (isFloat32NaN(a) || isFloat32NaN(b) 56 if (isFloat32NaN(a) || isFloat32NaN(b)) { 52 57 /* TODO: fix SigNaNs */ 53 58 if (isFloat32SigNaN(a)) { … … 55 60 result.parts.exp = a.parts.exp; 56 61 return result; 57 } ;62 } 58 63 if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */ 59 64 result.parts.fraction = b.parts.fraction; 60 65 result.parts.exp = b.parts.exp; 61 66 return result; 62 } ;67 } 63 68 /* set NaN as result */ 64 69 result.binary = FLOAT32_NAN; 65 70 return result; 66 } ;71 } 67 72 68 73 if (isFloat32Infinity(a)) { … … 98 103 result.parts.sign = a.parts.sign ^ b.parts.sign; 99 104 return result; 100 } ;105 } 101 106 102 107 if (exp < 0) { … … 106 111 result.parts.exp = 0x0; 107 112 return result; 108 } ;113 } 109 114 110 115 frac1 = a.parts.fraction; … … 113 118 } else { 114 119 ++exp; 115 } ;120 } 116 121 117 122 frac2 = b.parts.fraction; … … 121 126 } else { 122 127 ++exp; 123 } ;128 } 124 129 125 130 frac1 <<= 1; /* one bit space for rounding */ 126 131 127 132 frac1 = frac1 * frac2; 128 /* round and return */ 129 130 while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 131 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left) */133 134 /* round and return */ 135 while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= (1 << (FLOAT32_FRACTION_SIZE + 2)))) { 136 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left) */ 132 137 ++exp; 133 138 frac1 >>= 1; 134 } ;139 } 135 140 136 141 /* rounding */ … … 141 146 ++exp; 142 147 frac1 >>= 1; 143 } ;144 145 if (exp >= FLOAT32_MAX_EXPONENT 148 } 149 150 if (exp >= FLOAT32_MAX_EXPONENT) { 146 151 /* TODO: fix overflow */ 147 152 /* return infinity*/ … … 159 164 frac1 >>= 1; 160 165 ++exp; 161 } ;166 } 162 167 if (frac1 == 0) { 163 168 /* FIXME : underflow */ 164 result.parts.exp = 0;165 result.parts.fraction = 0;166 return result;167 } ;168 } ;169 result.parts.exp = 0; 170 result.parts.fraction = 0; 171 return result; 172 } 173 } 169 174 result.parts.exp = exp; 170 result.parts.fraction = frac1 & ( 175 result.parts.fraction = frac1 & ((1 << FLOAT32_FRACTION_SIZE) - 1); 171 176 172 177 return result; 173 174 178 } 175 179 176 /** Multiply two 64 bit float numbers 177 * 180 /** 181 * Multiply two double-precision floats. 182 * 183 * @param a First input operand. 184 * @param b Second input operand. 185 * @return Result of multiplication. 178 186 */ 179 187 float64 mulFloat64(float64 a, float64 b) … … 185 193 result.parts.sign = a.parts.sign ^ b.parts.sign; 186 194 187 if (isFloat64NaN(a) || isFloat64NaN(b) 195 if (isFloat64NaN(a) || isFloat64NaN(b)) { 188 196 /* TODO: fix SigNaNs */ 189 197 if (isFloat64SigNaN(a)) { … … 191 199 result.parts.exp = a.parts.exp; 192 200 return result; 193 } ;201 } 194 202 if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */ 195 203 result.parts.fraction = b.parts.fraction; 196 204 result.parts.exp = b.parts.exp; 197 205 return result; 198 } ;206 } 199 207 /* set NaN as result */ 200 208 result.binary = FLOAT64_NAN; 201 209 return result; 202 } ;210 } 203 211 204 212 if (isFloat64Infinity(a)) { … … 233 241 } else { 234 242 ++exp; 235 } ;243 } 236 244 237 245 frac2 = b.parts.fraction; … … 241 249 } else { 242 250 ++exp; 243 } ;251 } 244 252 245 253 frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1); 246 254 frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2); 247 255 248 mul64 integers(frac1, frac2, &frac1, &frac2);249 250 frac 2 |= (frac1!= 0);251 if (frac 2& (0x1ll << 62)) {252 frac 2<<= 1;256 mul64(frac1, frac2, &frac1, &frac2); 257 258 frac1 |= (frac2 != 0); 259 if (frac1 & (0x1ll << 62)) { 260 frac1 <<= 1; 253 261 exp--; 254 262 } 255 263 256 result = finishFloat64(exp, frac 2, result.parts.sign);264 result = finishFloat64(exp, frac1, result.parts.sign); 257 265 return result; 258 266 } 259 267 260 /** Multiply two 64 bit numbers and return result in two parts 261 * @param a first operand 262 * @param b second operand 263 * @param lo lower part from result 264 * @param hi higher part of result 265 */ 266 void mul64integers(uint64_t a,uint64_t b, uint64_t *lo, uint64_t *hi) 268 /** 269 * Multiply two quadruple-precision floats. 270 * 271 * @param a First input operand. 272 * @param b Second input operand. 273 * @return Result of multiplication. 274 */ 275 float128 mulFloat128(float128 a, float128 b) 267 276 { 268 uint64_t low, high, middle1, middle2; 269 uint32_t alow, blow; 270 271 alow = a & 0xFFFFFFFF; 272 blow = b & 0xFFFFFFFF; 273 274 a >>= 32; 275 b >>= 32; 276 277 low = ((uint64_t)alow) * blow; 278 middle1 = a * blow; 279 middle2 = alow * b; 280 high = a * b; 281 282 middle1 += middle2; 283 high += (((uint64_t)(middle1 < middle2)) << 32) + (middle1 >> 32); 284 middle1 <<= 32; 285 low += middle1; 286 high += (low < middle1); 287 *lo = low; 288 *hi = high; 289 290 return; 277 float128 result; 278 uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo; 279 int32_t exp; 280 281 result.parts.sign = a.parts.sign ^ b.parts.sign; 282 283 if (isFloat128NaN(a) || isFloat128NaN(b)) { 284 /* TODO: fix SigNaNs */ 285 if (isFloat128SigNaN(a)) { 286 result.parts.frac_hi = a.parts.frac_hi; 287 result.parts.frac_lo = a.parts.frac_lo; 288 result.parts.exp = a.parts.exp; 289 return result; 290 } 291 if (isFloat128SigNaN(b)) { /* TODO: fix SigNaN */ 292 result.parts.frac_hi = b.parts.frac_hi; 293 result.parts.frac_lo = b.parts.frac_lo; 294 result.parts.exp = b.parts.exp; 295 return result; 296 } 297 /* set NaN as result */ 298 result.binary.hi = FLOAT128_NAN_HI; 299 result.binary.lo = FLOAT128_NAN_LO; 300 return result; 301 } 302 303 if (isFloat128Infinity(a)) { 304 if (isFloat128Zero(b)) { 305 /* FIXME: zero * infinity */ 306 result.binary.hi = FLOAT128_NAN_HI; 307 result.binary.lo = FLOAT128_NAN_LO; 308 return result; 309 } 310 result.parts.frac_hi = a.parts.frac_hi; 311 result.parts.frac_lo = a.parts.frac_lo; 312 result.parts.exp = a.parts.exp; 313 return result; 314 } 315 316 if (isFloat128Infinity(b)) { 317 if (isFloat128Zero(a)) { 318 /* FIXME: zero * infinity */ 319 result.binary.hi = FLOAT128_NAN_HI; 320 result.binary.lo = FLOAT128_NAN_LO; 321 return result; 322 } 323 result.parts.frac_hi = b.parts.frac_hi; 324 result.parts.frac_lo = b.parts.frac_lo; 325 result.parts.exp = b.parts.exp; 326 return result; 327 } 328 329 /* exp is signed so we can easy detect underflow */ 330 exp = a.parts.exp + b.parts.exp - FLOAT128_BIAS - 1; 331 332 frac1_hi = a.parts.frac_hi; 333 frac1_lo = a.parts.frac_lo; 334 335 if (a.parts.exp > 0) { 336 or128(frac1_hi, frac1_lo, 337 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 338 &frac1_hi, &frac1_lo); 339 } else { 340 ++exp; 341 } 342 343 frac2_hi = b.parts.frac_hi; 344 frac2_lo = b.parts.frac_lo; 345 346 if (b.parts.exp > 0) { 347 or128(frac2_hi, frac2_lo, 348 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 349 &frac2_hi, &frac2_lo); 350 } else { 351 ++exp; 352 } 353 354 lshift128(frac2_hi, frac2_lo, 355 128 - FLOAT128_FRACTION_SIZE, &frac2_hi, &frac2_lo); 356 357 tmp_hi = frac1_hi; 358 tmp_lo = frac1_lo; 359 mul128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, 360 &frac1_hi, &frac1_lo, &frac2_hi, &frac2_lo); 361 add128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo); 362 frac2_hi |= (frac2_lo != 0x0ll); 363 364 if ((FLOAT128_HIDDEN_BIT_MASK_HI << 1) <= frac1_hi) { 365 frac2_hi >>= 1; 366 if (frac1_lo & 0x1ll) { 367 frac2_hi |= (0x1ull < 64); 368 } 369 rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo); 370 ++exp; 371 } 372 373 result = finishFloat128(exp, frac1_hi, frac1_lo, result.parts.sign, frac2_hi); 374 return result; 291 375 } 292 376
Note:
See TracChangeset
for help on using the changeset viewer.