使用成对求和，我需要多少项才能得出明显错误的结果？

Question

使用给定的fp数种类，例如float16，直接构造具有完全错误结果的和。例如，使用python / numpy：

import numpy as np

one = np.float16(1)
ope = np.nextafter(one,one+one)

np.array((ope,one,-one,-one)).cumsum()
# array([1.001, 2.   , 1.   , 0.   ], dtype=float16)

在这里，我们习惯于cumsum强制天真的求和。留给自己的设备numpy使用不同的求和顺序，会得到更好的答案：

np.array((ope,one,-one,-one)).sum()
# 0.000977

以上是基于取消。为了排除此类示例，让我们仅允许使用非否定术语。对于幼稚的求和，给出具有非常错误的求和的示例仍然很容易。以下求和10 ^ 4个相同的项，每个项等于10 ^ -4：

np.full(10**4,10**-4,np.float16).cumsum()
# array([1.0e-04, 2.0e-04, 3.0e-04, ..., 2.5e-01, 2.5e-01, 2.5e-01],
  dtype=float16)

最后一项相差4倍。

同样，允许numpy使用成对求和给出更好的结果：

np.full(10**4,10**-4,np.float16).sum()
# 1.0

可以构造超过成对求和的和。选择低于分辨率为1的eps时，我们可以使用1，eps，0，eps，3x0，eps，7x0，eps，15x0，eps，...，但这涉及到疯狂的术语数量。

我的问题：仅使用float16和非否定项，就需要多少项来从成对求和中得出至少相差2倍的结果。

奖励：同样的问题是“积极”而不是“非消极”。可能吗？

Answer 1

深度1432（因此2 ^ 1432项）足以使真实总和超出计算总和两倍。

我对如何确定所需的术语数少于两个的想法有所了解。

我们使用动态编程来回答以下问题：给定深度d和目标浮点和s，具有成对和的2^d非负float16s 的最大真和是s多少？

让那个数量成为T(d, s)。我们复发

T(0, s) = s,    for all s.
T(d, s) =            max            (T(d-1, a) + T(d-1, b)),    for all d, s.
          a, b : float16(a + b) = s

重复执行的每个步骤都涉及遍历大约2^29组合（因为我们可以假设a ? b，并且负浮点数和特殊值超出了限制），并且所需的深度不会超过10^4Hans和您的答案。对我来说似乎可行。

DP代码：

#include <algorithm>
#include <cstdio>
#include <vector>

using Float16 = int;
using Fixed = unsigned long long;

static constexpr int kExponentBits = 5;
static constexpr int kFractionBits = 10;
static constexpr Float16 kInfinity = ((1 << kExponentBits) - 1)
                                     << kFractionBits;

Fixed FixedFromFloat16(Float16 a) {
  int exponent = a >> kFractionBits;
  if (exponent == 0) {
    return a;
  }
  Float16 fraction = a - (exponent << kFractionBits);
  Float16 significand = (1 << kFractionBits) + fraction;
  return static_cast<Fixed>(significand) << (exponent - 1);
}

bool Plus(Float16 a, Float16 b, Float16* c) {
  Fixed exact_sum = FixedFromFloat16(a) + FixedFromFloat16(b);
  int exponent = 64 - kFractionBits - __builtin_clzll(exact_sum);
  if (exponent <= 0) {
    *c = static_cast<Float16>(exact_sum);
    return true;
  }
  Fixed ulp = Fixed{1} << (exponent - 1);
  Fixed remainder = exact_sum & (ulp - 1);
  Fixed rounded_sum = exact_sum - remainder;
  if (2 * remainder > ulp ||
      (2 * remainder == ulp && (rounded_sum & ulp) != 0)) {
    rounded_sum += ulp;
  }
  exponent = 64 - kFractionBits - __builtin_clzll(rounded_sum);
  if (exponent >= (1 << kExponentBits) - 1) {
    return false;
  }
  Float16 significand = rounded_sum >> (exponent - 1);
  Float16 fraction = significand - (Float16{1} << kFractionBits);
  *c = (exponent << kFractionBits) + fraction;
  return true;
}

int main() {
  std::vector<Fixed> greatest0(kInfinity);
  for (Float16 a = 0; a < kInfinity; a++) {
    greatest0[a] = FixedFromFloat16(a);
  }
  for (int depth = 1; true; depth++) {
    auto greatest1 = greatest0;
    for (Float16 a = 1; a < kInfinity; a++) {
      Fixed greatest0_a = greatest0[a];
      for (Float16 b = a; b < kInfinity; b++) {
        Float16 c;
        if (!Plus(a, b, &c)) {
          continue;
        }
        Fixed& value = greatest1[c];
        value = std::max(value, greatest0_a + greatest0[b]);
      }
    }

    std::vector<double> ratios;
    ratios.reserve(kInfinity - 1);
    for (Float16 a = 1; a < kInfinity; a++) {
      ratios.push_back(greatest1[a] / static_cast<double>(FixedFromFloat16(a)));
    }
    std::printf("depth %d, ratio = %.17g\n", depth,
                *std::max_element(ratios.begin(), ratios.end()));
    greatest0.swap(greatest1);
  }
}

我将运行它并在完成后发布更新。

Answer 2

它将花费大量的术语，以至于实际上是不可能的（如果允许零）或实际上是不可能的（如果由于溢出而不允许零）。维基百科总结了由于尼古拉斯·海厄姆（Nicolas Higham）造成的一些错误范围。由于所有项均为非负项，因此条件数为1，因此n个项的相对误差定为| E _n | / | S _n |。？εlog ₂ n /（1-εlog ₂ n），其中ε是机器ε。要偏离两倍，我们需要| E _n |。？| S _n |，仅当εlog ₂ n？1/2，等于n？float16的2 ^{1 /（} 2ε ^） = 2 ¹⁰²⁴。