Changes in uspace/lib/softfloat/generic/div.c [c67aff2:750636a] in mainline
- File:
-
- 1 edited
-
uspace/lib/softfloat/generic/div.c (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
uspace/lib/softfloat/generic/div.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 Division functions.32 /** @file 34 33 */ 35 34 … … 41 40 #include <common.h> 42 41 43 /**44 * Divide two single-precision floats.45 *46 * @param a Nominator.47 * @param b Denominator.48 * @return Result of division.49 */50 42 float32 divFloat32(float32 a, float32 b) 51 43 { … … 108 100 return result; 109 101 } 102 110 103 111 104 afrac = a.parts.fraction; … … 117 110 if (aexp == 0) { 118 111 if (afrac == 0) { 119 result.parts.exp = 0; 120 result.parts.fraction = 0; 121 return result; 122 } 123 112 result.parts.exp = 0; 113 result.parts.fraction = 0; 114 return result; 115 } 124 116 /* normalize it*/ 117 125 118 afrac <<= 1; 126 /* afrac is nonzero => it must stop */127 while (! (afrac & FLOAT32_HIDDEN_BIT_MASK)) {119 /* afrac is nonzero => it must stop */ 120 while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) { 128 121 afrac <<= 1; 129 122 aexp--; … … 133 126 if (bexp == 0) { 134 127 bfrac <<= 1; 135 /* bfrac is nonzero => it must stop */136 while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK)) {128 /* bfrac is nonzero => it must stop */ 129 while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) { 137 130 bfrac <<= 1; 138 131 bexp--; … … 140 133 } 141 134 142 afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1);143 bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE);144 145 if ( bfrac <= (afrac << 1)) {135 afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 ); 136 bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE ); 137 138 if ( bfrac <= (afrac << 1) ) { 146 139 afrac >>= 1; 147 140 aexp++; … … 151 144 152 145 cfrac = (afrac << 32) / bfrac; 153 if (( cfrac & 0x3F) == 0) {154 cfrac |= ( bfrac * cfrac != afrac << 32);146 if (( cfrac & 0x3F ) == 0) { 147 cfrac |= ( bfrac * cfrac != afrac << 32 ); 155 148 } 156 149 … … 158 151 159 152 /* find first nonzero digit and shift result and detect possibly underflow */ 160 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {153 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { 161 154 cexp--; 162 155 cfrac <<= 1; 163 /* TODO: fix underflow */164 } 156 /* TODO: fix underflow */ 157 }; 165 158 166 159 cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ … … 169 162 ++cexp; 170 163 cfrac >>= 1; 171 }164 } 172 165 173 166 /* check overflow */ 174 if (cexp >= FLOAT32_MAX_EXPONENT ) {167 if (cexp >= FLOAT32_MAX_EXPONENT ) { 175 168 /* FIXME: overflow, return infinity */ 176 169 result.parts.exp = FLOAT32_MAX_EXPONENT; … … 188 181 cfrac >>= 1; 189 182 while (cexp < 0) { 190 cexp ++;183 cexp ++; 191 184 cfrac >>= 1; 192 } 185 } 186 193 187 } else { 194 result.parts.exp = (uint32_t) 188 result.parts.exp = (uint32_t)cexp; 195 189 } 196 190 … … 200 194 } 201 195 202 /**203 * Divide two double-precision floats.204 *205 * @param a Nominator.206 * @param b Denominator.207 * @return Result of division.208 */209 196 float64 divFloat64(float64 a, float64 b) 210 197 { … … 213 200 uint64_t afrac, bfrac, cfrac; 214 201 uint64_t remlo, remhi; 215 uint64_t tmplo, tmphi;216 202 217 203 result.parts.sign = a.parts.sign ^ b.parts.sign; 218 204 219 205 if (isFloat64NaN(a)) { 206 220 207 if (isFloat64SigNaN(b)) { 221 208 /*FIXME: SigNaN*/ … … 275 262 } 276 263 264 277 265 afrac = a.parts.fraction; 278 266 aexp = a.parts.exp; … … 287 275 return result; 288 276 } 289 290 277 /* normalize it*/ 278 291 279 aexp++; 292 /* afrac is nonzero => it must stop */293 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK)) {280 /* afrac is nonzero => it must stop */ 281 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) { 294 282 afrac <<= 1; 295 283 aexp--; … … 299 287 if (bexp == 0) { 300 288 bexp++; 301 /* bfrac is nonzero => it must stop */302 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK)) {289 /* bfrac is nonzero => it must stop */ 290 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) { 303 291 bfrac <<= 1; 304 292 bexp--; … … 306 294 } 307 295 308 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2);309 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);310 311 if ( bfrac <= (afrac << 1)) {296 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 ); 297 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1); 298 299 if ( bfrac <= (afrac << 1) ) { 312 300 afrac >>= 1; 313 301 aexp++; … … 316 304 cexp = aexp - bexp + FLOAT64_BIAS - 2; 317 305 318 cfrac = div128est(afrac, 0x0ll, bfrac); 319 320 if ((cfrac & 0x1FF) <= 2) { 321 mul64(bfrac, cfrac, &tmphi, &tmplo); 322 sub128(afrac, 0x0ll, tmphi, tmplo, &remhi, &remlo); 306 cfrac = divFloat64estim(afrac, bfrac); 307 308 if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */ 309 mul64integers( bfrac, cfrac, &remlo, &remhi); 310 /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ 311 remhi = afrac - remhi - ( remlo > 0); 312 remlo = - remlo; 323 313 324 314 while ((int64_t) remhi < 0) { 325 315 cfrac--; 326 add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo); 327 } 328 cfrac |= (remlo != 0); 316 remlo += bfrac; 317 remhi += ( remlo < bfrac ); 318 } 319 cfrac |= ( remlo != 0 ); 329 320 } 330 321 … … 332 323 result = finishFloat64(cexp, cfrac, result.parts.sign); 333 324 return result; 325 334 326 } 335 327 336 /** 337 * Divide two quadruple-precision floats. 338 * 339 * @param a Nominator. 340 * @param b Denominator. 341 * @return Result of division. 342 */ 343 float128 divFloat128(float128 a, float128 b) 328 uint64_t divFloat64estim(uint64_t a, uint64_t b) 344 329 { 345 float128 result; 346 int64_t aexp, bexp, cexp; 347 uint64_t afrac_hi, afrac_lo, bfrac_hi, bfrac_lo, cfrac_hi, cfrac_lo; 348 uint64_t shift_out; 349 uint64_t rem_hihi, rem_hilo, rem_lohi, rem_lolo; 350 uint64_t tmp_hihi, tmp_hilo, tmp_lohi, tmp_lolo; 351 352 result.parts.sign = a.parts.sign ^ b.parts.sign; 353 354 if (isFloat128NaN(a)) { 355 if (isFloat128SigNaN(b)) { 356 /*FIXME: SigNaN*/ 357 return b; 358 } 359 360 if (isFloat128SigNaN(a)) { 361 /*FIXME: SigNaN*/ 362 } 363 /*NaN*/ 364 return a; 365 } 366 367 if (isFloat128NaN(b)) { 368 if (isFloat128SigNaN(b)) { 369 /*FIXME: SigNaN*/ 370 } 371 /*NaN*/ 372 return b; 373 } 374 375 if (isFloat128Infinity(a)) { 376 if (isFloat128Infinity(b) || isFloat128Zero(b)) { 377 /*FIXME: inf / inf */ 378 result.binary.hi = FLOAT128_NAN_HI; 379 result.binary.lo = FLOAT128_NAN_LO; 380 return result; 381 } 382 /* inf / num */ 383 result.parts.exp = a.parts.exp; 384 result.parts.frac_hi = a.parts.frac_hi; 385 result.parts.frac_lo = a.parts.frac_lo; 386 return result; 387 } 388 389 if (isFloat128Infinity(b)) { 390 if (isFloat128Zero(a)) { 391 /* FIXME 0 / inf */ 392 result.parts.exp = 0; 393 result.parts.frac_hi = 0; 394 result.parts.frac_lo = 0; 395 return result; 396 } 397 /* FIXME: num / inf*/ 398 result.parts.exp = 0; 399 result.parts.frac_hi = 0; 400 result.parts.frac_lo = 0; 401 return result; 402 } 403 404 if (isFloat128Zero(b)) { 405 if (isFloat128Zero(a)) { 406 /*FIXME: 0 / 0*/ 407 result.binary.hi = FLOAT128_NAN_HI; 408 result.binary.lo = FLOAT128_NAN_LO; 409 return result; 410 } 411 /* FIXME: division by zero */ 412 result.parts.exp = 0; 413 result.parts.frac_hi = 0; 414 result.parts.frac_lo = 0; 415 return result; 416 } 417 418 afrac_hi = a.parts.frac_hi; 419 afrac_lo = a.parts.frac_lo; 420 aexp = a.parts.exp; 421 bfrac_hi = b.parts.frac_hi; 422 bfrac_lo = b.parts.frac_lo; 423 bexp = b.parts.exp; 424 425 /* denormalized numbers */ 426 if (aexp == 0) { 427 if (eq128(afrac_hi, afrac_lo, 0x0ll, 0x0ll)) { 428 result.parts.exp = 0; 429 result.parts.frac_hi = 0; 430 result.parts.frac_lo = 0; 431 return result; 432 } 433 434 /* normalize it*/ 435 aexp++; 436 /* afrac is nonzero => it must stop */ 437 and128(afrac_hi, afrac_lo, 438 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 439 &tmp_hihi, &tmp_lolo); 440 while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) { 441 lshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo); 442 aexp--; 443 } 444 } 445 446 if (bexp == 0) { 447 bexp++; 448 /* bfrac is nonzero => it must stop */ 449 and128(bfrac_hi, bfrac_lo, 450 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 451 &tmp_hihi, &tmp_lolo); 452 while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) { 453 lshift128(bfrac_hi, bfrac_lo, 1, &bfrac_hi, &bfrac_lo); 454 bexp--; 455 } 456 } 457 458 or128(afrac_hi, afrac_lo, 459 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 460 &afrac_hi, &afrac_lo); 461 lshift128(afrac_hi, afrac_lo, 462 (128 - FLOAT128_FRACTION_SIZE - 1), &afrac_hi, &afrac_lo); 463 or128(bfrac_hi, bfrac_lo, 464 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 465 &bfrac_hi, &bfrac_lo); 466 lshift128(bfrac_hi, bfrac_lo, 467 (128 - FLOAT128_FRACTION_SIZE - 1), &bfrac_hi, &bfrac_lo); 468 469 if (le128(bfrac_hi, bfrac_lo, afrac_hi, afrac_lo)) { 470 rshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo); 471 aexp++; 472 } 473 474 cexp = aexp - bexp + FLOAT128_BIAS - 2; 475 476 cfrac_hi = div128est(afrac_hi, afrac_lo, bfrac_hi); 477 478 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_hi, 479 &tmp_lolo /* dummy */, &tmp_hihi, &tmp_hilo, &tmp_lohi); 480 481 /* sub192(afrac_hi, afrac_lo, 0, 482 * tmp_hihi, tmp_hilo, tmp_lohi 483 * &rem_hihi, &rem_hilo, &rem_lohi); */ 484 sub128(afrac_hi, afrac_lo, tmp_hihi, tmp_hilo, &rem_hihi, &rem_hilo); 485 if (tmp_lohi > 0) { 486 sub128(rem_hihi, rem_hilo, 0x0ll, 0x1ll, &rem_hihi, &rem_hilo); 487 } 488 rem_lohi = -tmp_lohi; 489 490 while ((int64_t) rem_hihi < 0) { 491 --cfrac_hi; 492 /* add192(rem_hihi, rem_hilo, rem_lohi, 493 * 0, bfrac_hi, bfrac_lo, 494 * &rem_hihi, &rem_hilo, &rem_lohi); */ 495 add128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo, &rem_hilo, &rem_lohi); 496 if (lt128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo)) { 497 ++rem_hihi; 498 } 499 } 500 501 cfrac_lo = div128est(rem_hilo, rem_lohi, bfrac_lo); 502 503 if ((cfrac_lo & 0x3FFF) <= 4) { 504 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_lo, 505 &tmp_hihi /* dummy */, &tmp_hilo, &tmp_lohi, &tmp_lolo); 506 507 /* sub192(rem_hilo, rem_lohi, 0, 508 * tmp_hilo, tmp_lohi, tmp_lolo, 509 * &rem_hilo, &rem_lohi, &rem_lolo); */ 510 sub128(rem_hilo, rem_lohi, tmp_hilo, tmp_lohi, &rem_hilo, &rem_lohi); 511 if (tmp_lolo > 0) { 512 sub128(rem_hilo, rem_lohi, 0x0ll, 0x1ll, &rem_hilo, &rem_lohi); 513 } 514 rem_lolo = -tmp_lolo; 515 516 while ((int64_t) rem_hilo < 0) { 517 --cfrac_lo; 518 /* add192(rem_hilo, rem_lohi, rem_lolo, 519 * 0, bfrac_hi, bfrac_lo, 520 * &rem_hilo, &rem_lohi, &rem_lolo); */ 521 add128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo, &rem_lohi, &rem_lolo); 522 if (lt128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo)) { 523 ++rem_hilo; 524 } 525 } 526 527 cfrac_lo |= ((rem_hilo | rem_lohi | rem_lolo) != 0 ); 528 } 529 530 shift_out = cfrac_lo << (64 - (128 - FLOAT128_FRACTION_SIZE - 1)); 531 rshift128(cfrac_hi, cfrac_lo, (128 - FLOAT128_FRACTION_SIZE - 1), 532 &cfrac_hi, &cfrac_lo); 533 534 result = finishFloat128(cexp, cfrac_hi, cfrac_lo, result.parts.sign, shift_out); 330 uint64_t bhi; 331 uint64_t remhi, remlo; 332 uint64_t result; 333 334 if ( b <= a ) { 335 return 0xFFFFFFFFFFFFFFFFull; 336 } 337 338 bhi = b >> 32; 339 result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32; 340 mul64integers(b, result, &remlo, &remhi); 341 342 remhi = a - remhi - (remlo > 0); 343 remlo = - remlo; 344 345 b <<= 32; 346 while ( (int64_t) remhi < 0 ) { 347 result -= 0x1ll << 32; 348 remlo += b; 349 remhi += bhi + ( remlo < b ); 350 } 351 remhi = (remhi << 32) | (remlo >> 32); 352 if (( bhi << 32) <= remhi) { 353 result |= 0xFFFFFFFF; 354 } else { 355 result |= remhi / bhi; 356 } 357 358 535 359 return result; 536 360 }
Note:
See TracChangeset
for help on using the changeset viewer.
