Наиболее точное скалярное произведение векторов типа double. Вычисление значения полинома28.11.2020 13:17

В прошлой статье речь шла о том, как можно сложить массив из чисел типа double наиболее точно, то есть получить такую сумму, как если бы мы выполняли вычисления в рамках арифметики с бесконечной точностью, а затем один раз округлили бы результат. Был показан алгоритм, который эквивалентен применению типа данных double-double, в котором сложение происходит сразу в двух переменных: основная сумма и хвостик-погрешность. Опытные читатели сразу догадались, что сложение хвостиков-погрешностей также допускает по отношению к себе рекурсивное применение того же алгоритма, что приводит не к удвоенной, а к утроенной точности, и вообще, можно организовать каскад сложений произвольного размера, получая любую наперёд заданную точность расчётов, поэтому фактически в прошлой статье была показана предпосылка к так называемой «дробной длинной арифметике». Опытный программист без труда разберётся как её реализовать, ну, а я обещал дать аналогичные фундаментальные основы для скалярного произведения и вычисления полинома в точке. Поскольку все базовые вводные слова уже были сказаны в двух предшествующих статьях, в этой будет меньше «воды» и «лишних», по мнению опытных математиков, сведений. Прошу под кат.

Как обычно, видео-презентация полностью дублирует содержание статьи, и создана для удобства тех посетителей Хабра, кому голос с иллюстрациями воспринять проще чем текст без таковых.

О сложении мы знаем уже достаточно: с помощью функции two_sum можно сложить N чисел наиболее точно, дело осталось за малым: научиться перемножать два числа типа double так, чтобы не терять отбрасываемый при округлении «хвостик», сохраняя его в отдельной переменной. Если в случае сложения сделать это можно очень просто при любом раскладе, лишь бы промежуточная сумма не переполнялась (кто не читал: «Сложение двух чисел с плавающей запятой без потери точности»), то в случае умножения есть дополнительное ограничение.

Обозначим через  минимальную экспоненту вашей системы чисел с плавающей запятой, в рамках которой выполняется умножение (для float: −126, для double: −1022), буквой обозначим количество битов мантиссы с учётом скрытой единицы (24 для float и 53 для double). Пусть и  — экспоненты наших исходных чисел и , тогда требуется, чтобы было выполнено условие:

Если данное условие не выполняется, то могут возникнуть случаи, когда нельзя представить в виде суммы , где и . Указанное неравенство является сложным способом сказать простую фразу: «хвостик» не должен быть денормализованным. И это вполне понятно: когда мы перемножаем два числа, получаем удвоенное количество битов, которые гарантированно влезают в две переменные того же типа данных. Если же одна из этих переменных попадает в денормализованный диапазон, происходит битовый сдвиг мантиссы с потерей младших битов. И кроме переполнения это единственная ситуация, при которой произведение двух чисел нельзя точно представить в форме суммы двух чисел того же типа.

Так ли существенно это ограничение? Для float нужно чтобы сумма экспонент оказалась не меньше −126+24−1=−103, то есть, грубо говоря, если каждый из множителей больше 2−51, то всё будет хорошо. Спросите себя: часто ли вы применяете числа такого порядка? Если да, то можете выполнить масштабирование и перейти к числам побольше. Можете применить double, и работать с числами больше 2−485. Но даже если плюнуть на всё и применить алгоритм в неподходящих условиях, в простых прикладных задачах он всё равно отработает правильно чуть чаще чем всегда.

Если указанное условие выполнено, то с помощью инструкции FMA (Fused Multiply-Add) точное умножение двух чисел выполняется следующим алгоритмом:

double __fastcall two_prod (double &t, double a, double b) {
  double p = a*b;
  t = fma (a, b, -p);  // t = a*b-p, не забудьте подключить 
  return p;
}

Результатом работы функции будет пара чисел , такая, что  является ближайшим к  числом, а  — хвостиком-погрешностью, причём в точности. Разумеется, по ходу расчётов не должно произойти переполнения, за этим придётся следить отдельно: нужно оставаться в рамках условия .

