Итак, в прошлый раз я представил базовую идею как можно реализовать 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.01₂ = 0.25₁₀ , в данном случае вес разряда равен 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

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