10.11.2020 00:17
dorex 14 1300 Источник
Продолжаем цикл лекций (часть 1 и часть 2). В части 2 мы посмотрели, что внутри у библиотеки libm и в данной работе попробуем немного переделать функцию do_sin, чтобы увеличить её точность и скорость работы. Приведу эту функцию ещё раз do_sin):

image

Как было показано в предыдущей статье часть 132-145. Выполняется для x в пределах [0.126, 0.855469]. Ну что. Попробуем написать функцию, которая в данных пределах будет точнее и, возможно, быстрее.

Способ, которым мы воспользуемся — достаточно очевидный. Надо расширить точность вычислений, чтобы включить больше знаков после запятой. Очевидным решением было бы выбрать тип long double, посчитать в нём, а потом сконвертировать обратно. С точки зрения точности решение должно быть хорошим, но с точки зрения производительности могут быть проблемы. Всё-таки long double достаточно экзотический вид данных и его поддержка в современных процессорах не является приоритетной. На x86_64 SSE/AVX инструкции с данным типом данных не работают. Будет «отдуваться» математический сопроцессор.

Что же тогда выбрать? Давайте посмотрим ещё раз внимательно пределы аргумента и функции.

Они находятся в районе 1.0. Т.е. по сути плавающая точка нам не нужна. Давайте при расчёте функции воспользуемся 64-битным целым. Это даст нам дополнительно 10-11 двоичных разрядов к исходной точности. Давайте разбираться, как же с этими числами работать. Число в данном формате представляется в виде a/d, где a это целое число, а d это делитель, который мы выбираем постоянным для всех переменных и храним «у себя в памяти», а не в памяти компьютера. Ниже приведены некоторые операции для таких чисел:

$\frac{c}{d}=\frac{a}{d}\pm\frac{b}{d}=\frac{a\pm b}{d}\\ \frac{c}{d}=\frac{a}{d}\cdot\frac{b}{d}=\frac{a \cdot b}{d^2}\\ \frac{c}{d}=\frac{a}{d}\cdot x=\frac{a \cdot x}{d}$


Как видите, ничего сложного в этом нет. Последняя формула показывает умножение на любое целое число. Заметим так же достаточно очевидную вещь, что результат умножения двух беззнаковых целых переменных размера N это чаще число размера до 2 * N включительно. Сложение же может вызвать переполнение максимум на 1 дополнительный бит.

Давайте попробуем выбрать делитель d. Очевидно, что в двоичном мире лучше всего его выбрать степенью двойки, чтобы не делить, а просто двигать регистр, например. Какую же степень двойки выбрать? Подсказку найдём в машинных инструкциях умножения. Например, стандартная инструкция MUL в системе x86 умножает 2 регистра и результат записывает тоже в 2 регистра, где 1 из регистров это «верхняя часть» результата, а второй — нижняя часть.

Например, если у нас есть 2 64-битных числа, то результатом будет 128-битное число, записанное в два 64-битных регистра. Назовём RH — «верхний» регистр, а RL — «нижний» регистр 1. Тогда, математически, результат можно записать в виде $R=R_H \cdot 2^{64}+R_L$. Теперь воспользуемся формулами сверху и напишем умножение для $d=2^{-64}$

$\frac{c}{d}=\frac{a}{2^{64}}\cdot\frac{b}{2^{64}}=\frac{a \cdot b}{2^{128}}=\frac{R_H \cdot 2^{64} + R_L}{2^{128}}=\frac{R_H + R_L \cdot 2^{-64}}{2^{64}}$


И получается, что результатом умножения этих двух чисел с фиксированной точкой является регистр $R=R_H$.

Для системы Aarch64 всё ещё проще. Команда «UMULH» умножает два регистра и записывает в 3-ий регистр «верхнюю» часть умножения.

Ну что же. Число с фиксированной точкой мы задали, но остаётся ещё одна проблема. Отрицательные числа. В ряду Тейлора разложение идёт с переменным знаком. Чтобы справится с этой проблемой преобразуем формулу расчёта полинома по методу Гонера, к следующему виду:

$\sin(x)\approx x(1-x^2(1/3!-x^2(1/5!-x^2(1/7!-x^2\cdot1/9!))))$


Проверьте, что она с математической точки зрения абсолютно такая же, как и исходная формула. Но в каждой скобке число вида $1/(2n + 1)! - x^2\cdot(\cdots)$ всегда положительно. Т.е. такое преобразование позволяет рассчитывать выражение в беззнаковых целых числах.

constexpr mynumber toint    = {{0x00000000, 0x43F00000}};  /*  18446744073709551616 = 2^64     */
constexpr mynumber todouble = {{0x00000000, 0x3BF00000}};  /*  ~5.42101086242752217003726400434E-20 = 2^-64     */

