Итак, в прошлый раз я представил базовую идею как можно реализовать Fixed-point arithmetic, а так же набросок кода на C++, в котором в комментариях нашли довольно много проблем (а я сам нашёл ещё больше). В этот раз хочется представить улучшенную реализацию, разбор тонких моментов в коде и провести более детальный анализ получаемых результатов.
Инициализация из double
В прошлый раз, в конструкторе explicit FixedPoint(double v)у меня были довольно наивные строки для работы с дробной частью:
uint32_t fraction_decimм al = std::abs(fraction)* ipow(10, FractionLength + 1); uint8_t count = fraction_decimal / fraction_unit_halv;Во-первых, при степень 10ки растёт довольно быстро и может переполнить uint32_t, во-вторых…это просто не нужно: манстисса double – уже двоичное число, всё что нужно – сдвинуть его в область целых чисел, используя функцию std::scalbln, получаем:
UnsignedBasetype count = static_cast<UnsignedBasetype>(std::scalbln(fraction_abs, FractionLength + 1));Операторы умножения и деления
Пропущена важная деталь: умножать(делить) надо модули и потом навешивать знак:
FixedPoint operator * (const FixedPoint& r) const { auto total_sign = isign(val) * isign(r.val); HelperType temp = HelperType(std::abs(val)) * HelperType(std::abs(r.val)); temp >>= (FractionLength - 1); uint8_t unit = temp & least_bit_mask; return makeFx(total_sign * ((temp >> 1) + unit)); }Конвертирование в строку
При конвертировании в строку дробная часть преобразуется в безнаковое целое число (а именно uint64_t), путём умножения дробной части (опять же это просто целое безнаковое), на вес разряда. Например, 0.012 = 0.2510 , в данном случае вес разряда равен 0.25, в коде это обобщено как:
static constexpr uint64_t fraction_unit = ipow(five, FractionLength);К сожалению, тут имеет место тот же эффект, что был в конструкторе explicit FixedPoint(double v): степень 5ки растёт довольно быстро, а именно т.к. 5=2^2.231, то переполнение грозит уже начиная с FractionLength равный 15. Переполнение – это не тоже самое, что потеря точности – это вывод мусора вместо актуального результата. Пришлось найти workaround для этой проблемы:
if (fraction) { uint64_t current_fraction = fraction; uint64_t current_unit = fraction_unit; while ((std::numeric_limits<uint64_t>::max() / current_fraction) < current_unit) { uint8_t least_bit = current_fraction & least_bit_mask; current_unit /= 5; current_fraction = (current_fraction + least_bit)/ 2; } res << '.' << current_fraction * current_unit; }Таким образом, мы округляем нашу дробную часть пока она не влезет в uint64_t, что, конечно, вызывает некоторую потерю точности.
Шаблонные типы и обобщения
По предыдущим snippets, читатели уже могли заметить, что теперь вместо того что б хранить поле в int8_t и результаты промежуточных операций в uint_16t я перешёл к неким алиасам: UnsignedBasetype, HelperType и т.д. Давайте посмотрим всё (скромное) «метапрограммирование», имеющее место быть в коде:
// FractionLength - how many bits are after a period template <uint8_t FractionLength, typename Basetype = int32_t, typename HelperType = uint64_t> class FixedPoint { public: static_assert(sizeof(HelperType) == 2 * sizeof(Basetype)); static_assert(FractionLength > 0); static_assert(FractionLength < sizeof(Basetype) * CHAR_BIT); static_assert(std::is_integral<Basetype>::value); static_assert(std::is_integral<HelperType>::value); static_assert(std::is_signed<Basetype>::value); static_assert(std::is_unsigned<HelperType>::value); // I've cut out some code here - I'll publish it later private: Basetype val; using UnsignedBasetype = std::make_unsigned_t<Basetype>; // fills and return given number of ones in the least significant bits of a byte static constexpr UnsignedBasetype mask(uint8_t num) { return num == 0 ? 0 : (mask(num - 1) << 1) | 1; } static constexpr uint64_t ipow(uint8_t num, unsigned int pow) { return (pow >= sizeof(unsigned int)*8) ? 0 : pow == 0 ? 1 : num * ipow(num, pow-1); } template<typename T> static int8_t isign(T v) { return v >= 0? +1 : -1; } static constexpr uint8_t full_length = sizeof(Basetype) * CHAR_BIT; static constexpr UnsignedBasetype decimal_lenght = full_length - FractionLength; static constexpr UnsignedBasetype max_dec = (1 << (decimal_lenght - 1)) - 1; static constexpr Basetype min_dec_abs = (1 << (decimal_lenght - 1)); static constexpr Basetype min_dec = -min_dec_abs; static constexpr HelperType five = 5; static constexpr uint64_t fraction_unit = ipow(five, FractionLength); static constexpr UnsignedBasetype full_fraction = ipow(2, FractionLength); static constexpr UnsignedBasetype sign_mask = 1 << (full_length - 1); static constexpr UnsignedBasetype fraction_mask = mask(FractionLength); static constexpr UnsignedBasetype decimal_mask = mask(full_length) & ~FractionLength; static constexpr uint8_t least_bit_mask = 0x01; };По сути, сколько именно байтов будет в переменной val не так важно – важно, что б компилятор умел делать над этим типом арифметические действия, и что б при умножении мы не теряли разряды, т.е. вдвое больший по разрядности тип (HelperType) тоже существует. Идущие ниже asserts илюстрируют наши требования к этим типам. Типы по умолачнию – самые подходящие, как минимум на платформе x86_64.
Примеры
Если вы дочитали до этой точки, подозреваю, что вы хотите уже увидеть примеры:)) Тут важно понимать почему я так упорно инициализирую всё из double, хотя это и довольно сложный путь, одновременно он обладает очень хорошей наглядность: можно сразу сравить расчёт в double с расчётом с FixedPoint, но…на самом деле, самый первый расчёт есть сама инициализация из double, посему начнём с неё:
// let's introduce aliases for convenience using FixedPoint_2 = FixedPoint<2>; using FixedPoint_4 = FixedPoint<4>; using FixedPoint_10 = FixedPoint<10>; using FixedPoint_12 = FixedPoint<12>; using FixedPoint_14 = FixedPoint<14>; using FixedPoint_16 = FixedPoint<16>; using FixedPoint_18 = FixedPoint<18>; using FixedPoint_19 = FixedPoint<18>; using FixedPoint_20 = FixedPoint<20>; using FixedPoint_24 = FixedPoint<24>; using FixedPoint_28 = FixedPoint<28>; using FixedPoint_30 = FixedPoint<30>; double x = 0.261799; std::cout << "Original number: " << x << std::endl; std::cout << "Its fixed point versions: " << FixedPoint_2{x} << ' ' << FixedPoint_4{x} << ' ' << FixedPoint_16{x} << ' ' << FixedPoint_20{x} << ' ' << FixedPoint_24{x} << ' ' << FixedPoint_30{x} << std::endl;Выведет:
Original number: 0.261799
Its fixed point versions: 0.25 0.2500 0.2617950439453125 0.2617988586425781250 0.2618007659912109375 0.2525831760799073280
Давайте проанализируем результаты: изначально, при увеличении числа разрядов в дробной части точность растёт, достигая своего максимума при числе разрядов равным 20, а именно погрешность будет в районе одной миллионной (10^(-6)), а вот далее, точность почему-то падает. Это – небольшой сюрприз, точность падает из-за огрубления результатов в operator std::string() const. На самом деле, уже при 20 цифрах она влияет: теоретическая погрешность должна была быть порядка 2^(-21), т.е. примерно 5*10^(-7).
Посмотрим результаты работы при вычислении содержательных математических функций, а именно, возьмем несколько первых членов ряда Тейлора (Маклорена) для sin(x), exp(x):
template <typename T> T sinus(T x, uint8_t n = 3) { T t = x; T res = t; for ( int i=1; i<n; ++i) { T mult = -x*x/((2*i+1)*(2*i)); t *= mult; res += t; } return res; } template <typename T> T exponenta(T x, uint8_t n = 5) { T res{1.0}; T num {1.0}; uint64_t den {1}; for ( int i=1; i<n; ++i) { num *= x; den *= i; res += num/den; } return res; } // somewhere in the main function double x = 0.261799; FixedPoint_16 fx{x}; std::cout << std::setprecision(20) << exponenta(x) << ' ' << exponenta(fx) << ' ' << std::exp(x) << std::endl; std::cout << std::setprecision(20) << sinus(x) << ' ' << sinus(fx) << ' ' << std::sin(x) << std::endl;Выведет:
1.2992546509215898709 1.2992553710937500 1.2992653638436806318
0.25881868722557693774 0.2588195800781250 0.25881867051728746354
Первое значение – результат нашей апроксимации над обычным double, второе – с использованием нашего класса, и третье – результат реализации из стандартной библиотеки. Третье значение нас интересует меньше всего: алгоритм там отличается от нашего наивного ряда (можно посмотреть в исходники), хотя в обеих случаях наша апроксимакция близка к библиотечному референсу, но сосредоточимся на сравнении первых двух значений.
В обеих случаях разница между вычислениями над doube и FixedPoint имеет величину порядка 10^(-6) – это очень хороший показатель, почему? Потому что, имея под дробную часть 16 разрядов, мы имеем ошибку округления порядка 2^(-17) примерно 7*10^(-7).
Выводы
Эта реализация имеет неопровержимые плюсы:
-
компилируется и неплохо работает
-
построена из «первых принципов»
-
не имеет никаких зависимостей, кроме стандартной библиотеки
-
занимает всего ~250 строк, из которых содержательна примерно половина
И минусы тоже:
-
мало тестов
-
при большом количестве знаков после запятой теряется точность при преобразовании в строку
Код
#include <iomanip> #include <stdexcept> #include <cmath> #include <climits> #include <iostream> #include <sstream> // FractionLength - how many bits are after a period template <uint8_t FractionLength, typename Basetype = int32_t, typename HelperType = uint64_t> class FixedPoint { public: static_assert(sizeof(HelperType) == 2 * sizeof(Basetype)); static_assert(FractionLength > 0); static_assert(FractionLength < sizeof(Basetype) * CHAR_BIT); static_assert(std::is_integral<Basetype>::value); static_assert(std::is_integral<HelperType>::value); static_assert(std::is_signed<Basetype>::value); static_assert(std::is_unsigned<HelperType>::value); explicit FixedPoint(int decimal, unsigned int fraction = 0): val(0) { unsigned int decimal_abs = std::abs(decimal); if (decimal_abs > max_dec || decimal < min_dec || (decimal == min_dec && fraction) || (fraction >= full_fraction)) { throw std::invalid_argument("It won't fit our so few bits"); } val |= fraction; val |= (decimal_abs << FractionLength); val *= isign(decimal); } explicit FixedPoint(double v): val(0) { double decimal_double = 0.0; double fraction = std::modf(v, &decimal_double); if (decimal_double > max_dec || decimal_double < min_dec) { throw std::invalid_argument("It won't fit our so few bits"); } Basetype decimal_part = static_cast<Basetype>(decimal_double); int8_t sign = isign(v); UnsignedBasetype decimal_part_abs = std::abs(decimal_part); double fraction_abs = std::abs(fraction); UnsignedBasetype count = static_cast<UnsignedBasetype>(std::scalbln(fraction_abs, FractionLength + 1)); count += count & least_bit_mask; count >>= 1; if (count && decimal_part == min_dec_abs) { throw std::invalid_argument("It won't fit our so few bits"); } if (count == full_fraction) { decimal_part_abs += 1; auto decimal_part_signed = sign * decimal_part; if (decimal_part_abs > max_dec || decimal_part_signed < min_dec) { throw std::invalid_argument("It won't fit our so few bits"); } } else { val |= count; } val |= (decimal_part_abs << FractionLength); val *= sign; } FixedPoint operator + (const FixedPoint& r) const { return makeFx(val + r.val); } FixedPoint operator - (const FixedPoint& r) const { return makeFx(val - r.val); } FixedPoint operator - () const { return makeFx(-val); } FixedPoint operator * (const FixedPoint& r) const { auto total_sign = isign(val) * isign(r.val); HelperType temp = HelperType(std::abs(val)) * HelperType(std::abs(r.val)); temp >>= (FractionLength - 1); uint8_t unit = temp & least_bit_mask; return makeFx(total_sign * ((temp >> 1) + unit)); } FixedPoint operator / (const FixedPoint& r) const { auto total_sign = isign(val) * isign(r.val); HelperType left = std::abs(val); left <<= (FractionLength + 1); UnsignedBasetype temp = left / std::abs(r.val); uint8_t unit = temp & least_bit_mask; return makeFx(total_sign * ((temp >> 1) + unit )); } operator std::string() const { std::stringstream res; UnsignedBasetype temp = val; auto fraction = temp & fraction_mask; if (sign_mask & temp) { res << '-'; temp = ~temp + 1; fraction = temp & fraction_mask; } res << (temp >> FractionLength); if (fraction) { uint64_t current_fraction = fraction; uint64_t current_unit = fraction_unit; while ((std::numeric_limits<uint64_t>::max() / current_fraction) < current_unit) { uint8_t least_bit = current_fraction & least_bit_mask; current_unit /= 5; current_fraction = (current_fraction + least_bit)/ 2; } res << '.' << current_fraction * current_unit; } return res.str(); } private: Basetype val; using UnsignedBasetype = std::make_unsigned_t<Basetype>; // fills and return given number of ones in the least significant bits of a byte static constexpr UnsignedBasetype mask(uint8_t num) { return num == 0 ? 0 : (mask(num - 1) << 1) | 1; } static constexpr uint64_t ipow(uint8_t num, unsigned int pow) { return (pow >= sizeof(unsigned int)*8) ? 0 : pow == 0 ? 1 : num * ipow(num, pow-1); } template<typename T> static int8_t isign(T v) { return v >= 0? +1 : -1; } static constexpr uint8_t full_length = sizeof(Basetype) * CHAR_BIT; static constexpr UnsignedBasetype decimal_lenght = full_length - FractionLength; static constexpr UnsignedBasetype max_dec = (1 << (decimal_lenght - 1)) - 1; static constexpr Basetype min_dec_abs = (1 << (decimal_lenght - 1)); static constexpr Basetype min_dec = -min_dec_abs; static constexpr HelperType five = 5; static constexpr uint64_t fraction_unit = ipow(five, FractionLength); static constexpr UnsignedBasetype full_fraction = ipow(2, FractionLength); static constexpr UnsignedBasetype sign_mask = 1 << (full_length - 1); static constexpr UnsignedBasetype fraction_mask = mask(FractionLength); static constexpr UnsignedBasetype decimal_mask = mask(full_length) & ~FractionLength; static constexpr uint8_t least_bit_mask = 0x01; explicit FixedPoint(): val(0) { } static FixedPoint makeFx(Basetype v) { FixedPoint fp; fp.val = v; return fp; } }; template <uint8_t FractionLength> std::ostream& operator << (std::ostream& out, const FixedPoint<FractionLength>& number) { out << std::string(number); return out; } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator + (const FixedPoint<FractionLength>& l, int r) { return l + FixedPoint<FractionLength>(r); } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator + (int l, const FixedPoint<FractionLength>& r) { return r + FixedPoint<FractionLength>(l); } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator - (const FixedPoint<FractionLength>& l, int r) { return l - FixedPoint<FractionLength>(r); } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator - (int l, const FixedPoint<FractionLength>& r) { return FixedPoint<FractionLength>(l) - r; } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator * (const FixedPoint<FractionLength>& l, int r) { return l * FixedPoint<FractionLength>(r); } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator * (int l, const FixedPoint<FractionLength>& r) { return r * FixedPoint<FractionLength>(l); } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator / (const FixedPoint<FractionLength>& l, int r) { return l / FixedPoint<FractionLength>(r); } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator *= (FixedPoint<FractionLength>& l, const FixedPoint<FractionLength>& r) { l = l * r; return l; } template <uint8_t FractionLength> const FixedPoint<FractionLength> operator += (FixedPoint<FractionLength>& l, const FixedPoint<FractionLength>& r) { l = l + r; return l; } // let's introduce aliases for convenience using FixedPoint_2 = FixedPoint<2>; using FixedPoint_4 = FixedPoint<4>; using FixedPoint_10 = FixedPoint<10>; using FixedPoint_12 = FixedPoint<12>; using FixedPoint_14 = FixedPoint<14>; using FixedPoint_16 = FixedPoint<16>; using FixedPoint_18 = FixedPoint<18>; using FixedPoint_19 = FixedPoint<18>; using FixedPoint_20 = FixedPoint<20>; using FixedPoint_24 = FixedPoint<24>; using FixedPoint_28 = FixedPoint<28>; using FixedPoint_30 = FixedPoint<30>; PS
Считаю, что тема полностью раскрыта, но если вам есть, что тут предложить – пишите в комментарии.
ссылка на оригинал статьи https://habr.com/ru/articles/832258/
Добавить комментарий