C++ big unsigned integer class

Question

Here is a big unsigned integer type I intend to use for learning and "fun" with cryptography.

It uses a std::vector for storing blocks of unsigned integers, the exact type of which can be changed with the template argument block_t. A double_block_t is used internally for the various math operations. This is selected automatically via a traits class (passing it as an individual template argument seemed unnecessary and error prone).

The last block in the vector must never be zero. A zero value is represented by an empty vector.

Some utility functions are omitted; the files below contain the most interesting parts of the code. Here are links to the full repository and a single-file online version.

rsa_math__big_uint.h: The class definition - constructors, utility functions, and operators.

#pragma once

#include "rsa__debug.h"
#include "rsa_math__utils.h"
#include "rsa_math_ops__operations.h"

#include <algorithm>
#include <cstdint>
#include <stdexcept>
#include <type_traits>
#include <vector>

namespace rsa
{

    namespace math
    {

        template<class block_t>
        struct block_traits;

        template<>
        struct block_traits<std::uint8_t> { using double_block_type = std::uint16_t; };

        template<>
        struct block_traits<std::uint16_t> { using double_block_type = std::uint32_t; };

        template<>
        struct block_traits<std::uint32_t> { using double_block_type = std::uint64_t; };


        template<class block_t>
        class big_uint
        {
        public:

            using block_type = block_t;
            using double_block_type = typename block_traits<block_type>::double_block_type;
            using data_type = std::vector<block_type>;
            using block_index_type = std::size_t;
            using bit_index_type = std::size_t;

            static_assert(utils::is_uint_v<block_type>, "`block_type` must be an unsigned integer.");
            static_assert(utils::is_uint_v<double_block_type>, "`double_block_type` must be an unsigned integer.");
            static_assert(utils::digits<double_block_type>() == 2 * utils::digits<block_type>(), "`double_block_type` have twice as many digits as `block_type`.");

#pragma region constructors

            big_uint();

            template<class uint_t>
            explicit big_uint(uint_t n);

            big_uint(std::initializer_list<block_type> block_values);

            explicit big_uint(block_index_type block_count, block_type block_value);

            template<class inputit_t>
            explicit big_uint(inputit_t first, inputit_t last);

            big_uint(big_uint const&) = default;
            big_uint(big_uint&&) = default;

#pragma endregion

#pragma region assignment

            big_uint& operator=(big_uint const&) = default;
            big_uint& operator=(big_uint&&) = default;

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator=(uint_t n);

#pragma endregion

#pragma region general

            template<class uint_t>
            uint_t to_uint() const;

            bool is_zero() const;

            bool get_bit(bit_index_type i) const;
            void set_bit(bit_index_type i, bool value);
            void flip_bit(bit_index_type i);

            bit_index_type get_most_significant_bit() const;

            data_type& data();
            data_type const& data() const;

#pragma endregion

#pragma region bitwise operators

