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.