floating-point - double 十进制表示形式中的连续零个数

Question

IEEE 754双精度数的精确十进制表示形式中，连续无前导无尾零(分别为9)的最大数量是多少？

语境

请考虑将double转换为十进制(向上舍入(向下))的问题，当您唯一可以使用的原语是转换为最接近的现有函数(正确舍入为任何所需的位数)时。

您可能会得到一些额外的数字，然后自己删除。例如，要将1.875舍入到小数点后一位，您可以将其转换为小数点后的小数点表示形式，并在小数点后两位或三位(1.87或1.875)，然后自己擦除数字以获得期望的答案， 1.8。

对于某些数字以及要打印的其他数字的选择，此方法将产生错误的结果。例如，对于最接近double的0.799999996，转换为十进制，四舍五入到最接近的点，在点产生0.80，0.800和0.8000后四舍五入至2、3或4位数字。当所需的结果是0.8时，转换后擦除附加数字将生成结果0.7。
double的数量有限，在初始转换中要打印出足够多的附加数字，以便在截取所获得的十进制表示形式后始终计算出正确的结果。此数字与double的精确十进制表示形式中可能出现的最大9或零有关。

有关

这个问题与this question about rounding down in the conversion of a double to decimal有关，是this question about the correctly rounded conversion of decimal representations to doubles的对偶。

Answer 1

[简短版本:答案为20。根据找到格式为2^e / 10^d的数字的良好有理近似值，重现该问题；然后使用连续分数求出每种合适的d和e的最佳近似值。]

答案似乎是20:也就是说，有一些IEEE 754 binary64浮点数的示例，其十进制扩展名具有20连续的零，但没有十进制扩展名的21连续零(不包括前导和尾随零)。九位数的字符串也是如此。

对于第一部分，我要做的就是展示这样一个 float 。值0x1.9527560bfbed8p-1000可以精确地表示为binary64浮点数，其十进制扩展包含20个零的字符串:

1.4770123739081015758322326613397693800319378788862225686396638475789157389044026850930817635789180868803699741668118826590044503912865915000931065333265410967343958956370955236330760696646247901278074331738806828003156818618589682432778455224012594723731303304343292224317331720902511661748324604219378419442700000000000000000000740694966568985212687104794747958616712153948337746429554804241586090095019654323133732729258896166004754316995632195371041441104566613036026346868128222593894931067078171989365490315525401375255259854894072456336393577718955037826961967325532389800834191597056333066925969522850884268136311674777047673845172073566950098844307658716553833345849153012040436628485227928616281881622762650607683099224232137203216552734375E-301

对于有关9的问题，0x1.c23142c9da581p-405的十进制扩展包含20个9的字符串:

2.1281879230726955335807850210217154063925201625883178484255611083143419771804363840555540649564561972915524003755585810639093316142142802802370643146138405668829829072583115539267860793180885129289357421479768187999999999999999999994516925425-753753577577772772626626800800190190274180125-125

为了说明如何找到上述数字，并说明没有21个连续零的示例，我们需要加倍努力。带有十进制扩展名的9s或0s长字符串的实数的形式为(a + eps)*10^d，表示某些整数a和d以及实数eps，a非零(我们还可以假设a为正数)，而eps非零且很小。例如，如果0 < abs(eps) < 10^-10，则a + eps在小数点后至少有10个零(如果eps为正)，或者在小数点后至少有10个9(如果eps为负)；乘以10^d可以移动零或9字符串的位置。

但是我们对上述形式的数字可以同时表示为IEEE 754 binary64 float有兴趣。换句话说，数字的形式也应为b*2^e和b满足e的整数2^52 <= b <= 2^53，并且e的范围受限制(一旦进入非正规范围，对b会有一些其他限制，但我们稍后会担心)。

因此，结合起来，我们正在寻找整数(a + eps) * 10^d = b * 2^e，a，b和d的e的解决方案，以使eps小，a为正和2^52 <= b <= 2^53(并且我们稍后会担心d和e的范围)。重新排列后，我们得到eps / b = 2^e / 10^d - a / b。换句话说，我们正在寻找分母有限的2^e / 10^d的良好有理近似值。那是连续分数的经典应用:给定d和e，可以有效地找到以2^53为边界的分母的最佳有理逼近。

因此，通常的解决方案策略是:

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

我们只有大约2000个值需要检查e，最糟糕的是每个e需要几百个d，因此整个事情在计算上是非常可行的。

现在来详细说明:“足够小”是什么意思？哪个 d和 e是“合适的”？

至于“足够小”:假设我们正在寻找至少19个零或九的字符串，因此我们正在寻找 0 < abs(eps) <= 10^-19的解决方案。因此，对于每个 d和 e，找到所有 a和 b这样的 abs(2^e / 10^d - a / b) <= 10^-19 * 2^-52就足够了。请注意，由于对 b的限制，只能有一个这样的分数 a / b；如果还有另一个这样的 a' / b'，那么我们就有 1 / 2^106 <= 1 / (b *b') <= abs(a / b - a' / b') <= 2 * 10^-19 * 2^-52，这是一个矛盾。因此，如果存在这样的分数，那么一定是有给定分母界限的最佳有理逼近。

对于 d和 e:要覆盖binary64范围(包括次常态)，我们希望 e的范围从 -1126到 971(包括两端)。如果 d太大，那么 2^e / 10^d将比 2^-53小得多，并且没有解决方案的希望； d <= 16 + floor(e*log10(2))是一个实际的界限。如果 d太小(或太负)，那么 2^e / 10^d将是一个整数，没有解决方案。为了避免这种情况，我们需要 d > min(e, 0)。

覆盖了所有内容之后，让我们编写一些代码。 Python解决方案非常简单明了，部分原因在于 Fraction.limit_deminator方法的存在，该方法确实可以在限制内找到最佳有理逼近。

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)

那里有很多性能改进的余地(不使用Python的 fractions模块将是一个不错的开始:-)；就目前而言，它需要几分钟才能完成。结果如下:

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