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

Если я хочу взять произведение списка чисел с плавающей запятой, то худшая ситуация/средний случай теряется, добавляя свои журналы, а затем беря exp суммы, а не просто умножая их. Есть ли когда-нибудь случай, когда это более точно?

Ответ 1

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

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

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

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

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

Обратите внимание, что ни один конечный double не имеет логарифма, абсолютное значение которого больше 710 или около того. Это означает, что наши суммы логарифмов, рассчитанные с использованием суммирования Кахана, имеют абсолютную погрешность не более 2130 Н эпсилон. Когда мы оцениваем наши суммы логарифмов, мы получаем что-то не более, чем напр. (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. Я уверен, что есть худшие примеры.