![](https://habrastorage.org/getpro/habr/upload_files/ed2/ed9/f1d/ed2ed9f1df411ec4e438528e1acca1d3.png)
Итак, в прошлый раз я представил базовую идею как можно реализовать 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
Считаю, что тема полностью раскрыта, но если вам есть, что тут предложить – пишите в комментарии.