            big_uint& operator&=(big_uint const& b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator&=(uint_t n);

            big_uint& operator|=(big_uint const& b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator|=(uint_t n);

            big_uint& operator^=(big_uint const& b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator^=(uint_t n);

            big_uint& operator<<=(bit_index_type n);
            big_uint& operator>>=(bit_index_type n);

#pragma endregion

#pragma region math operators

            big_uint& operator+=(big_uint const& b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator+=(uint_t n);

            big_uint& operator-=(big_uint const& b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator-=(uint_t n);

            big_uint& operator*=(big_uint const& b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator*=(uint_t n);

            big_uint& operator/=(big_uint b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator/=(uint_t n);

            big_uint& operator%=(big_uint b);

            template<class uint_t, typename = utils::enable_if_uint_t<uint_t>>
            big_uint& operator%=(uint_t n);

            big_uint& operator++();
            big_uint operator++(int);

            big_uint& operator--();
            big_uint operator--(int);

#pragma endregion

        private:

            data_type m_data;
        };

        using big_uint_8 = big_uint<std::uint8_t>;
        using big_uint_16 = big_uint<std::uint16_t>;
        using big_uint_32 = big_uint<std::uint32_t>;

#pragma region members - construct

        template<class block_t>
        big_uint<block_t>::big_uint()
        {

        }

        template<class block_t>
        template<class uint_t>
        big_uint<block_t>::big_uint(uint_t n):
            big_uint()
        {
            static_assert(utils::is_uint_v<uint_t>, "`uint_t` must be an unsigned integer.");

            // shifting by >= the number digits in the type is undefined behaviour.
            constexpr bool can_rshift = (utils::digits<block_type>() < utils::digits<uint_t>());

            while (n != uint_t{ 0 })
            {
                // integer promotion, conversion to greater rank, implicit conversion to block_type
                m_data.push_back(utils::max<block_type>() & n);

                if (can_rshift)
                    n >>= utils::digits<block_type>();
                else
                    n = uint_t{ 0 };
            }
        }

        template<class block_t>
        big_uint<block_t>::big_uint(std::initializer_list<block_type> block_values):
            m_data(block_values)
        {
            utils::trim(*this);
        }

        template<class block_t>
        big_uint<block_t>::big_uint(block_index_type block_count, block_type block_value):
            m_data(block_count, block_value)
        {
            utils::trim(*this);
        }

        template<class block_t>
        template<class inputit_t>
        big_uint<block_t>::big_uint(inputit_t first, inputit_t last):
            m_data(first, last)
        {
            utils::trim(*this);
        }

#pragma endregion

#pragma region members - general

        template<class block_t>
        template<class uint_t>
        uint_t big_uint<block_t>::to_uint() const
        {
            static_assert(utils::is_uint_v<uint_t>, "`uint_t` must be an unsigned integer.");

            // it's much easier to static_assert / throw here if uint_t may be too small.
            // checking the actual value would be a lot more work.
            {
                static_assert(utils::digits<block_t>() <= utils::digits<uint_t>(), "uint_t may be too small to represent this number.");

                if (m_data.size() * utils::digits<block_type>() > utils::digits<uint_t>())
                    throw std::range_error("uint_t may be too small to represent this number.");
            }

            auto result = uint_t{ 0 };

            if (m_data.empty())
                return result;

            for (auto i = std::size_t{ 0 }; i != m_data.size(); ++i)
                result |= (uint_t{ m_data[i] } << (i * utils::digits<block_type>()));

            return result;
        }

        template<class block_t>
        bool big_uint<block_t>::is_zero() const
        {
            return m_data.empty();
        }

        template<class block_t>
        bool big_uint<block_t>::get_bit(bit_index_type i) const
        {
            auto block_index = i / utils::digits<block_type>();

            if (m_data.size() <= block_index)
                return false;

            auto block_bit = i - (block_index * utils::digits<block_type>());
            debug::die_if(block_bit >= utils::digits<block_type>());

            return bool((m_data[block_index] >> block_bit) & 1u);
        }

        template<class block_t>
        void big_uint<block_t>::set_bit(bit_index_type i, bool value)
        {
            auto block_index = i / utils::digits<block_type>();

            if (m_data.size() <= block_index)
                m_data.resize(block_index + 1u);

            auto block_bit = i - (block_index * utils::digits<block_type>());
            debug::die_if(block_bit >= utils::digits<block_type>());

            auto mask = block_type(block_type{ 1u } << block_bit);
            m_data[block_index] |= mask & block_type(block_type{ value } << block_bit);
        }

        template<class block_t>
        void big_uint<block_t>::flip_bit(bit_index_type i)
        {
            auto block_index = i / utils::digits<block_type>();

            if (m_data.size() <= block_index)
                m_data.resize(block_index + 1u);

            auto block_bit = i - (block_index * utils::digits<block_type>());
            debug::die_if(block_bit >= utils::digits<block_type>());

            auto mask = block_type(block_type{ 1u } << block_bit);
            m_data[block_index] ^= mask;
        }

        template<class block_t>
        typename big_uint<block_t>::bit_index_type big_uint<block_t>::get_most_significant_bit() const
        {
            if (is_zero())
                throw std::logic_error("number must not be zero.");

            auto block = m_data.back();
            auto count = std::uint8_t{ 0u };

            while (block != block_type{ 1u })
            {
                ++count;
                block >>= 1u;
            }

            return bit_index_type{ count + (m_data.size() - 1u) * utils::digits<block_type>() };
        }

        template<class block_t>
        typename big_uint<block_t>::data_type& big_uint<block_t>::data()
        {
            return m_data;
        }

        template<class block_t>
        typename big_uint<block_t>::data_type const& big_uint<block_t>::data() const
        {
            return m_data;
        }

#pragma endregion

#pragma region members - assignment

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator=(uint_t n)
        {
            return (*this = big_uint(n));
        }

#pragma endregion

#pragma region members - bitwise operators

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator&=(big_uint const& b)
        {
            ops::bit_and_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator&=(uint_t n)
        {
            return operator&=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator|=(big_uint const& b)
        {
            ops::bit_or_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator|=(uint_t n)
        {
            return operator|=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator^=(big_uint const& b)
        {
            ops::bit_xor_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator^=(uint_t n)
        {
            return operator^=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator<<=(bit_index_type n)
        {
            ops::lshift_assign(*this, n);
            return *this;
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator>>=(bit_index_type n)
        {
            ops::rshift_assign(*this, n);
            return *this;
        }

#pragma endregion

#pragma region members - math operators

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator+=(big_uint const& b)
        {
            ops::add_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator+=(uint_t n)
        {
            return operator+=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator-=(big_uint const& b)
        {
            ops::sub_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator-=(uint_t n)
        {
            return operator-=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator*=(big_uint const& b)
        {
            ops::mul_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator*=(uint_t n)
        {
            return operator*=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator/=(big_uint b)
        {
            ops::div_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator/=(uint_t n)
        {
            return operator/=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator%=(big_uint b)
        {
            ops::mod_assign(*this, b);
            return *this;
        }

        template<class block_t>
        template<class uint_t, typename>
        big_uint<block_t>& big_uint<block_t>::operator%=(uint_t n)
        {
            return operator%=(big_uint(n));
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator++()
        {
            return operator+=(1u);
        }

        template<class block_t>
        big_uint<block_t> big_uint<block_t>::operator++(int)
        {
            auto temp = *this;
            operator++();
            return temp;
        }

        template<class block_t>
        big_uint<block_t>& big_uint<block_t>::operator--()
        {
            return operator-=(1u);
        }

        template<class block_t>
        big_uint<block_t> big_uint<block_t>::operator--(int)
        {
            auto temp = *this;
            operator--();
            return temp;
        }

#pragma endregion

#pragma region comparison operators

        template<class block_t>
        bool operator==(big_uint<block_t> const& a, big_uint<block_t> const& b)
        {
            return std::equal(a.data().begin(), a.data().end(), b.data().begin(), b.data().end());
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator==(big_uint<block_t> const& a, uint_t b)
        {
            return (a == big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator==(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) == b);
        }

        template<class block_t>
        bool operator!=(big_uint<block_t> const& a, big_uint<block_t> const& b)
        {
            return !(a == b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator!=(big_uint<block_t> const& a, uint_t b)
        {
            return (a != big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator!=(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) != b);
        }

        template<class block_t>
        bool operator<(big_uint<block_t> const& a, big_uint<block_t> const& b)
        {
            if (a.data().size() < b.data().size())
                return true;

            if (b.data().size() < a.data().size())
                return false;

            return std::lexicographical_compare(a.data().rbegin(), a.data().rend(), b.data().rbegin(), b.data().rend());
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator<(big_uint<block_t> const& a, uint_t b)
        {
            return (a < big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator<(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) < b);
        }

        template<class block_t>
        bool operator>(big_uint<block_t> const& a, big_uint<block_t> const& b)
        {
            return (b < a);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator>(big_uint<block_t> const& a, uint_t b)
        {
            return (a > big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator>(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) > b);
        }

        template<class block_t>
        bool operator<=(big_uint<block_t> const& a, big_uint<block_t> const& b)
        {
            return !(b < a);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator<=(big_uint<block_t> const& a, uint_t b)
        {
            return (a <= big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator<=(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) <= b);
        }

        template<class block_t>
        bool operator>=(big_uint<block_t> const& a, big_uint<block_t> const& b)
        {
            return !(a < b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator>=(big_uint<block_t> const& a, uint_t b)
        {
            return (a >= big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        bool operator>=(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) >= b);
        }

#pragma endregion

#pragma region bitwise operators

        template<class block_t>
        big_uint<block_t> operator&(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a &= b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator&(big_uint<block_t> a, uint_t b)
        {
            return (a &= big_uint<block_t>(b));
        }

        template<class block_t>
        big_uint<block_t> operator|(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a |= b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator|(big_uint<block_t> a, uint_t b)
        {
            return (a |= big_uint<block_t>(b));
        }

        template<class block_t>
        big_uint<block_t> operator^(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a ^= b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator^(big_uint<block_t> a, uint_t b)
        {
            return (a ^= big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator<<(big_uint<block_t> a, uint_t b)
        {
            return (a <<= b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator>>(big_uint<block_t> a, uint_t b)
        {
            return (a >>= b);
        }

#pragma endregion

#pragma region math operators

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator+(big_uint<block_t> a, uint_t b)
        {
            return (a += big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator+(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) += b);
        }

        template<class block_t>
        big_uint<block_t> operator+(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a += b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator-(big_uint<block_t> a, uint_t b)
        {
            return (a -= big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator-(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) -= b);
        }

        template<class block_t>
        big_uint<block_t> operator-(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a -= b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator*(big_uint<block_t> a, uint_t b)
        {
            return (a *= big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator*(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) *= b);
        }

        template<class block_t>
        big_uint<block_t> operator*(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a *= b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator/(big_uint<block_t> a, uint_t b)
        {
            return (a /= big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator/(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) /= b);
        }

        template<class block_t>
        big_uint<block_t> operator/(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a /= b);
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator%(big_uint<block_t> a, uint_t b)
        {
            return (a %= big_uint<block_t>(b));
        }

        template<class block_t, class uint_t, typename = utils::enable_if_uint_t<uint_t>>
        big_uint<block_t> operator%(uint_t a, big_uint<block_t> const& b)
        {
            return (big_uint<block_t>(a) %= b);
        }

        template<class block_t>
        big_uint<block_t> operator%(big_uint<block_t> a, big_uint<block_t> const& b)
        {
            return (a %= b);
        }

#pragma endregion

    } // math

} // rsa

rsa_math_ops__operations.h - The core of the various bitwise and math operations. The bitwise operations simply apply the operation to each block (with some fiddling for different sizes of vector). The math operations are implemented with a standard "apply and carry" approach. The division / modulus is based on the Hacker's Delight implementations of Knuth's Algorithm D.

#pragma once

#include "rsa__debug.h"
#include "rsa_math__utils.h"

#include <algorithm>
#include <stdexcept>

namespace rsa
{

    namespace math
    {

        template<class block_t>
        class big_uint;

        namespace ops
        {

            template<class block_t>
            void bit_and_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
            {
                const auto min_size = std::min(a.data().size(), b.data().size());

                for (auto i = std::size_t{ 0 }; i != min_size; ++i)
                    a.data()[i] &= b.data()[i];

                a.data().resize(min_size);
            }

            template<class block_t>
            void bit_or_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
            {
                const auto min_size = std::min(a.data().size(), b.data().size());

                for (auto i = std::size_t{ 0 }; i != min_size; ++i)
                    a.data()[i] |= b.data()[i];

                std::copy(b.data().begin() + min_size, b.data().end(), std::back_inserter(a.data()));
            }

            template<class block_t>
            void bit_xor_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
            {
                const auto min_size = std::min(a.data().size(), b.data().size());

                for (auto i = std::size_t{ 0 }; i != min_size; ++i)
                    a.data()[i] ^= b.data()[i];

                std::copy(b.data().begin() + min_size, b.data().end(), std::back_inserter(a.data()));

                utils::trim(a);
            }

            template<class block_t>
            void lshift_assign(big_uint<block_t>& a, typename big_uint<block_t>::bit_index_type n)
            {
                using bit_index_t = typename big_uint<block_t>::bit_index_type;
                constexpr auto block_digits = utils::digits<block_t>();

                if (n == bit_index_t{ 0 })
                    return;

                if (a.is_zero())
                    return;

                // shift by whole blocks
                if (n >= block_digits)
                {
                    auto blocks = n / block_digits;
                    a.data().insert(a.data().begin(), blocks, block_t{ 0 });

                    n -= (blocks * block_digits);

                    if (n == bit_index_t{ 0 })
                        return;
                }

                debug::die_if(n >= block_digits);
                const auto carry_shift = block_t(block_digits - n);
                auto carry = block_t{ 0 };

                // shift by partial blocks
                for (auto& block : a.data())
                {
                    // take high bits, shift them to low bits for next block (cast to fix type from integer promotion)
                    const auto carry_out = block_t(block >> carry_shift);

                    // shift low bits to high, apply carry bits
                    block = (block << n) | carry;

                    carry = carry_out;
                }

                if (carry != block_t{ 0 })
                    a.data().push_back(carry);

                debug::die_if(utils::has_extra_empty_blocks(a));
            }

            template<class block_t>
            void rshift_assign(big_uint<block_t>& a, typename big_uint<block_t>::bit_index_type n)
            {
                using bit_index_t = typename big_uint<block_t>::bit_index_type;
                constexpr auto block_digits = utils::digits<block_t>();

                if (n == bit_index_t{ 0 })
                    return;

                if (a.is_zero())
                    return;

                // shift by whole blocks
                if (n >= block_digits)
                {
                    auto blocks = n / block_digits;
                    a.data().erase(a.data().begin(), a.data().begin() + std::min<std::size_t>(blocks, a.data().size()));

                    if (a.is_zero())
                        return;

                    n -= (blocks * block_digits);

                    if (n == bit_index_t{ 0 })
                        return;
                }

                debug::die_if(n >= block_digits);
                const auto carry_shift = block_t(block_digits - n);
                auto carry = block_t{ 0 };

                // shift by partial blocks
                for (auto i_block = a.data().rbegin(); i_block != a.data().rend(); ++i_block)
                {
                    auto& block = *i_block;

                    // take low bits, shift them to high bits for the next block (cast to fix type from integer promotion)
                    const auto carry_out = block_t(block << carry_shift);

                    // shift high bits to low, apply carry bits
                    block = (block >> n) | carry;

                    carry = carry_out;
                }

                utils::trim(a);
            }

            template<class block_t>
            void add_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
            {
                using double_block_t = typename big_uint<block_t>::double_block_type;

                if (b.is_zero())
                    return;

                if (a.is_zero())
                {
                    a = b;
                    return;
                }

                auto& a_data = a.data();
                auto const& b_data = b.data();
                const auto min_size = std::min(a_data.size(), b_data.size());

                auto carry = double_block_t{ 0 };

                // both a and b have data
                for (auto i = std::size_t{ 0 }; i != min_size; ++i)
                {
                    carry += static_cast<double_block_t>(a_data[i]) + static_cast<double_block_t>(b_data[i]);
                    a_data[i] = static_cast<block_t>(carry);
                    carry >>= utils::digits<block_t>();
                }

                // ran out of data in a, copy over the rest of b
                a_data.insert(a_data.end(), b_data.begin() + min_size, b_data.end());

                // add carry
                for (auto i = min_size; i != a_data.size() && (carry != double_block_t{ 0 }); ++i)
                {
                    carry += static_cast<double_block_t>(a_data[i]);
                    a_data[i] = static_cast<block_t>(carry);
                    carry >>= utils::digits<block_t>();
                }

                // extend a if necessary
                if (carry)
                    a_data.push_back(static_cast<block_t>(carry));
            }

            template<class block_t>
            void sub_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
            {
                using double_block_t = typename big_uint<block_t>::double_block_type;

                if (b.is_zero())
                    return;

                if (b > a)
                    throw std::invalid_argument("cannot subtract larger value from smaller one.");

                debug::die_if(a.data().size() < b.data().size());

                auto& a_data = a.data();
                auto const& b_data = b.data();

                auto borrow = double_block_t{ 0 };

                // both a and b have data
                for (auto i = std::size_t{ 0 }; i != b_data.size(); ++i)
                {
                    borrow = static_cast<double_block_t>(a_data[i]) - static_cast<double_block_t>(b_data[i]) - borrow;
                    a_data[i] = static_cast<block_t>(borrow);
                    borrow = (borrow >> utils::digits<block_t>()) & double_block_t { 1 };
                }

                // ran out of data in b, subtract borrow
                for (auto i = b_data.size(); i != a_data.size() && (borrow != double_block_t{ 0 }); ++i)
                {
                    borrow = static_cast<double_block_t>(a_data[i]) - borrow;
                    a_data[i] = static_cast<block_t>(borrow);
                    borrow = (borrow >> utils::digits<block_t>()) & double_block_t { 1 };
                }

                utils::trim(a);
            }

            template<class block_t>
            void mul_assign(big_uint<block_t>& lhs, big_uint<block_t> const& rhs)
            {
                using double_block_t = typename big_uint<block_t>::double_block_type;

                if (lhs.is_zero()) return;
                if (rhs.is_zero()) { lhs.data().clear(); return; }

                if (rhs == 1u) return;
                if (lhs == 1u) { lhs = rhs; return; }

                // note: long multiplication relies on:
                // double_block_t holding (max<block_t>() * max<block_t>() + 2 * max<block_t>())
                // which is exactly the case if digits<double_block_t>() == 2 * digits<block_t>();

                {
                    auto b = rhs; // TODO: we only need to copy this if lhs and rhs refer to the same object
                    auto a = std::move(lhs);
                    auto& c = lhs;

                    auto const& a_data = a.data();
                    auto const& b_data = b.data();
                    auto& c_data = c.data();

                    c_data.resize(a_data.size() + b_data.size());

                    for (auto i = std::size_t{ 0 }; i != a_data.size(); ++i)
                    {
                        auto carry = double_block_t{ 0 };

                        for (auto j = std::size_t{ 0 }; j != b_data.size(); ++j)
                        {
                            carry += static_cast<double_block_t>(a_data[i]) * static_cast<double_block_t>(b_data[j]);
                            carry += c_data[i + j];
                            c_data[i + j] = static_cast<block_t>(carry);
                            carry >>= utils::digits<block_t>();
                        }

                        // c_data[i + b_data.size()] is always zero
                        if (carry)
                            c_data[i + b_data.size()] = static_cast<block_t>(carry);
                    }

                    utils::trim(c);
                }
            }

            template<class block_t>
            void divmod(big_uint<block_t>& quotient, big_uint<block_t>& remainder, big_uint<block_t> dividend, big_uint<block_t> divisor)
            {
                using double_block_t = typename big_uint<block_t>::double_block_type;

                quotient.data().clear();
                remainder.data().clear();

                debug::die_if(divisor.is_zero());
                debug::die_if(dividend < divisor);
                debug::die_if(dividend == divisor);
                debug::die_if(dividend.data().size() == 1u && divisor.data().size() == 1u);

                auto const get_num_leading_zeros = [] (block_t x)
                {
                    debug::die_if(x == 0);

                    auto count = std::size_t{ 0 };

                    while (x != 0)
                    {
                        ++count;
                        x >>= 1;
                    }

                    return utils::digits<block_t>() - count;
                };

                auto const promote = [] (double_block_t b) { return b << utils::digits<block_t>(); };
                auto const demote = [] (double_block_t b) { return b >> utils::digits<block_t>(); };
                auto const checked_sub = [] (block_t& out, block_t a, block_t b) { return ((out = a - b) > a); };

                {
                    auto& d = divisor;
                    auto& n = remainder;
                    remainder = std::move(dividend);
                    auto& q = quotient;

                    q.data().resize(n.data().size() - d.data().size() + 1);

                    // single digit divisor
                    if (d.data().size() == 1)
                    {
                        auto k = double_block_t{ 0 };
                        auto const v = d.data()[0];

                        for (auto i = n.data().size(); i != 0; --i)
                        {
                            auto const index = i - 1;
                            k = (k << utils::digits<block_t>()) + n.data()[index];
                            q.data()[index] = static_cast<block_t>(k / v);
                            k -= static_cast<double_block_t>(q.data()[index]) * v;
                        }

                        n.data().clear();

                        if (k != 0)
                            n.data().push_back(static_cast<block_t>(k));

                        utils::trim(q);

                        return;
                    }

                    auto const b = double_block_t{ 1 } << utils::digits<block_t>();
                    auto const ns = n.data().size(); // m
                    auto const ds = d.data().size(); // n

                    auto shift = get_num_leading_zeros(d.data().back());
                    d <<= shift;
                    n <<= shift;

                    if (n.data().size() == ns)
                        n.data().push_back(block_t{ 0 });

                    for (auto i_outer = (ns - ds) + 1; i_outer != 0; --i_outer)
                    {
                        auto const j = i_outer - 1;

                        // take the top two blocks of n, divide by top block of d, calc remainder
                        auto v = d.data()[ds - 1];
                        auto n_block = static_cast<double_block_t>(promote(n.data()[j + ds]) | n.data()[j + ds - 1]);
                        auto qhat = static_cast<double_block_t>(n_block / v);
                        auto rhat = static_cast<double_block_t>(n_block - qhat * v);

                        // q is too big or (looking at next block) remainder is smaller than what will be taken away
                        while (qhat >= b || (qhat * d.data()[ds - 2]) > (promote(rhat) + n.data()[j + ds - 2]))
                        {
                            qhat -= 1; rhat += v;
                            if (rhat >= b) break;
                        }

                        // qhat is now correct, or 1 too high (extremely rare)

                        // multiply divisor by qhat and subtract from n
                        auto underflow = false;
                        {
                            auto k = double_block_t{ 0 };
                            for (auto i = std::size_t{ 0 }; i != ds; ++i)
                            {
                                auto const p = static_cast<double_block_t>(qhat * d.data()[i]);
                                auto const t = static_cast<double_block_t>(n.data()[i + j] - k - static_cast<block_t>(p));
                                n.data()[i + j] = static_cast<block_t>(t);
                                k = static_cast<double_block_t>(demote(p) - (static_cast<std::make_signed_t<double_block_t>>(t) >> utils::digits<block_t>()));
                            }

                            if (k != 0)
                                underflow |= checked_sub(n.data()[j + ds], n.data()[j + ds], static_cast<block_t>(k));
                        }

                        // set quotient
                        q.data()[j] = static_cast<block_t>(qhat);

                        // underflow! (qhat was 1 too high)
                        // decrement q and add back one divisor to the remainder
                        if (underflow)
                        {
                            q.data()[j] = q.data()[j] - 1;

                            auto k = double_block_t{ 0 };
                            for (auto i = std::size_t{ 0 }; i != ds; ++i)
                            {
                                auto const t = double_block_t{ n.data()[i + j] } + d.data()[i] + k;
                                n.data()[i + j] = static_cast<block_t>(t);
                                k = static_cast<double_block_t>(t >> utils::digits<block_t>());
                            }

                            n.data()[j + ds] += static_cast<block_t>(k);
                        }
                    }

                    utils::trim(q);

                    // shift remainder back
                    utils::trim(n);
                    n >>= shift;
                }
            }

            template<class block_t>
            void div_assign(big_uint<block_t>& lhs, big_uint<block_t> const& rhs)
            {
                using double_block_t = typename big_uint<block_t>::double_block_type;

                if (rhs.is_zero())
                    throw std::invalid_argument("divisor cannot be zero.");

                if (lhs < rhs) { lhs.data().clear(); return; }
                if (lhs == rhs) { lhs.data().clear(); lhs.data().push_back(1); return; }

                if (lhs.data().size() == 1u && rhs.data().size() == 1u)
                {
                    lhs = static_cast<block_t>(lhs.data()[0] / rhs.data()[0]);
                    return;
                }

                {
                    auto q = big_uint<block_t>();
                    auto r = big_uint<block_t>();
                    divmod(q, r, lhs, rhs);
                    lhs = std::move(q);
                }
            }

            template<class block_t>
            void mod_assign(big_uint<block_t>& lhs, big_uint<block_t> const& rhs)
            {
                using double_block_t = typename big_uint<block_t>::double_block_type;

                if (rhs.is_zero())
                    throw std::invalid_argument("divisor cannot be zero.");

                if (lhs < rhs) { return; }
                if (lhs == rhs) { lhs.data().clear(); return; }

                if (lhs.data().size() == 1u && rhs.data().size() == 1u)
                {
                    lhs = static_cast<block_t>(lhs.data()[0] % rhs.data()[0]);
                    utils::trim(lhs);
                    return;
                }

                {
                    auto q = big_uint<block_t>();
                    auto r = big_uint<block_t>();
                    divmod(q, r, lhs, rhs);
                    lhs = std::move(r);
                }
            }

            // utility for testing (when we need both quotient and remainder)
            template<class block_t>
            void divmod_test(big_uint<block_t>& quotient, big_uint<block_t>& remainder, big_uint<block_t> const& dividend, big_uint<block_t> const& divisor)
            {
                quotient.data().clear();
                remainder.data().clear();

                if (divisor.is_zero())
                    throw std::invalid_argument("divisor cannot be zero.");

                if (dividend < divisor) { remainder = dividend; return; }
                if (dividend == divisor) { quotient.data().push_back(1); return; }

                if (dividend.data().size() == 1u && divisor.data().size() == 1u)
                {
                    quotient = static_cast<block_t>(dividend.data()[0] / divisor.data()[0]);
                    remainder = static_cast<block_t>(dividend.data()[0] % divisor.data()[0]);

                    utils::trim(remainder);
                    return;
                }

                divmod(quotient, remainder, dividend, divisor);
            }

        } // ops

    } // math

} // rsa

rsa_math__utils.h

#pragma once

#include <algorithm>
#include <cstdint>
#include <limits>
#include <type_traits>

namespace rsa
{

    namespace math
    {

        template<class block_t>
        class big_uint;

        namespace utils
        {

            template<class uint_t>
            constexpr bool is_uint_v = (std::is_integral_v<uint_t> && std::is_unsigned_v<uint_t> && !std::is_same_v<uint_t, bool>);

            template<class uint_t>
            using enable_if_uint_t = std::enable_if_t<is_uint_v<uint_t>>;


            template<class t>
            constexpr std::uint32_t digits()
            {
                return std::uint32_t(std::numeric_limits<t>::digits);
            }

            template<class t>
            constexpr t max()
            {
                return std::numeric_limits<t>::max();
            }


            template<class block_t>
            bool has_extra_empty_blocks(big_uint<block_t> const& a)
            {
                return
                    (std::find_if(a.data().rbegin(), a.data().rend(),
                        [] (block_t b) { return b != block_t{ 0 }; }).base() !=
                        a.data().end());
            }

            template<class block_t>
            void trim(big_uint<block_t>& a)
            {
                a.data().erase(
                    std::find_if(a.data().rbegin(), a.data().rend(),
                        [] (block_t b) { return b != block_t{ 0 }; }).base(),
                    a.data().end());
            }

        } // utils

    } // math

} // rsa

Any advice / feedback is welcome.

Any tips for neatening the math::ops functions are especially appreciated. The divmod function is quite fiddly and contains a lot of casts, several of which are only necessary for smaller block sizes (i.e. std::uint8_ts and integer promotion don't mix well).

(Templating the class on block size seemed like a good idea at the time, but now it just looks overcomplicated, so I'll probably just fix the block size to std::uint32_t in the class with a typedef.)

Is there anything missing / dodgy looking / surprising?

Suggested usage is just like a normal numeric type. (Note, however, that the various math operators require unsigned types, and don't work with signed ones).

#include <limits>
#include <iostream>

int main()
{
    auto a = rsa::math::big_uint_32(std::numeric_limits<std::uint64_t>::max());
    auto b = (a * a) + 2u * a;

    std::cout << std::hex;
    for (auto block : b.data())
        std::cout << block << " ";
    std::cout << std::endl;
}

Full repository.

Online compiler.

Toby Speight · Accepted Answer · 2018-11-20 18:24:02Z

Nice code; I didn't find any serious issues, so my comments are mostly limited to mere nitpicking.

namespace `math::util`

In has_extra_empty_blocks(), there's no need to continue searching if we don't have a match to begin with:

        template<class block_t>
        bool has_extra_empty_blocks(big_uint<block_t> const& a)
        {
            return
                (std::find_if(a.data().rbegin(), a.data().rend(),
                    [] (block_t b) { return b != block_t{ 0 }; }).base() !=
                    a.data().end());
        }

(Also, it's simpler to just compare the reverse iterator against a.data().rbegin() than to convert to forward iterator). We just want to check that the last element (if there is one) is zero:

        template<class block_t>
        bool has_extra_empty_blocks(big_uint<block_t> const& a)
        {
            return !a.empty() && a.back() == block_t{0};
        }

That looks much more readable, as well as being more efficient.

I think that has_extra_empty_blocks() and trim() would probably make more sense as member functions of big_uint.

namespace `math::ops`

These functions all look like they should be members of big_uint (generally, inlined into the respective operators); that should obviate the need to expose data() publicly.

Loops like this:

            for (auto i = std::size_t{ 0 }; i != min_size; ++i)
                a.data()[i] &= b.data()[i];

look like they would be more naturally written using std::transform():

        template<class block_t>
        void bit_and_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
        {
            if (a.data().size() > b.data().size())
                a.data().resize(b.data().size());

            std::transform(a.data().begin(), a.data().end(), b.data().begin(),
                           a.data().begin(), std::bit_and<block_t>{});
        }

        template<class block_t>
        void bit_or_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
        {
            if (a.data().size() < b.data().size())
                a.data().resize(b.data().size());

            std::transform(b.data().begin(), b.data().end(), a.data().begin(),
                           a.data().begin(), std::bit_or<block_t>{});
        }

        template<class block_t>
        void bit_xor_assign(big_uint<block_t>& a, big_uint<block_t> const& b)
        {
            if (a.data().size() < b.data().size())
                a.data().resize(b.data().size());

            std::transform(b.data().begin(), b.data().end(), a.data().begin(),
                           a.data().begin(), std::bit_xor<block_t>{});

            utils::trim(a);
        }

The shift operators repeat the test of n == 0 before and after shifting by units of block_digits. This could be simplified:

            if (a.is_zero())
                return;

            if (n >= block_digits)
            {
                 ... shift whole blocks
            }

            if (n == bit_index_t{ 0 })
                return;

Addition and subtraction don't need twice the width of block_t. Carry is only ever one bit, and unsigned overflow is well-defined:

            for (auto i = std::size_t{ 0 }; i != min_size; ++i)
            {
                bool carry_out = a_data[i] >= ~b_data[i] + !carry;
                a_data[i] += b_data[i] + carry;
                carry = carry_out;
            }

It is possible to implement multiplication without double_block_t, but I'm not sure that it's worth the effort. Something to consider, perhaps only when double_block_t isn't defined?

Both div_assign and mod_assign declare double_block_t, but never use it.

class `math::big_uint`

I'm not convinced that double_block_type ought to be public.

I'd prefer constraints to static_assert for the template constructors. For example:

        template<class uint_t, typename = std::enable_if_t<utils::is_uint_v<uint_t>>>
        explicit big_uint(uint_t n);

with

    template<class block_t>
        template<class uint_t, typename>
    big_uint<block_t>::big_uint(uint_t n):
        big_uint()
    {

I'm not convinced that this should be explicit - it's a natural, expected conversion. That said, allowing it to be an implicit conversion wouldn't allow us to reduce the overloads of the binary operators, given that they are all template functions.

Consider adding an overload that accepts block_type:

    template<class block_t>
    big_uint<block_t>::big_uint(block_t n)
    {
        if (n)
            m_data.push_back(n);
    }

With that in place, the template constructor could be constrained further, accepting only wider types, and letting the block_t constructor be used for narrower types - which can be achieved by adding && !std::is_assignable_v<block_type, uint_t> to the std::enable_if condition.

It might be useful to also have a converting constructor from big_uint with different block type.

Missing functionality

For RSA operations, you'll want to add a modular exponentiation function. That's fairly easily done using the existing functions, but may be more efficient with a modular multiply.

Thanks. This is really helpful. Is there a particular reason for preferring enable_if for the constructor? static_assert seemed like it would give a more understandable error for what might be a common mistake (big_uint_32(23) vs big_uint_32(23u)). — user673679, Commented Nov 19, 2018 at 13:16
Btw, It seems that bit_or_assign and bit_xor_assign need to use b as the main input sequence for std::transform, because b can be shorter than a. It also looks like I forgot to trim the result after bit_and_assign. :S — user673679, Commented Nov 19, 2018 at 13:27
I suggest the enable_if because it completely removes the invalid versions from participating in overload resolution (which might make the compiler's job easier). I agree that it's nice to be able to offer a better error message! — Toby Speight, Commented Nov 19, 2018 at 13:27
Yes, you're right about the (x)or_assign() - I copied from the and method without thinking! Now fixed; thanks. — Toby Speight, Commented Nov 19, 2018 at 13:28
Yep, the truncation is correct. But if, e.g., a contains blocks 1, 0, 0 and b is 1, 0 then we truncate and do the & and get a result of 0, 0. I've been trying to ensure that the last block is always non-zero, so these empty blocks need to be erased by calling utils::trim. I've been wondering if this "no trailing empty blocks" policy is actually a good idea though, when it comes to encoding byte strings for encryption I presumably need to preserve null values. — user673679, Commented Nov 20, 2018 at 17:51

Stack Exchange Network

C++ big unsigned integer class

1 Answer 1

namespace `math::util`

namespace `math::ops`

class `math::big_uint`

Missing functionality

Your Answer

Not the answer you're looking for? Browse other questions tagged
c++
integer
or ask your own question.

Hot Network Questions

C++ big unsigned integer class

1 Answer 1

namespace math::util

namespace math::ops

class math::big_uint

Missing functionality

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged c++integer or ask your own question.

Related

Hot Network Questions

namespace `math::util`

namespace `math::ops`

class `math::big_uint`

Not the answer you're looking for? Browse other questions tagged
c++
integer
or ask your own question.