Changes in uspace/lib/softfloat/generic/mul.c [c67aff2:750636a] in mainline
- File:
-
- 1 edited
-
uspace/lib/softfloat/generic/mul.c (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
uspace/lib/softfloat/generic/mul.c
rc67aff2 r750636a 1 1 /* 2 2 * Copyright (c) 2005 Josef Cejka 3 * Copyright (c) 2011 Petr Koupy4 3 * All rights reserved. 5 4 * … … 31 30 * @{ 32 31 */ 33 /** @file Multiplication functions.32 /** @file 34 33 */ 35 34 … … 39 38 #include <common.h> 40 39 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. 40 /** Multiply two 32 bit float numbers 41 * 47 42 */ 48 43 float32 mulFloat32(float32 a, float32 b) … … 54 49 result.parts.sign = a.parts.sign ^ b.parts.sign; 55 50 56 if (isFloat32NaN(a) || isFloat32NaN(b) ) {51 if (isFloat32NaN(a) || isFloat32NaN(b) ) { 57 52 /* TODO: fix SigNaNs */ 58 53 if (isFloat32SigNaN(a)) { … … 60 55 result.parts.exp = a.parts.exp; 61 56 return result; 62 } 57 }; 63 58 if (isFloat32SigNaN(b)) { /* TODO: fix SigNaN */ 64 59 result.parts.fraction = b.parts.fraction; 65 60 result.parts.exp = b.parts.exp; 66 61 return result; 67 } 62 }; 68 63 /* set NaN as result */ 69 64 result.binary = FLOAT32_NAN; 70 65 return result; 71 } 66 }; 72 67 73 68 if (isFloat32Infinity(a)) { … … 103 98 result.parts.sign = a.parts.sign ^ b.parts.sign; 104 99 return result; 105 } 100 }; 106 101 107 102 if (exp < 0) { … … 111 106 result.parts.exp = 0x0; 112 107 return result; 113 } 108 }; 114 109 115 110 frac1 = a.parts.fraction; … … 118 113 } else { 119 114 ++exp; 120 } 115 }; 121 116 122 117 frac2 = b.parts.fraction; … … 126 121 } else { 127 122 ++exp; 128 } 123 }; 129 124 130 125 frac1 <<= 1; /* one bit space for rounding */ 131 126 132 127 frac1 = frac1 * frac2; 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) 128 /* round and return */ 129 130 while ((exp < FLOAT32_MAX_EXPONENT) && (frac1 >= ( 1 << (FLOAT32_FRACTION_SIZE + 2)))) { 131 /* 23 bits of fraction + one more for hidden bit (all shifted 1 bit left)*/ 137 132 ++exp; 138 133 frac1 >>= 1; 139 } 134 }; 140 135 141 136 /* rounding */ … … 146 141 ++exp; 147 142 frac1 >>= 1; 148 } 149 150 if (exp >= FLOAT32_MAX_EXPONENT ) {143 }; 144 145 if (exp >= FLOAT32_MAX_EXPONENT ) { 151 146 /* TODO: fix overflow */ 152 147 /* return infinity*/ … … 164 159 frac1 >>= 1; 165 160 ++exp; 166 } 161 }; 167 162 if (frac1 == 0) { 168 163 /* FIXME : underflow */ 169 170 171 172 } 173 } 164 result.parts.exp = 0; 165 result.parts.fraction = 0; 166 return result; 167 }; 168 }; 174 169 result.parts.exp = exp; 175 result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1);170 result.parts.fraction = frac1 & ( (1 << FLOAT32_FRACTION_SIZE) - 1); 176 171 177 172 return result; 173 178 174 } 179 175 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. 176 /** Multiply two 64 bit float numbers 177 * 186 178 */ 187 179 float64 mulFloat64(float64 a, float64 b) … … 193 185 result.parts.sign = a.parts.sign ^ b.parts.sign; 194 186 195 if (isFloat64NaN(a) || isFloat64NaN(b) ) {187 if (isFloat64NaN(a) || isFloat64NaN(b) ) { 196 188 /* TODO: fix SigNaNs */ 197 189 if (isFloat64SigNaN(a)) { … … 199 191 result.parts.exp = a.parts.exp; 200 192 return result; 201 } 193 }; 202 194 if (isFloat64SigNaN(b)) { /* TODO: fix SigNaN */ 203 195 result.parts.fraction = b.parts.fraction; 204 196 result.parts.exp = b.parts.exp; 205 197 return result; 206 } 198 }; 207 199 /* set NaN as result */ 208 200 result.binary = FLOAT64_NAN; 209 201 return result; 210 } 202 }; 211 203 212 204 if (isFloat64Infinity(a)) { … … 241 233 } else { 242 234 ++exp; 243 } 235 }; 244 236 245 237 frac2 = b.parts.fraction; … … 249 241 } else { 250 242 ++exp; 251 } 243 }; 252 244 253 245 frac1 <<= (64 - FLOAT64_FRACTION_SIZE - 1); 254 246 frac2 <<= (64 - FLOAT64_FRACTION_SIZE - 2); 255 247 256 mul64 (frac1, frac2, &frac1, &frac2);257 258 frac 1 |= (frac2!= 0);259 if (frac 1& (0x1ll << 62)) {260 frac 1<<= 1;248 mul64integers(frac1, frac2, &frac1, &frac2); 249 250 frac2 |= (frac1 != 0); 251 if (frac2 & (0x1ll << 62)) { 252 frac2 <<= 1; 261 253 exp--; 262 254 } 263 255 264 result = finishFloat64(exp, frac 1, result.parts.sign);256 result = finishFloat64(exp, frac2, result.parts.sign); 265 257 return result; 266 258 } 267 259 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) 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) 276 267 { 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; 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; 375 291 } 376 292
Note:
See TracChangeset
for help on using the changeset viewer.
