347{
350
351 if (ieeeExponent == 0)
352 {
353 /* We subtract 2 so that the bounds computation has 2 additional bits. */
355 m2 = ieeeMantissa;
356 }
357 else
358 {
361 }
362
363#if STRICTLY_SHORTEST
364 const bool even = (m2 & 1) == 0;
365 const bool acceptBounds = even;
366#else
367 const bool acceptBounds = false;
368#endif
369
370 /* Step 2: Determine the interval of legal decimal representations. */
372
373 /* Implicit bool -> int conversion. True is 1, false is 0. */
374 const uint32 mmShift = ieeeMantissa != 0 || ieeeExponent <= 1;
375
376 /* We would compute mp and mm like this: */
377 /* uint64 mp = 4 * m2 + 2; */
378 /* uint64 mm = mv - 1 - mmShift; */
379
380 /* Step 3: Convert to a decimal power base using 128-bit arithmetic. */
382 vp,
383 vm;
385 bool vmIsTrailingZeros = false;
386 bool vrIsTrailingZeros = false;
387
388 if (e2 >= 0)
389 {
390 /*
391 * I tried special-casing q == 0, but there was no effect on
392 * performance.
393 *
394 * This expr is slightly faster than max(0, log10Pow2(e2) - 1).
395 */
398 const int32 i = -e2 + q + k;
399
400 e10 = q;
401
403
404 if (q <= 21)
405 {
406 /*
407 * This should use q <= 22, but I think 21 is also safe. Smaller
408 * values may still be safe, but it's more difficult to reason
409 * about them.
410 *
411 * Only one of mp, mv, and mm can be a multiple of 5, if any.
412 */
414
415 if (mvMod5 == 0)
416 {
418 }
419 else if (acceptBounds)
420 {
421 /*----
422 * Same as min(e2 + (~mm & 1), pow5Factor(mm)) >= q
423 * <=> e2 + (~mm & 1) >= q && pow5Factor(mm) >= q
424 * <=> true && pow5Factor(mm) >= q, since e2 >= q.
425 *----
426 */
428 }
429 else
430 {
431 /* Same as min(e2 + 1, pow5Factor(mp)) >= q. */
433 }
434 }
435 }
436 else
437 {
438 /*
439 * This expression is slightly faster than max(0, log10Pow5(-e2) - 1).
440 */
445
446 e10 = q + e2;
447
449
450 if (q <= 1)
451 {
452 /*
453 * {vr,vp,vm} is trailing zeros if {mv,mp,mm} has at least q
454 * trailing 0 bits.
455 */
456 /* mv = 4 * m2, so it always has at least two trailing 0 bits. */
457 vrIsTrailingZeros = true;
458 if (acceptBounds)
459 {
460 /*
461 * mm = mv - 1 - mmShift, so it has 1 trailing 0 bit iff
462 * mmShift == 1.
463 */
464 vmIsTrailingZeros = mmShift == 1;
465 }
466 else
467 {
468 /*
469 * mp = mv + 2, so it always has at least one trailing 0 bit.
470 */
471 --vp;
472 }
473 }
474 else if (q < 63)
475 {
476 /* TODO(ulfjack):Use a tighter bound here. */
477 /*
478 * We need to compute min(ntz(mv), pow5Factor(mv) - e2) >= q - 1
479 */
480 /* <=> ntz(mv) >= q - 1 && pow5Factor(mv) - e2 >= q - 1 */
481 /* <=> ntz(mv) >= q - 1 (e2 is negative and -e2 >= q) */
482 /* <=> (mv & ((1 << (q - 1)) - 1)) == 0 */
483
484 /*
485 * We also need to make sure that the left shift does not
486 * overflow.
487 */
489 }
490 }
491
492 /*
493 * Step 4: Find the shortest decimal representation in the interval of
494 * legal representations.
495 */
497 uint8 lastRemovedDigit = 0;
499
500 /* On average, we remove ~2 digits. */
501 if (vmIsTrailingZeros || vrIsTrailingZeros)
502 {
503 /* General case, which happens rarely (~0.7%). */
504 for (;;)
505 {
508
509 if (vpDiv10 <= vmDiv10)
510 break;
511
515
516 vmIsTrailingZeros &= vmMod10 == 0;
517 vrIsTrailingZeros &= lastRemovedDigit == 0;
518 lastRemovedDigit = (
uint8) vrMod10;
519 vr = vrDiv10;
520 vp = vpDiv10;
521 vm = vmDiv10;
522 ++removed;
523 }
524
525 if (vmIsTrailingZeros)
526 {
527 for (;;)
528 {
531
532 if (vmMod10 != 0)
533 break;
534
538
539 vrIsTrailingZeros &= lastRemovedDigit == 0;
540 lastRemovedDigit = (
uint8) vrMod10;
541 vr = vrDiv10;
542 vp = vpDiv10;
543 vm = vmDiv10;
544 ++removed;
545 }
546 }
547
548 if (vrIsTrailingZeros && lastRemovedDigit == 5 && vr % 2 == 0)
549 {
550 /* Round even if the exact number is .....50..0. */
551 lastRemovedDigit = 4;
552 }
553
554 /*
555 * We need to take vr + 1 if vr is outside bounds or we need to round
556 * up.
557 */
558 output = vr + ((vr == vm && (!acceptBounds || !vmIsTrailingZeros)) || lastRemovedDigit >= 5);
559 }
560 else
561 {
562 /*
563 * Specialized for the common case (~99.3%). Percentages below are
564 * relative to this.
565 */
566 bool roundUp = false;
569
570 if (vpDiv100 > vmDiv100)
571 {
572 /* Optimization:remove two digits at a time(~86.2 %). */
575
576 roundUp = vrMod100 >= 50;
577 vr = vrDiv100;
578 vp = vpDiv100;
579 vm = vmDiv100;
580 removed += 2;
581 }
582
583 /*----
584 * Loop iterations below (approximately), without optimization
585 * above:
586 *
587 * 0: 0.03%, 1: 13.8%, 2: 70.6%, 3: 14.0%, 4: 1.40%, 5: 0.14%,
588 * 6+: 0.02%
589 *
590 * Loop iterations below (approximately), with optimization
591 * above:
592 *
593 * 0: 70.6%, 1: 27.8%, 2: 1.40%, 3: 0.14%, 4+: 0.02%
594 *----
595 */
596 for (;;)
597 {
600
601 if (vpDiv10 <= vmDiv10)
602 break;
603
606
607 roundUp = vrMod10 >= 5;
608 vr = vrDiv10;
609 vp = vpDiv10;
610 vm = vmDiv10;
611 ++removed;
612 }
613
614 /*
615 * We need to take vr + 1 if vr is outside bounds or we need to round
616 * up.
617 */
618 output = vr + (vr == vm || roundUp);
619 }
620
621 const int32 exp = e10 + removed;
622
624
628}
#define DOUBLE_POW5_BITCOUNT
#define DOUBLE_POW5_INV_BITCOUNT
static uint64 mulShiftAll(uint64 m, const uint64 *const mul, const int32 j, uint64 *const vp, uint64 *const vm, const uint32 mmShift)
static bool multipleOfPowerOf5(const uint64 value, const uint32 p)
#define DOUBLE_MANTISSA_BITS
static bool multipleOfPowerOf2(const uint64 value, const uint32 p)
static const uint64 DOUBLE_POW5_SPLIT[326][2]
static const uint64 DOUBLE_POW5_INV_SPLIT[292][2]
static uint64 div10(const uint64 x)
static uint64 div5(const uint64 x)
static uint64 div100(const uint64 x)
static int fd(const char *x, int i)
static uint32 pow5bits(const int32 e)
static int32 log10Pow5(const int32 e)
static int32 log10Pow2(const int32 e)