浮動小数点数をテンプレートメタプログラミングで扱う

IEEE 754に準拠している倍精度浮動小数点数は符号が1ビット、指数部が11ビット、仮数部が1+52ビットという構造になっています。

\begin{eqnarray}
x&=&(-1)^s\times 2^{e}\times 1.a_1a_2...a_{52}\\
&=&(-1)^s\times 2^{e}\times 1a_1a_2...a_{52}.0\times 2^{-52}\\
&=&(-1)^s\times 2^{e-52}\times 1a_1a_2...a_{52}.0
\end{eqnarray}

このように仮数部が64ビット以下の整数になるように指数部を調整して、元の値をstatic constexprな値として保持するようにすると、double型をテンプレートメタプログラミングで扱うことができます。

main.cpp

#include <iostream>

#include <cstdint>
#include <limits>
#include <tuple>
#include <type_traits>

#include <sprout/math/isnan.hpp>
#include <sprout/math/signbit.hpp>
#include <sprout/math/float2_exponent.hpp>
#include <sprout/math/isnormal.hpp>
#include <sprout/math/fabs.hpp>

namespace a{
    constexpr double pow2(int n)
    {
        double x = 1.0;
    
        if(n > 0){
            while(n-- != 0)
                x *= 2.0;
        }else if(n < 0){
            while(n++ != 0)
                x /= 2.0;
        }
    
        return x;
    }
    
    template <bool Negative, ::std::uint64_t Significand, int Exponent>
    struct double_{
        static constexpr double value =
            (Negative ? -1.0 : 1.0) * static_cast<double>(Significand) * pow2(Exponent);
    
        constexpr operator double() const
        {
            return value;
        }
    };
    
    template <bool Negative>
    struct zero{
        static constexpr double value = (Negative ? -0.0 : 0.0);
    
        constexpr operator double() const
        {
            return value;
        }
    };
    
    template <bool Negative>
    struct inf{
        static constexpr double value = Negative
           ? -::std::numeric_limits<double>::infinity()
           : ::std::numeric_limits<double>::infinity();
    
        constexpr operator double() const
        {
            return value;
        }
    };
    
    struct qnan{
        static constexpr double value = ::std::numeric_limits<double>::quiet_NaN();
    
        constexpr operator double() const
        {
            return value;
        }
    };
    
    struct snan{
        static constexpr double value = ::std::numeric_limits<double>::signaling_NaN();
    
        constexpr operator double() const
        {
            return value;
        }
    };
    
    enum class encode_type{
        number, zero, inf, qnan, snan
    };
    
    constexpr ::std::tuple<::std::uint64_t, int> encode_impl(double x)
    {
        x = ::sprout::fabs(x);
    
        if(x < 1.0){
            int i = 0;
    
            while(::std::int64_t(x) != x){
                x *= 2.0;
                i--;
            }
    
            return ::std::make_tuple(::std::uint64_t(x), i);
        }else{
            int i = 0;
            
            while(x > 1.0){
                x /= 2.0;
                i++;
            }
    
            return ::std::make_tuple(::std::uint64_t(x * pow2(53)), i - 53);
        }
    }
    
    constexpr ::std::tuple<encode_type, bool, ::std::uint64_t, int> encode(double x)
    {
        if(x == 0.0)
            return ::std::make_tuple(encode_type::zero, ::sprout::signbit(x), 0, 0);
        else if(x == inf<false>::value)
            return ::std::make_tuple(encode_type::inf, false, 0, 0);
        else if(x == inf<true>::value)
            return ::std::make_tuple(encode_type::inf, true, 0, 0);
        else if(::sprout::isnan(x))
            return ::std::make_tuple(encode_type::qnan, false, 0, 0);
        else
            return ::std::make_tuple(
                encode_type::number, x < 0.0,
                ::std::get<0>(encode_impl(x)),
                ::std::get<1>(encode_impl(x)));
    }
    
    template <encode_type T, bool Negative, ::std::uint64_t Significand, int Exponent>
    struct get_encoded_type{
        using type = double_<Negative, Significand, Exponent>;
    };
    
    template <bool Negative, ::std::uint64_t Significand, int Exponent>
    struct get_encoded_type<encode_type::zero, Negative, Significand, Exponent>{
        using type = zero<Negative>;
    };
    
    template <bool Negative, ::std::uint64_t Significand, int Exponent>
    struct get_encoded_type<encode_type::inf, Negative, Significand, Exponent>{
        using type = inf<Negative>;
    };
    
    template <bool Negative, ::std::uint64_t Significand, int Exponent>
    struct get_encoded_type<encode_type::qnan, Negative, Significand, Exponent>{
        using type = qnan;
    };
    
    template <bool Negative, ::std::uint64_t Significand, int Exponent>
    struct get_encoded_type<encode_type::snan, Negative, Significand, Exponent>{
        using type = snan;
    };
    
    template <encode_type T, bool Negative, ::std::uint64_t Significand, int Exponent>
    using get_encoded_type_t = typename get_encoded_type<T, Negative, Significand, Exponent>::type;
    
    template <typename T>
    struct is_encoded_double : ::std::false_type{
    };

    template <bool Negative, ::std::uint64_t Significand, int Exponent>
    struct is_encoded_double<double_<Negative, Significand, Exponent>>
        : ::std::true_type
    {
    };

    template <bool Negative>
    struct is_encoded_double<zero<Negative>> : ::std::true_type{
    };

    template <bool Negative>
    struct is_encoded_double<inf<Negative>> : ::std::true_type{
    };

    template <>
    struct is_encoded_double<qnan> : ::std::true_type{
    };

    template <>
    struct is_encoded_double<snan> : ::std::true_type{
    };

    template <typename T>
    constexpr bool is_encoded_double_v = is_encoded_double<T>::value;
}

#define D_T(x)\
    ::a::get_encoded_type_t<\
        ::std::get<0>(::a::encode(x)),\
        ::std::get<1>(::a::encode(x)),\
        ::std::get<2>(::a::encode(x)),\
        ::std::get<3>(::a::encode(x))>

#define D_(x) D_T(x){}

namespace a{
    template <
        typename T, typename U,
        ::std::enable_if_t<is_encoded_double_v<T> && is_encoded_double_v<U>>* = nullptr
    >
    constexpr auto operator+(T, U)
    {
        return D_(T::value + U::value);
    }
}

int main()
{
    constexpr auto x = D_(0.1234567890123456);
    constexpr auto y = D_(0.01234567890123456);
    
    std::cout.setf(std::ios::scientific);
    std::cout.precision(16);
    
    std::cout << 0.1234567890123456 + 0.01234567890123456 << std::endl;
    std::cout << decltype(x + y)::value << std::endl;
    std::cout << x + y << std::endl << std::endl;
    
    std::cout << std::numeric_limits<double>::denorm_min() << std::endl;
    std::cout << D_T(std::numeric_limits<double>::denorm_min())::value << std::endl << std::endl;
    
    std::cout << 1e100 << std::endl;
    std::cout << D_T(1e100)::value << std::endl << std::endl;
    
    std::cout << std::numeric_limits<double>::max() << std::endl;
    std::cout << D_T(std::numeric_limits<double>::max())::value << std::endl;
    
    return 0;
}

実行結果

1.3580246791358017e-01
1.3580246791358017e-01
1.3580246791358017e-01

4.9406564584124654e-324
4.9406564584124654e-324

1.0000000000000000e+100
1.0000000000000000e+100

1.7976931348623157e+308
1.7976931348623157e+308