floating-point - 为什么 π 和 π/10 在表示为 binary64 时似乎具有相同的相对误差？

Question

假设您想快速确定当以 IEEE 754 binary64 格式表示时，π 或 π/10 中哪一个具有最大的相对误差。此外，您手头只有一个 C 编译器。

您可以编写下面的 C 程序，或更紧凑的版本:

#include <stdio.h>
#include <math.h>

volatile long double pil = 3.14159265358979323846L;
volatile double pi = 3.14159265358979323846;

volatile long double tpil = 0.314159265358979323846L;
volatile double tpi = 0.314159265358979323846;

int main() {

  volatile long double abs = pil - pi;

  printf("%La\n%La\n%La\n", pil, (long double)pi, abs);

  printf("pi:   abs err %La -> rel %La\n", abs, abs / pil);

  volatile long double abst = tpil - tpi;
  printf("pi/10: abs err %La -> rel %La\n", abst, abst / tpil);

}

有趣的是，这个程序显示 π 和 π/10 的相对误差是相同的:

0xc.90fdaa22168c235p-2
0xc.90fdaa22168cp-2
0x8.d4p-56
圆周率:绝对错误 0x8.d4p-56 -> rel 0xb.3d85789395e215bp-58
pi/10:绝对错误 0xe.2p-60 -> rel 0xb.3d85789395e215bp-58

我添加了 volatile限定符和中间计算以查看生成的程序集并确保这不是编译器错误。显然不是。

这很奇怪，因为许多其他值 v 不具有 v 的相对误差与 v/10 的相对误差相同的性质。我检查了 0.1 和 0.3。此外，这也不是数量级的问题，因为 3 和 3.5 的相对误差明显不同于它们各自的十分之一。虽然 5, 10, 15, ... 确实与它们各自的十分之一具有相同的相对误差，但这些应该被视为异常(exception)。

现在，程序不会计算相对误差的确切值。它只有 12 位来表示它们(一个符号位加上 64 位有效数和 53 位有效数之间的 11 位差异)。因此，在 4096 年中，π 的绝对误差可能会变成 π/10 的绝对误差的 10 倍。

这似乎仍然是一个不太可能的巧合。当表示为 double 时，在 3 和 3.5 之间取的实常数具有与其十分之一相同的相对误差的先验机会是什么？比我的直觉说的大吗？或者是否有另一种方式来看待它，例如“只要 (double)(π/10) 的有效数有足够的尾随零以使其精确地乘以十”，这似乎更频繁(接近于1/8)？

Answer 1

看待这个问题的更好方法是注意 double 的有效数。最近的 π 是 5 的倍数:

#include <stdio.h>
#include <math.h>
#include <stdint.h>
#include <inttypes.h>
#include <string.h>

volatile long double pil = 3.14159265358979323846L;
volatile double pi = 3.14159265358979323846;

volatile long double tpil = 0.314159265358979323846L;
volatile double tpi = 0.314159265358979323846;

void print_significand(double d) {
  uint64_t significand;

  memcpy(&significand, &d, 8);
  significand &= ((uint64_t)1<<52) - 1;
  significand |= (uint64_t)1<<52;

  printf("%" PRIx64 " %" PRIu64 "\n",
     significand, significand);
}

int main() {

  printf("Significand of (double)pi: ");
  print_significand(pi);

  printf("Significand of (double)(pi/10): ");
  print_significand(tpi);
  …

由此可见:

(double)pi 的有效数:1921fb54442d18 7074237752028440
(double)(pi/10) 的有效数:141b2f769cf0e0 5659390201622752

整数 5659390201622752 正好是 7074237752028440/5 * 4。乘以四用于保持 double表示 π/10 归一化。

因此，给定一个在 2 和 4 之间统一选择的任意实数 r，近似为最近的 double d , d 的有效数有五分之一的机会是 5 的倍数。当发生这种情况时， d/10.0是准确的。这使得 d/10.0 “最接近 r/10 的双倍”这个标题的一个很好的候选者。

候选人有多好？

如果除以 10 将我们从 binade 的顶部(略低于 4)带到 binade 的底部(略高于 0.25)，那么可表示的 double 在 r/10 附近比在 r 附近的密度要小，并且有没有机会找到比 d/10.0 更接近 r/10 的两倍.

如果除以 10 将我们从 binade 的底部(略高于 2)带到 binade 的顶部(略低于 0.25)，那么 double 在 r/10 附近比在 r 附近更密集。 r/10 的双近似值有可能比 d/10.0 更好存在，特别是如果 d已经是 r 的一个不好的近似值。

截止点是 2.5 .实数 π 高于此截止点，因此知道其 double 的有效数近似值可被 5 整除足以推断 π/10 将近似于 π 的近似值的十分之一，并且两个近似值的相对误差将相同。