Changes in uspace/lib/softfloat/generic/div.c [750636a:c67aff2] in mainline
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
uspace/lib/softfloat/generic/div.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 Division functions. 33 34 */ 34 35 … … 40 41 #include <common.h> 41 42 43 /** 44 * Divide two single-precision floats. 45 * 46 * @param a Nominator. 47 * @param b Denominator. 48 * @return Result of division. 49 */ 42 50 float32 divFloat32(float32 a, float32 b) 43 51 { … … 100 108 return result; 101 109 } 102 103 110 104 111 afrac = a.parts.fraction; … … 110 117 if (aexp == 0) { 111 118 if (afrac == 0) { 112 result.parts.exp = 0; 113 result.parts.fraction = 0; 114 return result; 115 } 119 result.parts.exp = 0; 120 result.parts.fraction = 0; 121 return result; 122 } 123 116 124 /* normalize it*/ 117 118 125 afrac <<= 1; 119 120 while (! (afrac & FLOAT32_HIDDEN_BIT_MASK)) {126 /* afrac is nonzero => it must stop */ 127 while (!(afrac & FLOAT32_HIDDEN_BIT_MASK)) { 121 128 afrac <<= 1; 122 129 aexp--; … … 126 133 if (bexp == 0) { 127 134 bfrac <<= 1; 128 129 while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK)) {135 /* bfrac is nonzero => it must stop */ 136 while (!(bfrac & FLOAT32_HIDDEN_BIT_MASK)) { 130 137 bfrac <<= 1; 131 138 bexp--; … … 133 140 } 134 141 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)) {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)) { 139 146 afrac >>= 1; 140 147 aexp++; … … 144 151 145 152 cfrac = (afrac << 32) / bfrac; 146 if (( cfrac & 0x3F) == 0) {147 cfrac |= ( bfrac * cfrac != afrac << 32);153 if ((cfrac & 0x3F) == 0) { 154 cfrac |= (bfrac * cfrac != afrac << 32); 148 155 } 149 156 … … 151 158 152 159 /* find first nonzero digit and shift result and detect possibly underflow */ 153 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 160 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)))) { 154 161 cexp--; 155 162 cfrac <<= 1; 156 157 } ;163 /* TODO: fix underflow */ 164 } 158 165 159 166 cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ … … 162 169 ++cexp; 163 170 cfrac >>= 1; 164 171 } 165 172 166 173 /* check overflow */ 167 if (cexp >= FLOAT32_MAX_EXPONENT 174 if (cexp >= FLOAT32_MAX_EXPONENT) { 168 175 /* FIXME: overflow, return infinity */ 169 176 result.parts.exp = FLOAT32_MAX_EXPONENT; … … 181 188 cfrac >>= 1; 182 189 while (cexp < 0) { 183 cexp 190 cexp++; 184 191 cfrac >>= 1; 185 } 186 192 } 187 193 } else { 188 result.parts.exp = (uint32_t) cexp;194 result.parts.exp = (uint32_t) cexp; 189 195 } 190 196 … … 194 200 } 195 201 202 /** 203 * Divide two double-precision floats. 204 * 205 * @param a Nominator. 206 * @param b Denominator. 207 * @return Result of division. 208 */ 196 209 float64 divFloat64(float64 a, float64 b) 197 210 { … … 200 213 uint64_t afrac, bfrac, cfrac; 201 214 uint64_t remlo, remhi; 215 uint64_t tmplo, tmphi; 202 216 203 217 result.parts.sign = a.parts.sign ^ b.parts.sign; 204 218 205 219 if (isFloat64NaN(a)) { 206 207 220 if (isFloat64SigNaN(b)) { 208 221 /*FIXME: SigNaN*/ … … 262 275 } 263 276 264 265 277 afrac = a.parts.fraction; 266 278 aexp = a.parts.exp; … … 275 287 return result; 276 288 } 289 277 290 /* normalize it*/ 278 279 291 aexp++; 280 281 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK)) {292 /* afrac is nonzero => it must stop */ 293 while (!(afrac & FLOAT64_HIDDEN_BIT_MASK)) { 282 294 afrac <<= 1; 283 295 aexp--; … … 287 299 if (bexp == 0) { 288 300 bexp++; 289 290 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK)) {301 /* bfrac is nonzero => it must stop */ 302 while (!(bfrac & FLOAT64_HIDDEN_BIT_MASK)) { 291 303 bfrac <<= 1; 292 304 bexp--; … … 294 306 } 295 307 296 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2);297 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK 298 299 if ( bfrac <= (afrac << 1)) {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)) { 300 312 afrac >>= 1; 301 313 aexp++; … … 304 316 cexp = aexp - bexp + FLOAT64_BIAS - 2; 305 317 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; 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); 313 323 314 324 while ((int64_t) remhi < 0) { 315 325 cfrac--; 316 remlo += bfrac; 317 remhi += ( remlo < bfrac ); 318 } 319 cfrac |= ( remlo != 0 ); 326 add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo); 327 } 328 cfrac |= (remlo != 0); 320 329 } 321 330 … … 323 332 result = finishFloat64(cexp, cfrac, result.parts.sign); 324 333 return result; 325 326 334 } 327 335 328 uint64_t divFloat64estim(uint64_t a, uint64_t b) 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) 329 344 { 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 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); 359 535 return result; 360 536 }
Note:
See TracChangeset
for help on using the changeset viewer.