Получение журналов и добавление против умножения

При отсутствии каких-либо махинаций с переполнением или недостатком, если a и b являются числами с плавающей запятой, то произведение a*b будет вычислено с относительной погрешностью в 1/2 ulp.

Таким образом, грубая оценка относительной ошибки после умножения цепочки из N doubles приводит к отклонению ответа не более чем в (1 - эпсилон/2)^-N, что примерно равно exp(эпсилон N/ 2). Я полагаю, что в среднем вы можете ожидать отклонение около эпсилон sqrt (N). (В первом порядке это примерно N эпсилон.)

Однако с этой стратегией более вероятно возникновение переполнения и потери экспоненты; вы, скорее всего, получите бесконечности, нули и NaN, а также неточные значения из-за округления субнормальных величин.

Другой подход более надежен в этом смысле, но он намного медленнее и допускает больше ошибок в случае, когда прямой подход не приводит к переполнению или потере значимости. Вот очень-очень грубый анализ для стандартных удвоений в случае, когда N по крайней мере на пару порядков меньше, чем 2⁵³:

Вы всегда можете взять логарифм конечного числа с плавающей запятой и получить конечное число с плавающей запятой, так что у нас все в порядке. Вы можете сложить N числа с плавающей запятой либо напрямую, чтобы получить N эпсилон наихудшего случая "относительную" ошибку и sqrt(N) эпсилон ожидаемую "относительную" ошибку, либо используя суммирование Кэхана, чтобы получить "относительную" ошибку примерно в 3 эпсилон в худшем случае. Страшные кавычки «относительны», потому что ошибка относится к сумме абсолютных значений вещей, которые вы суммируете.

Обратите внимание, что никакое конечное double не имеет логарифма, абсолютное значение которого больше 710 или около того. Это означает, что наша сумма логарифмов, вычисленная с помощью суммирования Кэхана, имеет абсолютную ошибку не более 2130 Н эпсилон. Когда мы возводим в степень нашу сумму логарифмов, мы получаем что-то не более чем в exp(2130 N эпсилон) от правильного ответа.

Патологические примеры для подхода log-sum-exp:

int main() {
  double foo[] = {0x1.000000000018cp1023, 0x1.0000000000072p-1023};
  double prod = 1;
  double sumlogs = 0;
  for (int i = 0; i < sizeof(foo) / sizeof(*foo); i++) {
    prod *= foo[i];
    sumlogs += log(foo[i]);
  }
  printf("%a %a\n", foo[0], foo[1]);
  printf("%a %a %a\n", prod, exp(sumlogs), prod - exp(sumlogs));
}

На моей платформе я получаю разницу 0x1.fep-44. Я уверен, что есть примеры и похуже.