Если операция FMA для вас недоступна, то ситуация сильно усложняется. Желаемое произведение можно получить путём разбиения исходных переменных и  на две части примерно посередине так, чтобы старшая часть операндов (ah и bh) содержала половину старших битов, а младшая (al и bl) — половину младших. Далее нужные операции выполняются именно с этими половинками, в общей сложности алгоритм потребует 17 операций сложения и умножения. Данный трюк в чём-то аналогичен известному алгоритмическому трюку для вычисления произведения двух целых чисел с сохранением результата в удвоенном типе данных (кто не знал, читайте Г. Уоррен «Алгоритмические трюки для программистов», в издании 2014 года это раздел 8.2). Нижеследующий алгоритм написан мною по мотивам книги Jean-Michel Muller, «Handbook of floating-point arithmetic», также как и остальная часть этой статьи.

    // Разбиение x на старшую и младшую часть
void __fastcall split (double &hi, double &lo, double x) {
  double gamma = 0x8000001*x; // gamma=(2**27+1)*x; 27 = ceil(53/2), где 53 - размер мантиссы double (со скрытым битом).
  double delta = x - gamma;
  hi = gamma + delta;
  lo = x - hi;
}

    // Произведение с разбиением, результат идентичен функции two_prod
double __fastcall two_prod_s (double &t, double a, double b) {
  double ah, al, bh, bl, p;
  split (ah, al, a);
  split (bh, bl, b);
  p = a*b;
  t = (ah*bh-p) + ah*bl + al*bh + al*bl;
  return p;
}

Данный алгоритм для двоичных типов данных из Стандарта IEEE-754 работает всегда, для десятичных типов он будет работать только если количество цифр мантиссы чётно, то для есть для decimal32 он может давать ошибочный результат. Если же вы изобретаете свой тип данных с плавающей запятой, читайте полное условие теоремы 4.8 в упомянутой выше книге в разделе 4.4.2.2 «Dekker«s multiplication algorithm». Там же, но разделом ранее разъясняется как подбирать магическую константу для правильного разбиения в функции split. В двоичной системе она подбирается равной 2s+1, при этом младшая часть lo будет содержать s битов, а старшая — все остальные.
Итак, теперь мы умеем умножать два числа точно. Стряпаем сумму произведений по образу и подобию предыдущей статьи. Данный алгоритм назван именами трёх авторов Ogita—Rump—Oishi:

    // Функция точного сложения двух чисел из предыдущей статьи
double __fastcall two_sum (double &t, double a, double b) {
  double s = a+b;
  double bs = s-a;
  double as = s-bs;
  t = (b-bs) + (a-as);
  return s;
}

using dvect_t = vector;

double __fastcall ogita (const dvect_t &X, const dvect_t &Y) {
  double s=0.0, c=0.0, p, pi, t;
  for (size_t i=0; i

double __fastcall naive (const dvect_t &X, const dvect_t &Y) {
  double s=0.0;
  for (size_t i=0; i

double __fastcall naive (const dvect_t &A, double x) {
  double s=0.0;
  for (auto a: A) {
    s *= x;
    s += a;
  }
  return s;
}

double __fastcall naive_fma (const dvect_t &A, double x) {
  double s=0.0;
  for (auto a: A)
    s = fma (s, x, a);
  return s;
}

    // Алгоритм по именам трёх авторов: Graillat–Langlois–Louvet
double __fastcall graillat (const dvect_t &A, double x) {
  double s=0.0, c=0.0, p, pi, t;
  for (auto a: A) {
    // Полином из основных значений
    p = two_prod (pi, s, x);
    s = two_sum (t, p, a);
    // Полином из «хвостиков»
    c *= x;
    c += pi+t;  // Можете применить здесь FMA, но не вижу особого смысла
  }
  return s+c;	// Фактически, здесь сумма значений двух полиномов.
}