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


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


Умножение двух чисел без потери точности (ну… почти)


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

Обозначим через $e_{min}$ минимальную экспоненту вашей системы чисел с плавающей запятой, в рамках которой выполняется умножение (для float: ?126, для double: ?1022), буквой $k$ обозначим количество битов мантиссы с учётом скрытой единицы (24 для float и 53 для double). Пусть $e_a$ и $e_b$ — экспоненты наших исходных чисел $a$ и $b$, тогда требуется, чтобы было выполнено условие:

$e_a+e_b \geq e_{min} +k-1.$


Если данное условие не выполняется, то могут возникнуть случаи, когда $ab$ нельзя представить в виде суммы $p+t$, где $p = RN(ab)$ и $t=ab-p$. Указанное неравенство является сложным способом сказать простую фразу: «хвостик» $t$ не должен быть денормализованным. И это вполне понятно: когда мы перемножаем два числа, получаем удвоенное количество битов, которые гарантированно влезают в две переменные того же типа данных. Если же одна из этих переменных попадает в денормализованный диапазон, происходит битовый сдвиг мантиссы с потерей младших битов. И кроме переполнения это единственная ситуация, при которой произведение двух чисел нельзя точно представить в форме суммы двух чисел того же типа.

Так ли существенно это ограничение? Для 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, не забудьте подключить <cmath>
  return p;
}

Результатом работы функции будет пара чисел $(p, t)$, такая, что $p$ является ближайшим к $ab$ числом, а $t$ — хвостиком-погрешностью, причём $p+t=ab$ в точности. Разумеется, по ходу расчётов не должно произойти переполнения, за этим придётся следить отдельно: нужно оставаться в рамках условия $RN(ab)\neq\infty$.

Если операция FMA для вас недоступна, то ситуация сильно усложняется. Желаемое произведение можно получить путём разбиения исходных переменных $a$ и $b$ на две части примерно посередине так, чтобы старшая часть операндов (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>;

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<X.size(); ++i) {
    p = two_prod (pi, X[i], Y[i]);	// p = X[i]*Y[i], но в pi остался «хвостик».
    s = two_sum (t, s, p);	// s += p, но в t остался «хвостик».
    c += pi+t;	// Собираем «хвостики»
  }
  return s+c;	// Общая сумма и «хвостик»
}

Понятно, что two_prod нужно заменить на two_prod_s, если вам недоступна FMA.

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

Помимо этой версии скалярного произведения есть ещё две наивные функции, разъяснение которых в тексте статьи будет моветоном.

double __fastcall naive (const dvect_t &X, const dvect_t &Y) {
  double s=0.0;
  for (size_t i=0; i<X.size(); ++i)
    s += X[i]*Y[i];
  return s;
}

double __fastcall naive_fma (const dvect_t &X, const dvect_t &Y) {
  double s=0.0;
  for (size_t i=0; i<X.size(); ++i)
    s = fma (X[i], Y[i], s);
  return s;
}

Тестирование выполнялось точно также как в статье про сумму N чисел. Напомню только, что верхнее число в ячейке — средняя погрешность в ULP на 100 тестах, а нижнее — максимальная погрешность на них же. Размер каждого теста N=106. В первом столбце указано выбранное распределение случайных чисел.

Naive Naive FMA Ogita—Rump—Oishi
U[1, 2) 95.31
275		 95.36
276		 0.00
0
±U[1, 2) 441.61
6878		 432.99
7489		 0.00
0
U[1e-10, 1e10) 2478.61
2740		 2478.55
2739		 0.00
0
±U[1e-10, 1e10) 392.50
11893		 394.86
12126		 0.00
0
exp[2] 183.70
513		 183.65
512		 0.00
0
±exp[2] 329.21
5962		 347.31
7251		 0.00
0
N(0, 1) 1554.07
109700		 1408.30
95896		 0.00
0

Ожидаемо, что алгоритм с применением double-double сработал безошибочно, и так будет происходить на всех типичных бытовых задачах. Таблица не требует особых комментариев кроме одного: наивный алгоритм с FMA может давать менее точный результат, чем ещё более наивный алгоритм с двумя округлениями на каждом шаге. При этом ДА, я выполнил дизассемблирование и убедился, что в этом наивном методе компилятор не всунул FMA в качестве оптимизации. Вообще говоря, он не имеет на это права без моего разрешения, но всё же осторожность не помешает.

Вычисление значения полинома


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

Когда мы применяем схему Горнера, то по сути имеем похожий по структуре цикл, в теле которого также одно умножение и сложение, поэтому никаких дополнительных пояснений к построению алгоритма не требуется, просто пишем код. Предполагается, что коэффициенты полинома записаны от старшего к младшему, а степень полинома на единицу меньше количества элементов в массиве A.

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;	// Фактически, здесь сумма значений двух полиномов.
}

Поскольку значение полинома на случайных тестах очень быстро возрастает и создаёт трудности для точного расчёта (эталонное значение ответа, как вы знаете, я вычисляю через дроби), для тестов взято ограничение N=100 элементов в массиве. Комментарии снова излишни, кроме одного: здесь применение FMA для наивного расчёта становится оправданным — точность заметно повышается.

Naive Naive FMA Graillat—Langlois—Louvet
U[1, 2) 3.15
11		 2.67
9		 0.00
0
±U[1, 2) 3.96
29		 3.03
27		 0.00
0
U[1/10, 10) 1.86
11		 1.47
7		 0.00
0
±U[1/10, 10) 2.20
61		 1.29
10		 0.00
0
exp[2] 0.59
5		 0.46
4		 0.00
0
±exp[2] 0.90
30		 0.71
6		 0.00
0
N(0, 1) 0.90
30		 0.71
6		 0.00
0

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

Одна из лучших наград автору — успешное применение описанных им знаний на практике для всеобщей пользы. Благодарю за внимание!

Список источников


[1] Разделы 4.4.2.2 и 5.4-5.5 книги Jean-Michel Muller, “Handbook of floating-point arithmetic”, 2018.

[2] Раздел 8.2 книги (по поводу умножения целых чисел) Генри С. Уоррен мл. «Алгоритмические трюки для программистов», 2014.