python - 什么时候使用浮点运算计算上限整数log2会失败？

Question

我很好奇区分这两个函数的第一个输入是什么:

from math import *

def ilog2_ceil_alt(i: int) -> int:
    # DO NOT USE THIS FUNCTION IT'S WRONG
    return ceil(log2(i))

def ilog2_ceil(i: int) -> int:
    # Correct
    if i <= 0:
        raise ValueError("math domain error")
    return (i-1).bit_length()

...

显然，当通过管道中的有限double填充无限大小的整数(其对数)时，由于舍入/截断错误，第一个输入将会失败，但是，我尝试运行此测试代码几分钟，但没有发现问题:

...

def test(i):
    if ilog2_ceil(i) != ilog2_ceil_alt(i):
        return i

def main(start=1):
    import multiprocessing, itertools
    p = multiprocessing.Pool()
    it = p.imap_unordered(test, itertools.count(start), 100)
    return next(i for i in it if i is not None)

if __name__ == '__main__':
    i = main()
    print("Failing case:", i)

我尝试测试各种大值，例如 2**32 和 2**64 + 9999，但在这些上没有失败。

“alt”函数失败的最小(正)整数是多少？

Answer 1

ceil(log2(i)) 的第一个问题是整数 i 首先会被转换为接近 0 的浮点类型(这是错误的方向!)，这里具有 53 位精度。例如，如果 i 为 2⁵³ + 1，它将转换为 2⁵³，并且您将得到 53 而不是 54。

但即使使用较小的值，也可能会出现另一个问题:要考虑的值是那些略大于 2 的幂的值，例如 2ⁿ + k，其中一个小整数 k > 0，因为log2 可能会舍入为整数 n(即使 log2 已正确舍入)，而此时您希望 FP 数略大于 n。因此这将给出ceil(n)，即n而不是n+1。

现在，让我们考虑 k 最坏的情况，即 k = 1。让 p 表示精度(在您的示例中，p = 53 表示 double ，但让我们概括一下)。其中 log2(2ⁿ + 1) = n + log2(1 + 2⁻ⁿ) ≈ n + 0.43·2⁻ⁿ。如果 n 的表示正好需要 q 位，则 ulp 将为 2^q−p。为了得到预期的结果，0.43·2⁻ⁿ 需要大于 1/2 ulp，即 2⁻ⁿ⁻² ⩾ 2^{q−p− 1}，即 n ⩽ p−q−1。

由于 q = ⌈log2(n)⌉，因此条件为 n + ⌈log2(n)⌉ ⩽ p − 1。

但是由于 ⌈log2(n)⌉ 与 n 相比很小，因此 n 的最大值将是 p 的量级，因此通过替换 ⌈log2(n) 可以得到一个近似条件)⌉ 乘以 ⌈log2(p)⌉，即 n ⩽ p − ⌈log2(p)⌉ − 1。

这是一个使用 GNU MPFR 的 C 程序对于每个精度 p，找到第一个失败的 n 值:

#include <stdio.h>
#include <stdlib.h>
#include <mpfr.h>

static void test (long p)
{
  mpfr_t x;

  mpfr_init2 (x, p);
  for (long n = 1; ; n++)
    {
      /* i = 2^n + 1 */
      mpfr_set_ui_2exp (x, 1, n, MPFR_RNDN);
      mpfr_add_ui (x, x, 1, MPFR_RNDN);
      mpfr_log2 (x, x, MPFR_RNDN);
      mpfr_ceil (x, x);
      long r = mpfr_get_si (x, MPFR_RNDN);
      if (r != n + 1)
        {
          printf ("p = %ld, fail for n = %ld\n", p, n);
          break;
        }
    }
  mpfr_clear (x);
}

int main (int argc, char **argv)
{
  long p, pmax;

  if (argc != 2)
    exit (1);
  pmax = strtol (argv[1], NULL, 0);
  for (p = 2; p <= pmax; p++)
    test (p);

  return 0;
}

使用参数 64，可以得到:

p = 2, fail for n = 2
p = 3, fail for n = 3
p = 4, fail for n = 4
p = 5, fail for n = 4
p = 6, fail for n = 5
p = 7, fail for n = 6
p = 8, fail for n = 7
p = 9, fail for n = 8
p = 10, fail for n = 8
p = 11, fail for n = 9
p = 12, fail for n = 10
p = 13, fail for n = 11
p = 14, fail for n = 12
p = 15, fail for n = 13
p = 16, fail for n = 14
p = 17, fail for n = 15
p = 18, fail for n = 16
p = 19, fail for n = 16
p = 20, fail for n = 17
p = 21, fail for n = 18
p = 22, fail for n = 19
p = 23, fail for n = 20
p = 24, fail for n = 21
p = 25, fail for n = 22
p = 26, fail for n = 23
p = 27, fail for n = 24
p = 28, fail for n = 25
p = 29, fail for n = 26
p = 30, fail for n = 27
p = 31, fail for n = 28
p = 32, fail for n = 29
p = 33, fail for n = 30
p = 34, fail for n = 31
p = 35, fail for n = 32
p = 36, fail for n = 32
p = 37, fail for n = 33
p = 38, fail for n = 34
p = 39, fail for n = 35
p = 40, fail for n = 36
p = 41, fail for n = 37
p = 42, fail for n = 38
p = 43, fail for n = 39
p = 44, fail for n = 40
p = 45, fail for n = 41
p = 46, fail for n = 42
p = 47, fail for n = 43
p = 48, fail for n = 44
p = 49, fail for n = 45
p = 50, fail for n = 46
p = 51, fail for n = 47
p = 52, fail for n = 48
p = 53, fail for n = 49
p = 54, fail for n = 50
p = 55, fail for n = 51
p = 56, fail for n = 52
p = 57, fail for n = 53
p = 58, fail for n = 54
p = 59, fail for n = 55
p = 60, fail for n = 56
p = 61, fail for n = 57
p = 62, fail for n = 58
p = 63, fail for n = 59
p = 64, fail for n = 60

编辑:因此，对于 double (p = 53)，如果 log2 正确舍入(或至少具有良好的精度)，则最小失败整数为 2⁴⁹ + 1.