Имплементация чисел с фиксированной точкой (часть 2)29.07.2024 16:15

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(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::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 
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::value);
    static_assert(std::is_integral::value);
    static_assert(std::is_signed::value);
    static_assert(std::is_unsigned::value);
  
// I've cut out some code here - I'll publish it later 
private:
    Basetype val;

    using UnsignedBasetype = std::make_unsigned_t;

    // 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
    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 
T sinus(T x, uint8_t n = 3)
{
    T t = x;
    T res = t;
    for ( int i=1; i
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

Выведет:

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 
#include 
#include 
#include 
#include 
#include 

// FractionLength - how many bits are after a period
template 
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::value);
    static_assert(std::is_integral::value);
    static_assert(std::is_signed::value);
    static_assert(std::is_unsigned::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(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(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::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;

    // 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
    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 
std::ostream& operator << (std::ostream& out, const FixedPoint& number)
{
    out << std::string(number);

    return out;
}

template 
const FixedPoint operator + (const FixedPoint& l, int r)
{
    return l + FixedPoint(r);
}

template 
const FixedPoint operator + (int l, const FixedPoint& r)
{
    return r + FixedPoint(l);
}

template 
const FixedPoint operator - (const FixedPoint& l, int r)
{
    return l - FixedPoint(r);
}

template 
const FixedPoint operator - (int l, const FixedPoint& r)
{
    return FixedPoint(l) - r;
}

template 
const FixedPoint operator * (const FixedPoint& l, int r)
{
    return l * FixedPoint(r);
}

template 
const FixedPoint operator * (int l, const FixedPoint& r)
{
    return r * FixedPoint(l);
}

template 
const FixedPoint operator / (const FixedPoint& l, int r)
{
    return l / FixedPoint(r);
}

template 
const FixedPoint operator *= (FixedPoint& l, const FixedPoint& r)
{
    l = l * r;
    return l;
}

template 
const FixedPoint operator += (FixedPoint& l, const FixedPoint& 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

Считаю, что тема полностью раскрыта, но если вам есть, что тут предложить — пишите в комментарии.

© Habrahabr.ru