| blob | 6948df7a43317123aed0d523f840f28fa4ec5c50 |
1 // Licensed to the .NET Foundation under one or more agreements.
2 // The .NET Foundation licenses this file to you under the MIT license.
3 // See the LICENSE file in the project root for more information.
7 namespace System
8 {
9 internal static partial class Number
10 {
11 // This is a port of the `Grisu3` implementation here: https://github.com/google/double-conversion/blob/a711666ddd063eb1e4b181a6cb981d39a1fc8bac/double-conversion/fast-dtoa.cc
12 // The backing algorithm and the proofs behind it are described in more detail here: http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf
13 // ========================================================================================================================================
14 //
15 // Overview:
16 //
17 // The general idea behind Grisu3 is to leverage additional bits and cached powers of ten to generate the correct digits.
18 // The algorithm is imprecise for some numbers. Fortunately, the algorithm itself can determine this scenario and gives us
19 // a result indicating success or failure. We must fallback to a different algorithm for the failing scenario.
20 internal static class Grisu3
21 {
27 // 1 / Log2(10)
30 // The minimal and maximal target exponents define the range of w's binary exponent,
31 // where w is the result of multiplying the input by a cached power of ten.
32 //
33 // A different range might be chosen on a different platform, to optimize digit generation,
34 // but a smaller range requires more powers of ten to be cached.
39 {
127 };
130 {
131 CachedPowersMinDecimalExponent,
217 CachedPowersPowerMaxDecimalExponent,
218 };
221 {
309 };
312 {
323 };
326 {
337 {
338 DiyFp w = DiyFp.CreateAndGetBoundaries(v, out DiyFp boundaryMinus, out DiyFp boundaryPlus).Normalize();
339 result = TryRunShortest(in boundaryMinus, in w, in boundaryPlus, number.Digits, out length, out decimalExponent);
340 }
341 else
342 {
345 }
348 {
354 }
357 }
360 {
371 {
372 DiyFp w = DiyFp.CreateAndGetBoundaries(v, out DiyFp boundaryMinus, out DiyFp boundaryPlus).Normalize();
373 result = TryRunShortest(in boundaryMinus, in w, in boundaryPlus, number.Digits, out length, out decimalExponent);
374 }
375 else
376 {
379 }
382 {
388 }
391 }
393 // The counted version of Grisu3 only generates requestedDigits number of digits.
394 // This version does not generate the shortest representation, and with enough requested digits 0.1 will at some point print as 0.9999999...
395 // Grisu3 is too imprecise for real halfway cases (1.5 will not work) and therefore the rounding strategy for halfway cases is irrelevant.
396 private static bool TryRunCounted(in DiyFp w, int requestedDigits, Span<byte> buffer, out int length, out int decimalExponent)
397 {
403 DiyFp tenMk = GetCachedPowerForBinaryExponentRange(tenMkMinimalBinaryExponent, tenMkMaximalBinaryExponent, out int mk);
408 // Note that tenMk is only an approximation of 10^-k.
409 // A DiyFp only contains a 64-bit significand and tenMk is thus only precise up to 64-bits.
411 // The DiyFp.Multiply procedure rounds its result and tenMk is approximated too.
412 // The variable scaledW (as well as scaledBoundaryMinus/Plus) are now off by a small amount.
413 //
414 // In fact, scaledW - (w * 10^k) < 1ulp (unit in last place) of scaledW.
415 // In other words, let f = scaledW.f and e = scaledW.e, then:
416 // (f - 1) * 2^e < (w * 10^k) < (f + 1) * 2^e
420 // We now have (double)(scaledW * 10^-mk).
421 //
422 // DigitGenCounted will generate the first requestedDigits of scaledW and return together with a kappa such that:
423 // scaledW ~= buffer * 10^kappa.
424 //
425 // It will not always be exactly the same since DigitGenCounted only produces a limited number of digits.
427 bool result = TryDigitGenCounted(in scaledW, requestedDigits, buffer, out length, out int kappa);
430 }
432 // Provides a decimal representation of v.
433 // Returns true if it succeeds; otherwise, the result cannot be trusted.
434 //
435 // There will be length digits inside the buffer (not null-terminated).
436 // If the function returns true then:
437 // v == (double)(buffer * 10^decimalExponent)
438 //
439 // The digits in the buffer are the shortest represenation possible (no 0.09999999999999999 instead of 0.1).
440 // The shorter representation will even be chosen if the longer one would be closer to v.
441 //
442 // The last digit will be closest to the actual v.
443 // That is, even if several digits might correctly yield 'v' when read again, the closest will be computed.
444 private static bool TryRunShortest(in DiyFp boundaryMinus, in DiyFp w, in DiyFp boundaryPlus, Span<byte> buffer, out int length, out int decimalExponent)
445 {
446 // boundaryMinus and boundaryPlus are the boundaries between v and its closest floating-point neighbors.
447 // Any number strictly between boundaryMinus and boundaryPlus will round to v when converted to a double.
448 // Grisu3 will never output representations that lie exactly on a boundary.
455 DiyFp tenMk = GetCachedPowerForBinaryExponentRange(tenMkMinimalBinaryExponent, tenMkMaximalBinaryExponent, out int mk);
460 // Note that tenMk is only an approximation of 10^-k.
461 // A DiyFp only contains a 64-bit significan and tenMk is thus only precise up to 64-bits.
463 // The DiyFp.Multiply procedure rounds its result and tenMk is approximated too.
464 // The variable scaledW (as well as scaledBoundaryMinus/Plus) are now off by a small amount.
465 //
466 // In fact, scaledW - (w * 10^k) < 1ulp (unit in last place) of scaledW.
467 // In other words, let f = scaledW.f and e = scaledW.e, then:
468 // (f - 1) * 2^e < (w * 10^k) < (f + 1) * 2^e
473 // In theory, it would be possible to avoid some recomputations by computing the difference between w
474 // and boundaryMinus/Plus (a power of 2) and to compute scaledBoundaryMinus/Plus by subtracting/adding
475 // from scaledW. However, the code becomes much less readable and the speed enhancements are not terrific.
480 // DigitGen will generate the digits of scaledW. Therefore, we have:
481 // v == (double)(scaledW * 10^-mk)
482 //
483 // Set decimalExponent == -mk and pass it to DigitGen and if scaledW is not an integer than it will be updated.
484 // For instance, if scaledW == 1.23 then the buffer will be filled with "123" and the decimalExponent will be decreased by 2.
486 bool result = TryDigitGenShortest(in scaledBoundaryMinus, in scaledW, in scaledBoundaryPlus, buffer, out length, out int kappa);
489 }
491 // Returns the biggest power of ten that is less than or equal to the given number.
492 // We furthermore receive the maximum number of bits 'number' has.
493 //
494 // Returns power == 10^(exponent) such that
495 // power <= number < power * 10
496 // If numberBits == 0, then 0^(0-1) is returned.
497 // The number of bits must be <= 32.
498 //
499 // Preconditions:
500 // number < (1 << (numberBits + 1))
502 {
503 // Inspired by the method for finding an integer log base 10 from here:
504 // http://graphics.stanford.edu/~seander/bithacks.html#IntegerLog10
508 // 1233/4096 is approximately 1/log2(10)
514 // We don't have any guarantees that 2^numberBits <= number
516 {
519 }
523 }
525 // Generates (at most) requestedDigits of input number w.
526 //
527 // w is a floating-point number (DiyFp), consisting of a significand and an exponent.
528 // Its exponent is bounded by MinimalTargetExponent and MaximalTargetExponent, hence:
529 // -60 <= w.e <= -32
530 //
531 // Returns false if it fails, in which case the generated digits in the buffer should not be used.
532 //
533 // Preconditions:
534 // w is correct up to 1 ulp (unit in last place). That is, its error must be strictly less than a unit of its last digit.
535 // MinimalTargetExponent <= w.e <= MaximalTargetExponent
536 //
537 // Postconditions:
538 // Returns false if the procedure fails; otherwise:
539 // * buffer is not null-terminated, but length contains the number of digits.
540 // * The representation in buffer is the most precise representation of requestedDigits digits.
541 // * buffer contains at most requestedDigits digits of w. If there are less than requestedDigits digits then some trailing '0's have been removed.
542 // * kappa is such that w = buffer * 10^kappa + eps with |eps| < 10^kappa / 2.
543 //
544 // This procedure takes into account the imprecision of its input numbers.
545 // If the precision is not enough to guarantee all the postconditions, then false is returned.
546 // This usually happens rarely, but the failure-rate increases with higher requestedDigits
547 private static bool TryDigitGenCounted(in DiyFp w, int requestedDigits, Span<byte> buffer, out int length, out int kappa)
548 {
554 // w is assumed to have an error less than 1 unit.
555 // Whenever w is scaled we also scale its error.
558 // We cut the input number into two parts: the integral digits and the fractional digits.
559 // We don't emit any decimal separator, but adapt kapp instead.
560 // For example: instead of writing "1.2", we put "12" into the buffer and increase kappa by 1.
563 // Division by one is a shift.
566 // Modulo by one is an and.
569 // We deviate from the original algorithm here and do some early checks to determine if we can satisfy requestedDigits.
570 // If we determine that we can't, we exit early and avoid most of the heavy lifting that the algorithm otherwise does.
571 //
572 // When fractionals is zero, we can easily determine if integrals can satisfy requested digits:
573 // If requestedDigits >= 11, integrals is not able to exhaust the count by itself since 10^(11 -1) > uint.MaxValue >= integrals.
574 // If integrals < 10^(requestedDigits - 1), integrals cannot exhaust the count.
575 // Otherwise, integrals might be able to exhaust the count and we need to execute the rest of the code.
576 if ((fractionals == 0) && ((requestedDigits >= 11) || (integrals < s_SmallPowersOfTen[requestedDigits - 1])))
577 {
582 }
587 // Loop invariant:
588 // buffer = w / 10^kappa (integer division)
589 // These invariants hold for the first iteration:
590 // kappa has been initialized with the divisor exponent + 1
591 // The divisor is the biggest power of ten that is smaller than integrals
593 {
598 length++;
599 requestedDigits--;
600 kappa--;
602 // Note that kappa now equals the exponent of the
603 // divisor and that the invariant thus holds again.
605 {
607 }
610 }
613 {
616 buffer,
617 length,
618 rest,
621 ref kappa
622 );
623 }
625 // The integrals have been generated and we are at the point of the decimal separator.
626 // In the following loop, we simply multiply the remaining digits by 10 and divide by one.
627 // We just need to pay attention to multiply associated data (the unit), too.
628 // Note that the multiplication by 10 does not overflow because:
629 // w.e >= -60 and thus one.e >= -60
636 {
640 // Integer division by one.
645 length++;
646 requestedDigits--;
647 kappa--;
649 // Modulo by one.
651 }
654 {
659 }
662 buffer,
663 length,
667 ref kappa
668 );
669 }
671 // Generates the digits of input number w.
672 //
673 // w is a floating-point number (DiyFp), consisting of a significand and an exponent.
674 // Its exponent is bounded by kMinimalTargetExponent and kMaximalTargetExponent, hence:
675 // -60 <= w.e() <= -32.
676 //
677 // Returns false if it fails, in which case the generated digits in the buffer should not be used.
678 //
679 // Preconditions:
680 // low, w and high are correct up to 1 ulp (unit in the last place). That is, their error must be less than a unit of their last digits.
681 // low.e() == w.e() == high.e()
682 // low < w < high, and taking into account their error: low~ <= high~
683 // kMinimalTargetExponent <= w.e() <= kMaximalTargetExponent
684 //
685 // Postconditions:
686 // Returns false if procedure fails; otherwise:
687 // * buffer is not null-terminated, but len contains the number of digits.
688 // * buffer contains the shortest possible decimal digit-sequence such that LOW < buffer * 10^kappa < HIGH, where LOW and HIGH are the correct values of low and high (without their error).
689 // * If more than one decimal representation gives the minimal number of decimal digits then the one closest to W (where W is the correct value of w) is chosen.
690 //
691 // This procedure takes into account the imprecision of its input numbers.
692 // If the precision is not enough to guarantee all the postconditions then false is returned.
693 // This usually happens rarely (~0.5%).
694 //
695 // Say, for the sake of example, that:
696 // w.e() == -48, and w.f() == 0x1234567890abcdef
697 //
698 // w's value can be computed by w.f() * 2^w.e()
699 //
700 // We can obtain w's integral digits by simply shifting w.f() by -w.e().
701 // -> w's integral part is 0x1234
702 // w's fractional part is therefore 0x567890abcdef.
703 //
704 // Printing w's integral part is easy (simply print 0x1234 in decimal).
705 //
706 // In order to print its fraction we repeatedly multiply the fraction by 10 and get each digit.
707 // For example, the first digit after the point would be computed by
708 // (0x567890abcdef * 10) >> 48. -> 3
709 //
710 // The whole thing becomes slightly more complicated because we want to stop once we have enough digits.
711 // That is, once the digits inside the buffer represent 'w' we can stop.
712 //
713 // Everything inside the interval low - high represents w.
714 // However we have to pay attention to low, high and w's imprecision.
715 private static bool TryDigitGenShortest(in DiyFp low, in DiyFp w, in DiyFp high, Span<byte> buffer, out int length, out int kappa)
716 {
725 // low, w, and high are imprecise, but by less than one ulp (unit in the last place).
726 //
727 // If we remove (resp. add) 1 ulp from low (resp. high) we are certain that the new numbers
728 // are outside of the interval we want the final representation to lie in.
729 //
730 // Inversely adding (resp. removing) 1 ulp from low (resp. high) would yield numbers that
731 // are certain to lie in the interval. We will use this fact later on.
732 //
733 // We will now start by generating the digits within the uncertain interval.
734 // Later, we will weed out representations that lie outside the safe interval and thus might lie outside the correct interval.
741 // tooLow and tooHigh are guaranteed to lie outside the interval we want the generated number in.
745 // We now cut the input number into two parts: the integral digits and the fractional digits.
746 // We will not write any decimal separator, but adapt kappa instead.
747 //
748 // Reminder: we are currently computing the digits (Stored inside the buffer) such that:
749 // tooLow < buffer * 10^kappa < tooHigh
750 //
751 // We use tooHigh for the digitGeneration and stop as soon as possible.
752 // If we stop early, we effectively round down.
756 // Division by one is a shift.
759 // Modulo by one is an and.
765 // Loop invariant:
766 // buffer = tooHigh / 10^kappa (integer division)
767 // These invariants hold for the first iteration:
768 // kappa has been initialized with the divisor exponent + 1
769 // The divisor is the biggest power of ten that is smaller than integrals
771 {
776 length++;
777 kappa--;
779 // Note that kappa now equals the exponent of the
780 // divisor and that the invariant thus holds again.
784 // Invariant: tooHigh = buffer * 10^kappa + DiyFp(rest, one.e)
785 // Reminder: unsafeInterval.e == one.e
788 {
789 // Rounding down (by not emitting the remaining digits)
790 // yields a number that lies within the unsafe interval
793 buffer,
794 length,
797 rest,
799 unit
800 );
801 }
804 }
806 // The integrals have been generated and we are at the point of the decimal separator.
807 // In the following loop, we simply multiply the remaining digits by 10 and divide by one.
808 // We just need to pay attention to multiply associated data (the unit), too.
809 // Note that the multiplication by 10 does not overflow because:
810 // w.e >= -60 and thus one.e >= -60
817 {
823 // Integer division by one.
828 length++;
829 kappa--;
831 // Modulo by one.
835 {
837 buffer,
838 length,
843 unit
844 );
845 }
846 }
847 }
849 // Returns a cached power-of-ten with a binary exponent in the range [minExponent; maxExponent] (boundaries included).
850 private static DiyFp GetCachedPowerForBinaryExponentRange(int minExponent, int maxExponent, out int decimalExponent)
851 {
865 }
867 // Rounds the buffer upwards if the result is closer to v by possibly adding 1 to the buffer.
868 // If the precision of the calculation is not sufficient to round correctly, return false.
869 //
870 // The rounding might shift the whole buffer, in which case, the kappy is adjusted.
871 // For example "99", kappa = 3 might become "10", kappa = 4.
872 //
873 // If (2 * rest) > tenKappa then the buffer needs to be round up.
874 // rest can have an error of +/- 1 unit.
875 // This function accounts for the imprecision and returns false if the rounding direction cannot be unambiguously determined.
876 //
877 // Preconditions:
878 // rest < tenKappa
879 private static bool TryRoundWeedCounted(Span<byte> buffer, int length, ulong rest, ulong tenKappa, ulong unit, ref int kappa)
880 {
883 // The following tests are done in a specific order to avoid overflows.
884 // They will work correctly with any ulong values of rest < tenKappa and unit.
885 //
886 // If the unit is too big, then we don't know which way to round.
887 // For example, a unit of 50 means that the real number lies within rest +/- 50.
888 // If 10^kappa == 40, then there is no way to tell which way to round.
889 //
890 // Even if unit is just half the size of 10^kappa we are already completely lost.
891 // And after the previous test, we know that the expression will not over/underflow.
893 {
895 }
897 // If 2 * (rest + unit) <= 10^kappa, we can safely round down.
899 {
901 }
903 // If 2 * (rest - unit) >= 10^kappa, we can safely round up.
904 if ((rest > unit) && (tenKappa <= (rest - unit) || ((tenKappa - (rest - unit)) <= (rest - unit))))
905 {
906 // Increment the last digit recursively until we find a non '9' digit.
910 {
912 {
914 }
918 }
920 // If the first digit is now '0'+10, we had a buffer with all '9's.
921 // With the exception of the first digit, all digits are now '0'.
922 // Simply switch the first digit to '1' and adjust the kappa.
923 // For example, "99" becomes "10" and the power (the kappa) is increased.
925 {
927 kappa++;
928 }
931 }
934 }
936 // Adjusts the last digit of the generated number and screens out generated solutions that may be inaccurate.
937 // A solution may be inaccurate if it is outside the safe interval or if we cannot provide that it is closer to the input than a neighboring representation of the same length.
938 //
939 // Input:
940 // buffer containing the digits of tooHigh / 10^kappa
941 // the buffer's length
942 // distanceTooHighW == (tooHigh - w).f * unit
943 // unsafeInterval == (tooHigh - tooLow).f * unit
944 // rest = (tooHigh - buffer * 10^kapp).f * unit
945 // tenKappa = 10^kappa * unit
946 // unit = the common multiplier
947 //
948 // Output:
949 // Returns true if the buffer is guaranteed to contain the closest representable number to the input.
950 //
951 // Modifies the generated digits in the buffer to approach (round towards) w.
952 private static bool TryRoundWeedShortest(Span<byte> buffer, int length, ulong distanceTooHighW, ulong unsafeInterval, ulong rest, ulong tenKappa, ulong unit)
953 {
957 // Let wLow = tooHigh - bigDistance, and wHigh = tooHigh - smallDistance.
958 //
959 // Note: wLow < w < wHigh
960 //
961 // The real w * unit must lie somewhere inside the interval
962 // ]w_low; w_high[ (often written as "(w_low; w_high)")
964 // Basically the buffer currently contains a number in the unsafe interval
965 // ]too_low; too_high[ with too_low < w < too_high
966 //
967 // tooHigh - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
968 // ^v 1 unit ^ ^ ^ ^
969 // boundaryHigh --------------------- . . . .
970 // ^v 1 unit . . . .
971 // - - - - - - - - - - - - - - - - - - - + - - + - - - - - - . .
972 // . . ^ . .
973 // . bigDistance . . .
974 // . . . . rest
975 // smallDistance . . . .
976 // v . . . .
977 // wHigh - - - - - - - - - - - - - - - - - - . . . .
978 // ^v 1 unit . . . .
979 // w --------------------------------------- . . . .
980 // ^v 1 unit v . . .
981 // wLow - - - - - - - - - - - - - - - - - - - - - . . .
982 // . . v
983 // buffer -------------------------------------------------+-------+--------
984 // . .
985 // safeInterval .
986 // v .
987 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - .
988 // ^v 1 unit .
989 // boundaryLow ------------------------- unsafeInterval
990 // ^v 1 unit v
991 // tooLow - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
992 //
993 //
994 // Note that the value of buffer could lie anywhere inside the range tooLow to tooHigh.
995 //
996 // boundaryLow, boundaryHigh and w are approximations of the real boundaries and v (the input number).
997 // They are guaranteed to be precise up to one unit.
998 // In fact the error is guaranteed to be strictly less than one unit.
999 //
1000 // Anything that lies outside the unsafe interval is guaranteed not to round to v when read again.
1001 // Anything that lies inside the safe interval is guaranteed to round to v when read again.
1002 //
1003 // If the number inside the buffer lies inside the unsafe interval but not inside the safe interval
1004 // then we simply do not know and bail out (returning false).
1005 //
1006 // Similarly we have to take into account the imprecision of 'w' when finding the closest representation of 'w'.
1007 // If we have two potential representations, and one is closer to both wLow and wHigh, then we know it is closer to the actual value v.
1008 //
1009 // By generating the digits of tooHigh we got the largest (closest to tooHigh) buffer that is still in the unsafe interval.
1010 // In the case where wHigh < buffer < tooHigh we try to decrement the buffer.
1011 // This way the buffer approaches (rounds towards) w.
1012 //
1013 // There are 3 conditions that stop the decrementation process:
1014 // 1) the buffer is already below wHigh
1015 // 2) decrementing the buffer would make it leave the unsafe interval
1016 // 3) decrementing the buffer would yield a number below wHigh and farther away than the current number.
1017 //
1018 // In other words:
1019 // (buffer{-1} < wHigh) && wHigh - buffer{-1} > buffer - wHigh
1020 //
1021 // Instead of using the buffer directly we use its distance to tooHigh.
1022 //
1023 // Conceptually rest ~= tooHigh - buffer
1024 //
1025 // We need to do the following tests in this order to avoid over- and underflows.
1029 while ((rest < smallDistance) && ((unsafeInterval - rest) >= tenKappa) && (((rest + tenKappa) < smallDistance) || ((smallDistance - rest) >= (rest + tenKappa - smallDistance))))
1030 {
1033 }
1035 // We have approached w+ as much as possible.
1036 // We now test if approaching w- would require changing the buffer.
1037 // If yes, then we have two possible representations close to w, but we cannot decide which one is closer.
1038 if ((rest < bigDistance) && ((unsafeInterval - rest) >= tenKappa) && (((rest + tenKappa) < bigDistance) || ((bigDistance - rest) > (rest + tenKappa - bigDistance))))
1039 {
1041 }
1043 // Weeding test.
1044 //
1045 // The safe interval is [tooLow + 2 ulp; tooHigh - 2 ulp]
1046 // Since tooLow = tooHigh - unsafeInterval this is equivalent to
1047 // [tooHigh - unsafeInterval + 4 ulp; tooHigh - 2 ulp]
1048 //
1049 // Conceptually we have: rest ~= tooHigh - buffer
1051 }
1052 }
1053 }
1054 }