double sin_e7(double xd) {
  uint64_t x = xd * toint.x;
  uint64_t xx = mul2(x, x);
  uint64_t res = tsx[19]; 
  for(int i = 17; i >= 3; i -= 2) {
    res = tsx[i] - mul2(res, xx);
  }
  res = mul2(res, xx);
  res = x - mul2(x, res);
  return res * todouble.x;
}

Значения tsx[i] 
constexpr array<uint64_t, 18> tsx = { // 2^64/i!
    0x0000000000000000LL,
    0x0000000000000000LL,
    0x8000000000000000LL,
    0x2aaaaaaaaaaaaaaaLL, // Change to 0x2aaaaaaaaaaaaaafLL and check.
    0x0aaaaaaaaaaaaaaaLL,
    0x0222222222222222LL,
    0x005b05b05b05b05bLL,
    0x000d00d00d00d00dLL,
    0x0001a01a01a01a01LL,
    0x00002e3bc74aad8eLL,
    0x0000049f93edde27LL,
    0x0000006b99159fd5LL,
    0x00000008f76c77fcLL,
    0x00000000b092309dLL,
    0x000000000c9cba54LL,
    0x0000000000d73f9fLL,
    0x00000000000d73f9LL,
    0x000000000000ca96LL
};


Где $tsx[i]=1/i!$ в формате с фиксированной точкой. В этот раз для удобства я выложил весь код на гитхаб fast_sine, избавился от quadmath для совместимости с clang и arm. И поменял немного методику расчёта погрешности.

Данные по сравнению стандартной функции синуса и fixed point даны в двух таблицах ниже. В первой таблице приведена точность расчётов (она полностью совпадает для x86_64 и ARM). Во второй таблице — сравнение производительности.

Функция Количество ошибок Максимальное значение ULP Среднее значение отклонения
sin_e7 0.0822187% 0.504787 7.10578e-20
sin_e7a 0.0560688% 0.503336 2.0985e-20
std::sin 0.234681% 0.515376 ---


При проведение тестирования «истинное» значение синуса считалось с помощью библиотеки MPFR. Максимальное значение ULP считалось, как максимальное отклонение от «истинного» значения. Процент ошибок — количество случаев, когда рассчитанное значение функции синуса нами или libm бинарно не совпадало с округлённым до double значением синуса. Среднее значение отклонения показывает «направление» ошибки расчётов: завышение или занижение значения. Как видно из таблицы наша функция склонна завышать значения синуса. Это можно поправить! Кто сказал, что значения tsx должны быть равны точно коэффициентам ряда Тейлора. Напрашивается достаточно очевидная идея проварьировать значения коэффициентов, чтобы улучшить точность аппроксимации и убрать постоянную составляющую ошибки. Грамотно сделать такую вариацию достаточно сложно. Но попробовать мы можем. Давайте возьмем, например, 4-ое значение из массива коэффициентов tsx (tsx[3]) и поменяем последнее число a на f. Перезапустим программу и посмотрим на точность (sin_e7a). Посмотрите, она немного, но увеличилась! Добавляем в нашу копилку и этот способ.

Теперь посмотрим, что с производительностью. Для тестирования я взял что было под рукой i5 mobile и немного разогнанную четвёртую малинку (Raspberry PI 4 8GB), GCC10 из дистрибутива Ubuntu 20.04 x64 для обоих систем.

Функция x86_64 время, c ARM время, c
sin_e7 0.174371 0.469210
std::sin 0.154805 0.447807


Я не претендую на большую точность данных измерений. Возможны вариации в несколько десятков процентов в зависимости от загруженности процессора. Основной вывод сделать можно и так. Переход на целочисленную арифметику не даёт выигрыш в производительности на современных процессорах2. Невообразимое количество транзисторов в современных процессорах делает возможным быстро делать сложные вычисления. Но, я думаю, что на таких процессорах, как Intel Atom, а так же на слабых контроллерах данный подход может дать существенный прирост производительности. А вы что думаете?

Хотя данный подход дал прирост точности, всё таки этот прирост точности кажется более любопытным, чем полезным. По производительности данный подход может себя найти, например, в IoT. Но для высокопроизводительных вычислений он уже совсем не маинстрим. В современном мире SSE/AVX/CUDA предпочитают использовать параллельный расчёт функций. Причем в арифметике с плавающей точкой. Параллельных аналогов функции MUL нет. Сама функция это уже скорее дань традиции.

В следующей главе я опишу, как можно эффективно использовать AVX для расчётов. Опять же залезем в код libm и попробуем его улучшить.

1 Регистров с такими названиями в известных мне процессорах нет. Имена были выбраны, для примера.
2 Тут надо отметить, что мой ARM оборудован последней версией математического сопроцессора. Если бы процессор эмулировал бы вычисления с плавающей точкой, то результаты могли бы разительно различаться.