Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
412 views
in Technique[技术] by (71.8m points)

floating point - Number of consecutive zeros in the decimal representation of a double

What is the maximum number of consecutive non-leading non-trailing zeros (resp. nines) in the exact decimal representation of an IEEE 754 double-precision number?

Context

Consider the problem of converting a double to decimal, rounding up (resp. down), when the only primitive you are able to use is an existing function that converts to the nearest (correctly rounded to any desired number of digits).

You could get a few additional digits and remove them yourself. For instance, to round 1.875 down to one digit after the dot, you could convert it to the nearest decimal representation with two or three digits after the dot (1.87 or 1.875) and then erase the digits yourself in order to obtain the expected answer, 1.8.

For some numbers and choices of an additional number of digits to print, this method produces the wrong result. For instance, for the double nearest to 0.799999996, converting to decimal, rounding to the nearest, to 2, 3 or 4 digits after the dot produces 0.80, 0.800 and 0.8000. Erasing the additional digits after the conversion produces the result 0.8, when the desired result was 0.7.

There being a finite number of double, there exists a number of additional digits that it is enough to print in the initial conversion in order to always compute the correct result after truncation of the obtained decimal representation. This number is related to the maximal number of nines or zeros that can occur in the exact decimal representation of a double.

Related

This question is related to this question about rounding down in the conversion of a double to decimal, and is dual of this question about the correctly rounded conversion of decimal representations to doubles.

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

[Short version: The answer is 20. Recast the problem in terms of finding good rational approximations to numbers of the form 2^e / 10^d; then use continued fractions to find the best such approximation for each suitable d and e.]

The answer appears to be 20: that is, there are examples of IEEE 754 binary64 floats whose decimal expansion has 20 consecutive zeros, but there are none with 21 consecutive zeros in their decimal expansion (excluding leading and trailing zeros). The same is true for strings of nines.

For the first part, all I need to do is exhibit such a float. The value 0x1.9527560bfbed8p-1000 is exactly representable as a binary64 float, and its decimal expansion contains a string of 20 zeros:

1.4770123739081015758322326613397693800319378788862225686396638475789157389044026850930817635789180868803699741668118826590044503912865915000931065333265410967343958956370955236330760696646247901278074331738806828003156818618589682432778455224012594723731303304343292224317331720902511661748324604219378419442700000000000000000000740694966568985212687104794747958616712153948337746429554804241586090095019654323133732729258896166004754316995632195371041441104566613036026346868128222593894931067078171989365490315525401375255259854894072456336393577718955037826961967325532389800834191597056333066925969522850884268136311674777047673845172073566950098844307658716553833345849153012040436628485227928616281881622762650607683099224232137203216552734375E-301

For the part of the question about nines, the decimal expansion of 0x1.c23142c9da581p-405 contains a string of 20 nines:

2.12818792307269553358078502102171540639252016258831784842556110831434197718043638405555406495645619729155240037555858106390933161420388023706431461384056688295540725831155392678607931808851292893574214797681879999999999999999999941026584542575391157788777223962620780080784703190447744595561259568772261019375946489162743091583251953125E-122

To explain how I found the numbers above, and to show that there are no examples with 21 consecutive zeros, we'll need to work a bit harder. A real number with a long string of 9s or 0s in its decimal expansion has the form (a + eps)*10^d for some integers a and d and real number eps, with a nonzero (we might as well assume a positive) and eps nonzero and small. For example, if 0 < abs(eps) < 10^-10 then a + eps has at least 10 zeros following the decimal point (if eps is positive), or 10 nines following the decimal point (if eps is negative); multiplying by 10^d allows for shifting the location of the string of zeros or nines.

But we're interested in numbers of the above form that are simultaneously representable as an IEEE 754 binary64 float; in other words, numbers that are also of the form b*2^e for integers b and e satisfying 2^52 <= b <= 2^53, with e limited in range (and with some additional restrictions on b once we get into the subnormal range, but we can worry about that later).

So combining this, we're looking for solutions to (a + eps) * 10^d = b * 2^e in integers a, b, d and e such that eps is small, a is positive and 2^52 <= b <= 2^53 (and we'll worry about ranges for d and e later). Rearranging, we get eps / b = 2^e / 10^d - a / b. In other words, we're looking for good rational approximations to 2^e / 10^d, with limited denominator. That's a classic application of continued fractions: given d and e, one can efficiently find the best rational approximation with denominator bounded by 2^53.

So the solution strategy in general is:

for each appropriate d and e:
    find the best rational approximation a / b to 2^e / 10^d with denominator <= 2^53
    if (the error in this rational approximation is small enough):
        # we've got a candidate
        examine the decimal expansion of b*2^e

