Точные и быстрые вычисления для чисел с плавающей точкой на примере функции синуса. Часть 3: fixed-point

    Продолжаем цикл лекций (часть 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 оборудован последней версией математического сопроцессора. Если бы процессор эмулировал бы вычисления с плавающей точкой, то результаты могли бы разительно различаться.
    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 14

      +1
      Кто сказал, что значения tsx должны быть равны точно коэффициентам ряда Тейлора. Напрашивается достаточно очевидная идея проварьировать значения коэффициентов, чтобы улучшить точность аппроксимации и убрать постоянную составляющую ошибки. Грамотно сделать такую вариацию достаточно сложно.
      А если воспользоваться математическим пакетом — то очень даже просто, потому что там обычно присутствует инструментал для аппроксимации. В WM, в частности, есть функция FindFit, которой можно аппроксимировать что угодно чем угодно (в том числе и рациональным полиномом, который для многих функций даёт большую точность).
        0
        всего 10 используеых коэффициентов, у каждого можно поменять только три-четыре младших бита максимум, полный перебор потребует 4^10- миллион вариантов. Примерно 2е суток на машине автора перебирать пришлось бы. Если не трогать младшие коэффициенты (например, до пятого включительно), то весь перебор ужмется до трех минут.
          0
          Я перебирал 30 коэффициентов и Вольфрам справлялся за несколько секунд, не минут. Методом наименьших квадратов так и ещё быстрее.
            +1
            Мне кажется, вы спутали МНК и градиентный спуск (МНК — это не про поиск максимумов и минимумов функций многих переменных), это раз, и второе- у нас целевая функция- это средняя ошибка округления на большом множестве чисел- она сама по себе вычисляется долго (ТС указал на х86- 0,174с на вычисление одного значения- быстро перебирать варианты не получится), при этом коэффициенты- де-факто целочисленные, и потому к ним градиентный спуск не применим (градиента для целых чисел нету), мы уже имеем решение, близкое к идеальному (где градиентный спуск неэффективен), это решение очень близкое к идеальному- отличающесеся в одном-двух младших битах, поэтому мы даже не сможем вычислить вектор оптимального приращения, и решение возле текущей аппроксимации осциллирует- и градиентный и покоординатный спуски на таких функциях ломаются. Вот и остается чистый перебор множетства значений медленно вычисляемой функции.
              0
              Спасибо за комментарий! Всё абсолютно правильно. Более того добавлю, что метод наименьших квадратов использовать нельзя. Нужно минимизировать максимальное значение отклонения. А градиент от этого посчитать ещё более невозможно. В следующих статьях я покажу, как можно аналитически аппроксимировать последний коэффициент. Последние коэффициенты, очевидно можно варьировать больше, чем последние биты. Там есть красивый математических подход.
                +1
                Мне кажется, вы спутали МНК и градиентный спуск
                Мне так не кажется, потому что я и не утверждал, что они дают идентичный результат.

                коэффициенты- де-факто целочисленные, и потому к ним градиентный спуск не применим (градиента для целых чисел нету)
                Так коэффициенты можно округлять уже после, а ещё лучше — рационализировать. Понятно, что таким образом (равно как и МНК) наилучшего результата не получить — но в качестве первичной аппроксимации сойти вполне может, и результат таки может оказаться лучше Тейлора.
          0
          «Всё-таки long double достаточно экзотический вид данных и его поддержка в современных процессорах не является приоритетной».
          А если в процессоре нет целочисленного умножения 64*64? :) Можно конечно через 32*32,
          но это в 10 раз медленнее чем double*double.
          Правильно ли я понимаю, что если в сложении-вычитании double можно производить вычисления большей точности за счет double.double, то с умножением такой фокус не проходит? Умножили и хвост безнадежно потерян. Кстати, в процессоре можно легко сделать вычисления double+double c сохранением обоих частей результата. Аналогично можно было бы сделать и с умножением. т.е. по аналогии с целыми числами int*int = long long.
          Был бы какой-то толк от подобных инструкций?
            0
            1) А если в процессоре нет аппаратной поддержки FP арифметики? Почему вы говорите, что
            но это в 10 раз медленнее чем double*double.

            Если у Вас есть какой-то конкретный процессор возьмите его и сделайте эксперимент.
            2) a) Есть такая суперштука как FMA. С помощью этих комманд можно проводить точное умножение.
            b) Если мантисса заполнена на половину, то результат умножения будет точным. Поэтому использование битовых масок решит проблему. Правда понадобится много комманд.
              0
              если в сложении-вычитании double можно производить вычисления большей точности за счет double.double, то с умножением такой фокус не проходит?
              Проходит, есть алгоритмы для умножения и double на double, и double-double на double-double. Тут в архиве qd-2.3.22.tar.gz есть и описание в pdf, и реализация на си и фортране.
                0
                Спасибо. Именно об этом и был мой вопрос. Попробую посмотреть, может чего и пойму :)
                А то жалко иметь хороший аппаратный double, а расчеты производить на счетах, потому что точнее выходит :)
              0
              1. я как раз про свой конкретный случай :) в моем процессоре можно 2 double умножения за такт, но ни одного int64*int64. Эксперимент с sin_e7 примером я провел. Реализовал команду умножения int64*int64 через последовательность умножений int32*int32 и сложений.
              2. насколько я помню, FMA избегает ненужного округления перед сложением. Результат будет немного точнее. Но у меня нет и FMA :) Только умножить и сложить.
              3. Да, я использовал такой подход в синусе. Тот еще геморрой :)
                0
                Если Вам интересна практическая реализация то почему бы Вам не адаптировать s_sin? Вам нужна большая точность? Если так то самым простым образом можно решить эту проблему увеличением числа членов в разложении синуса и косинуса:
                ...
                  s = x + (dx + x * xx * (sn3 + xx * sn5)); // здесь
                  c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); // и здесь
                  SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
                ...
                


                  0
                  P.S. По моим прикидкам следующего члена будет достаточно (cs7 и cs8) соответственно. Вопрос, разумеется остаётся к качеству самой таблички. Я, если честно, не проверял.
                  P.P.S. Проверьте существующие значения cs и sn. Они могут отличатся от ряда Тейлора. Если так, то рассчитайте их сами точно.
                    0
                    Спасибо. C s_sin я уже активно работаю :) У меня имеется своя реализация sin-cos (некоторые подробности реализации я изложил в части 2). Интерес к s_sin у меня как к генератору эталонного результата. Также в s_sin мне понравились некоторые приемы работы с double. Моя реализация очень мала по обьему кода и по размеру таблицы (80 байт) и имеет очень хорошую скорость. Длина полинома у меня на два члена длиннее чем в s-sin. Цель — улучшить точность, незначительно потеряв в скорости. Вполне возможно, что я реализую на ассемблере и s_sin, как альтернативный вариант в библиотеке.

              Only users with full accounts can post comments. Log in, please.