С парным суммированием, сколько терминов мне нужно, чтобы получить заметно неправильный результат?

Ответ 1

Глубина 1432 (т.е. 2 ^ 1432 слагаемых) достаточна для того, чтобы истинная сумма превысила вычисленную сумму в два раза.

У меня была идея, как определить количество терминов, которое должно быть меньше, чем в два раза.

Мы используем динамическое программирование, чтобы ответить на следующий вопрос: учитывая глубину d и целевую сумму с плавающей точкой s, какая наибольшая истинная сумма неотрицательного 2^d float16 с попарной суммой 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^4 или около того Ганс и ваш ответ. Кажется возможным для меня.

Код 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);
  }
}

Я запустлю это и опубликую обновление, когда это будет сделано.

Ответ 2

Потребовалось бы так много терминов, что это фактически невозможно (если разрешены нули) или фактически невозможно (если нули не разрешены из-за переполнения). В Википедии обобщены некоторые границы ошибок из-за Николаса Хайама. Поскольку все слагаемые неотрицательны, номер условия равен 1, поэтому относительная погрешность для n слагаемых ограничена как | E _n|/| S _n| ≤ & epsilon; log ₂ n/(1 - & epsilon; log ₂ n), где & epsilon; это машина эпсилон. Чтобы быть вдвое меньше, нам нужно | E _n| ≥ | S _n|, что возможно только в том случае, если & epsilon; log ₂ n ≥ 1/2, что эквивалентно n ≥ 2 ^{1/(2 ε)} = 2 ¹⁰²⁴ для float16.

Ответ 3

Остается вопрос, является ли сумма настолько точной, что вы можете получить относительную ошибку 2 при парном суммировании, если разрешите ноль в сумме (*).

Простой ответ - да, добавив неверную последовательность для cum-суммы с экспоненциальным числом нулей следующим образом (где a1, a2, a3,... an проблематично для нормальной суммы):

a1,
a2,
a3, 0,
a4, 0, 0, 0,
a5, 0, 0, 0, 0, 0, 0, 0,
a6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
...

Он сгенерирует ту же сумму с той же ошибкой округления для попарного суммирования, и вам нужно "только" слагаемые 2**(n-1) вместо n. Таким образом, поскольку члены 10**4 могут генерировать коэффициент 4 для нормального суммирования, то члены 2**(10**4-1) могут давать коэффициент 4 для попарного суммирования.

*: Ответ дэвида Эйстенстата показывает, что при запрете на ноль сумма будет переполнена, прежде чем возникнет такая проблема. (Я предполагаю, что парное суммирование повторяется до конца.)