diyfp.hpp

Go to the documentation of this file.

 
#ifndef JOIN_DATA_DIYFP_HPP
#define JOIN_DATA_DIYFP_HPP
 
// C++.
#include <limits>
 
// C.
#include <cstdint>
#include <cstring>
 
namespace join
{
    class DiyFp
    {
    public:
        constexpr DiyFp () noexcept = default;
 
        constexpr DiyFp (const DiyFp& other) noexcept = default;
 
        constexpr DiyFp& operator= (const DiyFp& other) noexcept = default;
 
        explicit DiyFp (double value) noexcept
        {
            uint64_t u64 = 0;
            memcpy (&u64, &value, sizeof (double));
 
            _mantissa = u64 & _mantissaMask;
            _exponent = static_cast<int> ((u64 & _exponentMask) >> _mantissaSize);
 
            if (_exponent)
            {
                _mantissa += _hiddenBit;
                _exponent -= _exponentBias;
            }
            else
            {
                _exponent = 1 - _exponentBias;
            }
        }
 
        constexpr DiyFp (uint64_t mantissa, int exponent) noexcept
        : _mantissa (mantissa)
        , _exponent (exponent)
        {
        }
 
        ~DiyFp () = default;
 
        inline constexpr DiyFp& normalize () noexcept
        {
            if (_mantissa == 0)
            {
                return *this;
            }
 
            int shift = __builtin_clzll (_mantissa);
            _mantissa <<= shift;
            _exponent -= shift;
 
            return *this;
        }
 
    private:
        inline constexpr DiyFp& normalizeBoundary () noexcept
        {
            if (_mantissa != 0)
            {
                int shift = __builtin_clzll (_mantissa) - (64 - _mantissaSize - 2);
                if (shift > 0)
                {
                    _mantissa <<= shift;
                    _exponent -= shift;
                }
            }
 
            constexpr int shift = _diyMantissaSize - _mantissaSize - 2;
            _mantissa <<= shift;
            _exponent -= shift;
 
            return *this;
        }
 
    public:
        constexpr void normalizedBoundaries (DiyFp& minus, DiyFp& plus) const noexcept
        {
            plus._mantissa = (_mantissa << 1) + 1;
            plus._exponent = _exponent - 1;
            plus.normalizeBoundary ();
 
            const bool special = __builtin_expect (_mantissa == _hiddenBit, 0);
            minus._mantissa = (_mantissa << (special ? 2 : 1)) - 1;
            minus._exponent = _exponent - (special ? 2 : 1);
 
            const int diff = minus._exponent - plus._exponent;
            minus._mantissa <<= diff;
            minus._exponent = plus._exponent;
        }
 
        inline constexpr DiyFp& operator-= (const DiyFp& rhs) noexcept
        {
            _mantissa -= rhs._mantissa;
            return *this;
        }
 
        inline constexpr DiyFp& operator*= (const DiyFp& rhs) noexcept
        {
#if defined(__SIZEOF_INT128__)
            __uint128_t product = static_cast<__uint128_t> (_mantissa) * static_cast<__uint128_t> (rhs._mantissa);
            _mantissa = static_cast<uint64_t> ((product >> 64) + ((product >> 63) & 1));
#else
            uint64_t M32 = 0xFFFFFFFFU;
 
            uint64_t a = _mantissa >> 32;
            uint64_t b = _mantissa & M32;
            uint64_t c = rhs._mantissa >> 32;
            uint64_t d = rhs._mantissa & M32;
 
            uint64_t ac = a * c;
            uint64_t bc = b * c;
            uint64_t ad = a * d;
            uint64_t bd = b * d;
 
            uint64_t tmp = (bd >> 32) + (ad & M32) + (bc & M32) + (1U << 31);
            _mantissa = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32);
#endif
            _exponent += rhs._exponent + 64;
 
            return *this;
        }
 
        static constexpr int _diyMantissaSize = std::numeric_limits<uint64_t>::digits;
 
        static constexpr int _mantissaSize = std::numeric_limits<double>::digits - 1;
 
        static constexpr int _exponentBias = 0x3FF + _mantissaSize;
 
        static constexpr uint64_t _mantissaMask = 0x000FFFFFFFFFFFFFLLU;
 
        static constexpr uint64_t _exponentMask = 0x7FF0000000000000LLU;
 
        static constexpr uint64_t _hiddenBit = 0x0010000000000000LLU;
 
        uint64_t _mantissa = 0;
 
        int _exponent = 0;
    };
 
    inline constexpr DiyFp operator- (const DiyFp& lhs, const DiyFp& rhs)
    {
        return DiyFp (lhs) -= rhs;
    }
 
    inline constexpr DiyFp operator* (const DiyFp& lhs, const DiyFp& rhs)
    {
        return DiyFp (lhs) *= rhs;
    }
}
 
#endif