We have only around 2 thousand values for e to check, and at worst a few hundred d for each such e, so the whole thing is computationally very feasible.

Now to details: what does "small enough" mean? Which d and e are "appropriate"?

As to "small enough": let's say that we're looking for strings of at least 19 zeros or nines, so we're looking for solutions with 0 < abs(eps) <= 10^-19. So it's enough to find, for each d and e, all a and b such that abs(2^e / 10^d - a / b) <= 10^-19 * 2^-52. Note that because of the limit on b there can be only one such fraction a / b; if there were another such a' / b' then we have 1 / 2^106 <= 1 / (b *b') <= abs(a / b - a' / b') <= 2 * 10^-19 * 2^-52, a contradiction. So if such a fraction exists it's necessarily the best rational approximation with the given denominator bound.

For d and e: to cover the binary64 range including subnormals, we want e to range from -1126 to 971 inclusive. If d is too large then 2^e / 10^d will be much smaller than 2^-53 and there's no hope of a solution; d <= 16 + floor(e*log10(2)) is a practical bound. If d is too small (or too negative) then 2^e / 10^d will be an integer and there's no solution; to avoid that, we want d > min(e, 0).

With all that covered, let's write some code. The Python solution is pretty straightforward, thanks in part to the existence of the Fraction.limit_deminator method, which does exactly the job of finding the best rational approximation within limits.

from fractions import Fraction
from itertools import groupby
from math import floor, log10

def longest_run(s, c):
    """Length of the longest run of a given character c in the string s."""
    runs = [list(g) for v, g in groupby(s, lambda k: k == c) if v]
    return max(len(run) for run in runs) if runs else 0

def closest_fraction(d, e):
    """Closest rational to 2**e/10**d with denominator at most 2**53."""
    f = Fraction(2**max(e-d, 0) * 5**max(-d, 0), 2**max(0, d-e) * 5**max(0, d))
    approx = f.limit_denominator(2**53)
    return approx.numerator, approx.denominator

seen = set()
emin = -1126
emax = 971
for e in range(emin, emax+1):
    dmin = min(e, 0) + 1
    dmax = int(floor(e*log10(2))) + 16
    for d in range(dmin, dmax+1):
        num, den = closest_fraction(d, e)
        x = float.fromhex('0x{:x}p{}'.format(den, e))
        # Avoid duplicates.
        if x in seen:
            continue
        seen.add(x)
        digits = '{:.1000e}'.format(x).split('e')[0].replace('.','').strip('0')
        zero_run = longest_run(digits, '0')
        if zero_run >= 20:
            print "{} has {} zeros in its expansion".format(x.hex(), zero_run)
        nine_run = longest_run(digits, '9')
        if nine_run >= 20:
            print "{} has {} nines in its expansion".format(x.hex(), nine_run)

There's plenty of scope for performance improvements there (not using Python's fractions module would be a good start :-); as it stands, it takes a few minutes to run to completion. And here are the results:

0x1.9527560bfbed8p-1000 has 20 zeros in its expansion
0x1.fa712b8efae8ep-997 has 20 zeros in its expansion
0x1.515476ae79b24p-931 has 20 nines in its expansion
0x1.a5a9945a181edp-928 has 20 nines in its expansion
0x1.86049d3311305p-909 has 20 zeros in its expansion
0x1.69c08f3dd8742p-883 has 20 zeros in its expansion
0x1.1b41d80091820p-861 has 20 zeros in its expansion
0x1.62124e00b5e28p-858 has 20 zeros in its expansion
0x1.ba96e180e35b2p-855 has 20 zeros in its expansion
0x1.31c5be6377c48p-786 has 20 zeros in its expansion
0x1.7e372dfc55b5ap-783 has 20 zeros in its expansion
0x1.7e89dc1c3860ap-555 has 20 nines in its expansion
0x1.7e89dc1c3860ap-554 has 20 nines in its expansion
0x1.7e89dc1c3860ap-553 has 20 nines in its expansion
0x1.7e89dc1c3860ap-552 has 20 nines in its expansion
0x1.30bd91ea994cbp-548 has 20 zeros in its expansion
0x1.4a5f9de9ee064p-468 has 20 nines in its expansion
0x1.9cf785646987dp-465 has 20 nines in its expansion
0x1.c23142c9da581p-408 has 20 nines in its expansion
0x1.c23142c9da581p-407 has 20 nines in its expansion
0x1.c23142c9da581p-406 has 20 nines in its expansion
0x1.c23142c9da581p-405 has 20 nines in its expansion
0x1.ba431f4e34be9p+738 has 20 nines in its expansion

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...