Подтвердить что ты не робот

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

Используя заданный вид чисел 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

Можно построить суммы, которые превосходят парное суммирование. Выбрав eps ниже разрешения 1, мы можем использовать 1, eps, 0, eps, 3x0, eps, 7x0, eps, 15x0, eps,..., но это включает безумное количество терминов.

Мой вопрос: используя float16 и только неотрицательные термины, сколько терминов требуется, чтобы получить из парного суммирования результат, который по меньшей мере в 2 раза больше.

Бонус: тот же вопрос с "положительным" вместо "неотрицательного". Это вообще возможно?

4b9b3361

ОТВЕТЫ

Ответ 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 2 n/(1 - & epsilon; log 2 n), где & epsilon; это машина эпсилон. Чтобы быть вдвое меньше, нам нужно | E n| ≥ | S n|, что возможно только в том случае, если & epsilon; log 2 n ≥ 1/2, что эквивалентно n ≥ 2 1/(2 ε) = 2 1024 для 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 для попарного суммирования.

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