Aeu.h

Go to the documentation of this file.

 
#ifndef AEU_MULTIPRECISION
#define AEU_MULTIPRECISION
 
#include <cstdint>
#include <cctype>
#include <cwctype>
#include <iostream>
#include <cassert>
 
#ifdef __CUDACC__
#define gpu __host__ __device__
    #include <cuda/std/utility>
    #include <cuda/std/array>
#else
#define gpu
    #include <utility>
    #include <array>
    #include <vector>
#endif
 
#if defined AESI_UNSAFE
    #warning Enabled nightly mode for the library. Functions and methods input arguments are not checked for validity. Be really gentle
#endif
 
 
#ifdef AESI_GMP_INTEGRATION
#include <gmpxx.h>
#endif
 
namespace {
    using byte = uint8_t;
    using block = uint32_t;
 
    constexpr inline std::size_t bitsInByte = 8;
    constexpr inline std::size_t blockBitLength = sizeof(block) * bitsInByte;
    constexpr inline uint64_t blockBase = 1ULL << blockBitLength;
    constexpr inline block blockMax = 0xff'ff'ff'ffu;
 
    enum class Comparison { equal = 0, less = 1, greater = 2, equivalent = 3 };
}
 
template <std::size_t bitness = 512> requires (bitness % blockBitLength == 0)
class Aeu final {
    static_assert(bitness > sizeof(uint64_t), "Use built-in types for numbers 64-bit or less.");
    
    static constexpr std::size_t blocksNumber = bitness / blockBitLength;
 
#ifdef __CUDACC__
    template <typename T1, typename T2>
    using pair = cuda::std::pair<T1, T2>;
    using blockLine = cuda::std::array<block, blocksNumber>;
#else
    template<typename T1, typename T2>
    using pair = std::pair<T1, T2>;
    using blockLine = std::array<block, blocksNumber>;
#endif
 
    /* -------------------------- @name Class members. ----------------------- */
    blockLine blocks;
    /* ----------------------------------------------------------------------- */
 
    /* ------------------------ @name Helper functions. ---------------------- */
    gpu static constexpr auto addLine(blockLine& dst, const blockLine& src) noexcept -> uint64_t {
        uint64_t carryOut = 0;
        for (std::size_t i = 0; i < blocksNumber; ++i) {
            uint64_t sum = static_cast<uint64_t>(dst[i]) + static_cast<uint64_t>(src[i]) + carryOut;
            carryOut = sum / blockBase; dst[i] = static_cast<block>(sum % blockBase);
        }
        return carryOut;
    }
 
    gpu static constexpr auto makeComplement(blockLine& line) noexcept -> uint64_t {
        uint64_t carryOut = 1;
        for(std::size_t i = 0; i < blocksNumber; ++i) {
            const uint64_t sum = blockBase - 1ULL - static_cast<uint64_t>(line[i]) + carryOut;
            carryOut = sum / blockBase; line[i] = static_cast<block>(sum % blockBase);
        }
        return carryOut;
    }
    /* ----------------------------------------------------------------------- */
 
public:
    /* --------------------- @name Different constructors. ------------------- */
    constexpr Aeu() noexcept = default;
 
    constexpr ~Aeu() noexcept = default;
 
    constexpr Aeu(const Aeu& copy) noexcept = default;
 
    constexpr Aeu& operator=(const Aeu& other) = default;
 
    template <typename Integral> requires (std::is_unsigned_v<Integral>)
    gpu constexpr Aeu(Integral value) noexcept {
        if(value != 0) {
            unsigned long long tValue = static_cast<unsigned long long>(value);
            for (std::size_t i = 0; i < blocksNumber; ++i) {
                blocks[i] = static_cast<block>(tValue % blockBase);
                tValue /= blockBase;
            }
        } else
            blocks = {};
    }
 
    template <typename Char> requires (std::is_same_v<Char, char> || std::is_same_v<Char, wchar_t>)
    gpu constexpr Aeu(const Char* data, std::size_t size) noexcept : Aeu {} {
        if(size == 0) return;
 
        std::size_t position = 0;
 
        if constexpr (std::is_same_v<Char, char>) {
            for(; position < size && !std::isalnum(data[position]); ++position) ;
        } else {
            for(; position < size && !std::iswalnum(static_cast<wint_t>(data[position])); ++position) ;
        }
 
        if(position == size)
            return;
 
        const auto base = [&data, &size, &position] {
            if (data[position] == Char('0') && size > position + 1) {
                switch (data[position + 1]) {
                    case Char('b'): case Char('B'):
                        position += 2; return 2u;
                    case Char('o'): case Char('O'):
                        position += 2; return 8u;
                    case Char('x'): case Char('X'):
                        position += 2; return 16u;
                    default:
                        return 10u;
                }
            } return 10u;
        } ();
 
        for(; position < size; ++position) {
            const auto digit = [] (Char ch) {
                if(Char('0') <= ch && ch <= Char('9'))
                    return static_cast<unsigned>(ch) - static_cast<unsigned>(Char('0'));
                if(Char('a') <= ch && ch <= Char('f'))
                    return static_cast<unsigned>(ch) - static_cast<unsigned>(Char('a')) + 10u;
                if(Char('A') <= ch && ch <= Char('F'))
                    return static_cast<unsigned>(ch) - static_cast<unsigned>(Char('A')) + 10u;
                return 99u;
            } (data[position]);
 
            if(digit < base) {
                *this *= base;
                *this += digit;
            }
        }
    }
 
    template <typename Char, std::size_t arrayLength> requires (arrayLength > 1 && (std::is_same_v<Char, char> || std::is_same_v<Char, wchar_t>))
    gpu constexpr Aeu(const Char (&literal)[arrayLength]) noexcept : Aeu(literal, arrayLength) {}
 
    template <typename String, typename Char = typename String::value_type> requires (std::is_same_v<std::basic_string<Char>,
        std::decay_t<String>> || std::is_same_v<std::basic_string_view<Char>, std::decay_t<String>>)
    gpu constexpr Aeu(String&& stringView) noexcept : Aeu(stringView.data(), stringView.size()) {}
 
    template <typename String, typename Char = typename String::value_type> requires (std::is_same_v<std::basic_string<Char>,
        std::decay_t<String>> || std::is_same_v<std::basic_string_view<Char>, std::decay_t<String>>)
    gpu constexpr Aeu(const String& stringView) noexcept : Aeu(stringView.data(), stringView.size()) {}
 
#ifdef AESI_GMP_INTEGRATION
    constexpr Aeu(const mpz_class& value) : Aeu {} {
        const auto bitLength = mpz_sizeinbase(value.get_mpz_t(), 2);
        auto tBlocksNumber = 1 + (bitLength - 1) / sizeof(block) * 8;
 
        std::vector<block> buffer (tBlocksNumber, 0u);
        mpz_export(buffer.data(), &tBlocksNumber, -1, sizeof(block), -1, 0, value.get_mpz_t());
 
        for (std::size_t i = 0; i < tBlocksNumber; ++i)
            setBlock(i, buffer[i]);
    }
#endif
    /* ----------------------------------------------------------------------- */
 
 
    /* --------------------- @name Arithmetic operators. --------------------- */
    /* ------------------------- @name Unary operators. -------------------------- */
        [[nodiscard]]
        gpu constexpr auto operator+() const noexcept -> Aeu { return *this; }
 
        [[nodiscard]]
        gpu constexpr auto operator-() const noexcept -> Aeu {
            Aeu result = *this; makeComplement(result.blocks); return result;
        }
 
        gpu constexpr auto operator++() noexcept -> Aeu& {
            for(std::size_t i = 0; i < blocksNumber; ++i) {
                if(blocks[i] < blockMax) {
                    ++blocks[i];
                    break;
                } else blocks[i] = 0u;
            }
            return *this;
        }
 
        gpu constexpr auto operator++(int) & noexcept -> Aeu { Aeu old = *this; operator++(); return old; }
 
        gpu constexpr auto operator--() noexcept -> Aeu& {
            for(std::size_t i = 0; i < blocksNumber; ++i) {
                if(blocks[i] > 0u) {
                    --blocks[i]; break;
                }
                blocks[i] = blockMax;
            }
            return *this;
        }
 
        gpu constexpr auto operator--(int) & noexcept -> Aeu { Aeu old = *this; operator--(); return old; }
    /* --------------------------------------------------------------------------- */
 
 
    /* ------------------------ @name Addition operators. ------------------------ */
        template <typename Unsigned> requires (std::is_unsigned_v<Unsigned>) [[nodiscard]]
        gpu constexpr friend auto operator+(const Aeu& addition, Unsigned addendum) noexcept -> Aeu {
            Aeu result = addition; result += addendum; return result;
        }
 
        [[nodiscard]]
        gpu constexpr friend auto operator+(const Aeu& addition, const Aeu& addendum) noexcept -> Aeu {
            Aeu result = addition; result += addendum; return result;
        }
 
        template <typename Unsigned> requires (std::is_unsigned_v<Unsigned>)
        gpu constexpr friend auto operator+=(Aeu& addition, Unsigned addendum) noexcept -> Aeu& {
            for(std::size_t i = 0; i < blocksNumber; ++i) {
                const auto currentSum = static_cast<uint64_t>(addition.blocks[i]) + static_cast<uint64_t>(addendum);
                addendum = static_cast<Unsigned>(currentSum / blockBase); addition.blocks[i] = static_cast<block>(currentSum % blockBase);
            }
            return addition;
        }
 
        gpu constexpr friend auto operator+=(Aeu& addition, const Aeu& addendum) noexcept -> Aeu& {
            addLine(addition.blocks, addendum.blocks); return addition;
        }
    /* --------------------------------------------------------------------------- */
 
 
    /* ----------------------- @name Subtraction operators. ---------------------- */
        [[nodiscard]]
        gpu constexpr friend auto operator-(const Aeu& subtraction, const Aeu& subtrahend) noexcept -> Aeu {
            Aeu result = subtraction; result -= subtrahend; return result;
        }
 
        gpu constexpr friend auto operator-=(Aeu& subtraction, const Aeu& subtrahend) noexcept -> Aeu& {
            return subtraction += -subtrahend;
        }
    /* --------------------------------------------------------------------------- */
 
 
    /* --------------------- @name Multiplication operators. --------------------- */
        template <typename Unsigned> requires (std::is_unsigned_v<Unsigned>) [[nodiscard]]
        gpu constexpr friend auto operator*(Aeu& multiplication, Unsigned factor) noexcept -> Aeu {
            Aeu result = multiplication; result *= factor; return result;
        }
 
        [[nodiscard]]
        gpu constexpr friend auto operator*(const Aeu& multiplication, const Aeu& factor) noexcept -> Aeu {
            Aeu result = multiplication; result *= factor; return result;
        }
 
        template <typename Unsigned> requires (std::is_unsigned_v<Unsigned>)
        gpu constexpr friend auto operator*=(Aeu& multiplication, Unsigned factor) noexcept -> Aeu& {
            if constexpr (std::is_same_v<Unsigned, uint64_t>) {
                const auto longerLength = multiplication.filledBlocksNumber();
                const auto smallerLength = (factor > blockMax ? 2UL : 1UL);
                blockLine buffer {};
 
                for(std::size_t i = 0; i < longerLength; ++i) {
                    const uint64_t tBlock = multiplication.blocks[i];
                    uint64_t carryOut = 0;
 
                    for(std::size_t j = 0; j < smallerLength && j < buffer.size() - i; ++j) {
                        const auto product = tBlock * (factor >> blockBitLength * j & 0x00'00'00'00'ff'ff'ff'ff) + carryOut;
                        const auto partial = static_cast<uint64_t>(buffer[i + j]) + product % blockBase;
                        carryOut = product / blockBase + partial / blockBase;
                        buffer[i + j] = static_cast<block>(partial % blockBase);
                    }
 
                    if(smallerLength + i < buffer.size())
                        buffer[smallerLength + i] += static_cast<block>(carryOut);
                }
 
                multiplication.blocks = buffer;
                return multiplication;
            } else {
                uint64_t carryOut = 0;
                for (std::size_t i = 0; i < blocksNumber; ++i) {
                    const auto product = static_cast<uint64_t>(factor) * static_cast<uint64_t>(multiplication.blocks[i]) + carryOut;
                    multiplication.blocks[i] = static_cast<block>(product % blockBase); carryOut = product / blockBase;
                }
                return multiplication;
            }
        }
 
        gpu constexpr friend auto operator*=(Aeu& multiplication, const Aeu& factor) noexcept -> Aeu& {
            constexpr auto multiplyLines = [] (const blockLine& longerLine, const std::size_t longerLength,
                                               const blockLine& smallerLine, const std::size_t smallerLength) {
                blockLine buffer {};
 
                for(std::size_t i = 0; i < longerLength; ++i) {
                    const uint64_t tBlock = longerLine[i];
                    uint64_t carryOut = 0;
 
                    for(std::size_t j = 0; j < smallerLength && j < buffer.size() - i; ++j) {
                        const auto product = tBlock * static_cast<uint64_t>(smallerLine[j]) + carryOut;
                        const auto partial = static_cast<uint64_t>(buffer[i + j]) + product % blockBase;
                        carryOut = product / blockBase + partial / blockBase;
                        buffer[i + j] = static_cast<block>(partial % blockBase);
                    }
 
                    if(smallerLength < blocksNumber && smallerLength + i < buffer.size())
                        buffer[smallerLength + i] += static_cast<block>(carryOut);
                }
 
                return buffer;
            };
 
            const std::size_t thisLength = multiplication.filledBlocksNumber();
            if(const std::size_t valueLength = factor.filledBlocksNumber(); thisLength > valueLength)
                multiplication.blocks = multiplyLines(multiplication.blocks, thisLength, factor.blocks, valueLength);
            else
                multiplication.blocks = multiplyLines(factor.blocks, valueLength, multiplication.blocks, thisLength);
 
            return multiplication;
        }
    /* --------------------------------------------------------------------------- */
 
 
    /* ------------------------ @name Division operators. ------------------------ */
        [[nodiscard]]
        gpu constexpr friend auto operator/(const Aeu& division, const Aeu& divisor) noexcept -> Aeu {
            Aeu quotient, _; divide(division, divisor, quotient, _); return quotient;
        }
 
        gpu constexpr friend auto operator/=(Aeu& division, const Aeu& divisor) noexcept -> Aeu& { return division = division / divisor; }
    /* --------------------------------------------------------------------------- */
 
 
    /* ------------------------- @name Modulo operators. ------------------------- */
        [[nodiscard]]
        gpu constexpr friend auto operator%(const Aeu& modulation, const Aeu& modulo) noexcept -> Aeu {
            Aeu _, remainder; divide(modulation, modulo, _, remainder); return remainder;
        }
 
        gpu constexpr friend auto operator%=(Aeu& modulation, const Aeu& modulo) noexcept -> Aeu& {
            return modulation = modulation % modulo;
        }
    /* --------------------------------------------------------------------------- */
    /* ----------------------------------------------------------------------- */
 
 
    /* ----------------------- @name Bitwise operators. ---------------------- */
    [[nodiscard]]
    gpu constexpr auto operator~() const noexcept -> Aeu {
        Aeu result;
        for(std::size_t i = 0; i < blocksNumber; ++i)
            result.blocks[i] = ~blocks[i];
        return result;
    }
 
    [[nodiscard]]
    gpu constexpr friend auto operator^(const Aeu& left, const Aeu& right) noexcept -> Aeu {
        Aeu result = left; result ^= right; return result;
    }
 
    gpu constexpr friend auto operator^=(Aeu& left, const Aeu& right) noexcept -> Aeu& {
        for(std::size_t i = 0; i < blocksNumber; ++i)
            left.blocks[i] ^= right.blocks[i];
        return left;
    }
 
    [[nodiscard]]
    gpu constexpr friend auto operator&(const Aeu& left, const Aeu& right) noexcept -> Aeu {
        Aeu result = left; result &= right; return result;
    }
 
    gpu constexpr friend auto operator&=(Aeu& left, const Aeu& right) noexcept -> Aeu& {
        for(std::size_t i = 0; i < blocksNumber; ++i)
            left.blocks[i] &= right.blocks[i];
        return left;
    }
 
    [[nodiscard]]
    gpu constexpr friend auto operator|(const Aeu& left, const Aeu& right) noexcept -> Aeu {
        Aeu result = left; result |= right; return result;
    }
 
    gpu constexpr friend auto operator|=(Aeu& left, const Aeu& right) noexcept -> Aeu& {
        for(std::size_t i = 0; i < blocksNumber; ++i)
            left.blocks[i] |= right.blocks[i];
        return left;
    }
 
    template <typename Unsigned> requires (std::is_integral_v<Unsigned> && std::is_unsigned_v<Unsigned>) [[nodiscard]]
    gpu constexpr friend auto operator<<(const Aeu& value, Unsigned bitShift) noexcept -> Aeu {
        Aeu result = value; result <<= bitShift; return result;
    }
 
    template <typename Unsigned> requires (std::is_integral_v<Unsigned> && std::is_unsigned_v<Unsigned>)
    gpu constexpr friend auto operator<<=(Aeu& value, Unsigned bitShift) noexcept -> Aeu& {
        if(bitShift >= bitness || bitShift == 0) return value;
 
        const std::size_t quotient = bitShift / blockBitLength, remainder = bitShift % blockBitLength;
        const block stamp = static_cast<block>((1ULL << (blockBitLength - remainder)) - 1);
 
        const std::size_t minI = quotient + (remainder ? 1U : 0U);
        for (std::size_t i = blocksNumber; i-- > minI;)
            value.blocks[i] = (value.blocks[i - quotient] & stamp) << remainder
                | ((value.blocks[i - quotient - (remainder ? 1U : 0U)] & ~stamp) >> (blockBitLength - remainder) % blockBitLength);
 
        value.blocks[quotient] = (value.blocks[0] & stamp) << remainder;
 
        for (std::size_t i = 0; i < quotient; ++i)
            value.blocks[i] = 0;
        return value;
    }
 
    template <typename Unsigned> requires (std::is_integral_v<Unsigned> && std::is_unsigned_v<Unsigned>) [[nodiscard]]
    gpu constexpr friend auto operator>>(const Aeu& value, Unsigned bitShift) noexcept -> Aeu {
        Aeu result = value; result >>= bitShift; return result;
    }
 
    template <typename Unsigned> requires (std::is_integral_v<Unsigned> && std::is_unsigned_v<Unsigned>)
    gpu constexpr friend auto operator>>=(Aeu& value, Unsigned bitShift) noexcept -> Aeu& {
        if(bitShift >= bitness || bitShift == 0) return value;
 
        const std::size_t quotient = bitShift / blockBitLength, remainder = bitShift % blockBitLength;
        const block stamp = static_cast<block>((1ULL << remainder) - 1);
 
        for(std::size_t i = 0; i < blocksNumber - (quotient + (remainder ? 1 : 0)); ++i)
            value.blocks[i] = ((value.blocks[i + quotient + (remainder ? 1 : 0)] & stamp) << (blockBitLength - remainder) % blockBitLength) | (value.blocks[i + quotient] & ~stamp) >> remainder;
 
        value.blocks[blocksNumber - 1 - quotient] = (value.blocks[blocksNumber - 1] & ~stamp) >> remainder;
 
        for(std::size_t i = blocksNumber - quotient; i < blocksNumber; ++i)
            value.blocks[i] = 0;
        return value;
    }
    /* ----------------------------------------------------------------------- */
 
 
    /* --------------------- @name Comparison operators. --------------------- */
    /* ------------------------ @name Equality operators. ------------------------ */
        template <typename Unsigned> requires (std::is_unsigned_v<Unsigned> && sizeof(Unsigned) < 8)
        gpu constexpr friend auto operator==(const Aeu& our, Unsigned other) noexcept -> bool {
            return our.filledBlocksNumber() <= 1 && static_cast<block>(other) == our.blocks[0];
        }
        gpu constexpr friend auto operator==(const Aeu& our, uint64_t other) noexcept -> bool {
            return our.filledBlocksNumber() <= 2 && (static_cast<uint64_t>(our.blocks[1]) << blockBitLength | static_cast<uint64_t>(our.blocks[0])) == other;
        }
 
        template <std::size_t otherBitness>
        gpu constexpr friend auto operator==(const Aeu& our, const Aeu<otherBitness>& other) noexcept -> bool {
            return our.compareTo(other) == Comparison::equal;
        }
 
#ifdef AESI_GMP_INTEGRATION
        template <typename T, typename U>
        constexpr friend auto operator==(const Aeu& our, const __gmp_expr<T, U>& other) noexcept -> bool {
            return our == Aeu(mpz_class(other));
        }
#endif
    /* --------------------------------------------------------------------------- */
 
 
    /* ----------------------- @name Comparison operators. ----------------------- */
        template <typename Unsigned> requires (std::is_unsigned_v<Unsigned> && sizeof(Unsigned) < 8) [[nodiscard]]
        gpu constexpr auto compareTo(Unsigned other) const noexcept -> Comparison {
            using enum Comparison;
 
            if(filledBlocksNumber() > 1) return greater;
            const auto cmp = static_cast<block>(other);
 
            if(blocks[0] > cmp)
                return greater;
 
            if(blocks[0] < cmp)
                return less;
 
            return equal;
        }
 
        [[nodiscard]]
        gpu constexpr auto compareTo(uint64_t other) const noexcept -> Comparison {
            using enum Comparison;
            if(filledBlocksNumber() > 2) return greater;
            const auto base = (static_cast<uint64_t>(blocks[1]) << blockBitLength) | static_cast<uint64_t>(blocks[0]);
            if(base > other)
                return greater; else if(base < other) return less; else return equal;
        }
 
        template <std::size_t otherBitness = bitness> /*requires (otherBitness != bitness)*/ [[nodiscard]]
        gpu constexpr auto compareTo(const Aeu<otherBitness>& other) const noexcept -> Comparison {
            using enum Comparison;
 
            const auto lowerBlockBorder = (blocksNumber < other.totalBlocksNumber() ? blocksNumber : other.totalBlocksNumber());
            for(std::size_t i = lowerBlockBorder; i-- > 0;) {
                const block thisBlock = blocks[i], otherBlock = other.getBlock(i);
                if(thisBlock != otherBlock)
                    return (thisBlock > otherBlock ? greater : less);
            }
 
            if constexpr (otherBitness != blocksNumber * blockBitLength) {
                using enum Comparison;
                if (other.totalBlocksNumber() > blocksNumber) {
                    for (std::size_t i = other.totalBlocksNumber(); i-- > lowerBlockBorder;)
                        if (other.getBlock(i) != 0)
                            return less;
                } else if (blocksNumber > other.totalBlocksNumber()) {
                    for (std::size_t i = blocksNumber; i-- > lowerBlockBorder;)
                        if (blocks[i] != 0)
                            return greater;
                }
            }
 
            return equal;
        }
    /* --------------------------------------------------------------------------- */
 
 
    /* ------------------------ @name Spaceship operators. ----------------------- */
#if (defined(__CUDACC__) || __cplusplus < 202002L || defined (PRE_CPP_20)) && !defined DOXYGEN_SKIP
        gpu constexpr auto operator!=(const Aeu& value) const noexcept -> bool { return !this->operator==(value); }
        gpu constexpr auto operator<(const Aeu& value) const noexcept -> bool { return this->compareTo(value) == Comparison::less; }
        gpu constexpr auto operator<=(const Aeu& value) const noexcept -> bool { return !this->operator>(value); }
        gpu constexpr auto operator>(const Aeu& value) const noexcept -> bool { return this->compareTo(value) == Comparison::greater; }
        gpu constexpr auto operator>=(const Aeu& value) const noexcept -> bool { return !this->operator<(value); }
#else
        gpu constexpr auto operator<=>(const Aeu& other) const noexcept -> std::strong_ordering {
            switch(this->compareTo(other)) {
                using enum Comparison;
                case less:
                    return std::strong_ordering::less;
                case greater:
                    return std::strong_ordering::greater;
                case equal:
                    return std::strong_ordering::equal;
                default:
                    return std::strong_ordering::equivalent;
            }
        }
 
        template <typename Unsigned>
        gpu constexpr auto operator<=>(const Unsigned& other) const noexcept -> std::strong_ordering {
            switch(this->compareTo(other)) {
                using enum Comparison;
                case less:
                    return std::strong_ordering::less;
                case greater:
                    return std::strong_ordering::greater;
                case equal:
                    return std::strong_ordering::equal;
                default:
                    return std::strong_ordering::equivalent;
            }
        }
#endif
    /* --------------------------------------------------------------------------- */
    /* ----------------------------------------------------------------------- */
 
 
    /* ---------------------- @name Supporting methods. ---------------------- */
    gpu constexpr auto setBit(std::size_t index, bool bit) noexcept -> void {
#ifndef AESI_UNSAFE
        if(index >= bitness) return;
#endif
        const std::size_t blockNumber = index / blockBitLength, bitNumber = index % blockBitLength;
        assert(blockNumber < blocksNumber && bitNumber < blockBitLength);
        if(bit)
            blocks[blockNumber] |= (1U << bitNumber);
        else
            blocks[blockNumber] &= (~(1U << bitNumber));
    }
 
    [[nodiscard]]
    gpu constexpr auto getBit(std::size_t index) const noexcept -> bool {
#ifndef AESI_UNSAFE
        if(index >= bitness) return false;
#endif
        const std::size_t blockNumber = index / blockBitLength, bitNumber = index % blockBitLength;
        assert(blockNumber < blocksNumber && bitNumber < blockBitLength);
        return blocks[blockNumber] & (1U << bitNumber);
    }
 
    gpu constexpr auto setByte(std::size_t index, byte byte) noexcept -> void {
#ifndef AESI_UNSAFE
        if(index > blocksNumber * sizeof(block)) return;
#endif
        const std::size_t blockNumber = index / sizeof(block), byteInBlock = index % sizeof(block), shift = byteInBlock * bitsInByte;
        assert(blockNumber < blocksNumber && byteInBlock < sizeof(block));
        blocks[blockNumber] &= ~(0xffU << shift); blocks[blockNumber] |= static_cast<block>(byte) << shift;
    }
 
    [[nodiscard]]
    gpu constexpr auto getByte(std::size_t index) const noexcept -> byte {
#ifndef AESI_UNSAFE
        if(index > blocksNumber * sizeof(block)) return 0;
#endif
        const std::size_t blockNumber = index / sizeof(block), byteInBlock = index % sizeof(block), shift = byteInBlock * bitsInByte;
        assert(blockNumber < blocksNumber && byteInBlock < sizeof(block));
        return static_cast<byte>((blocks[blockNumber] & (0xffU << shift)) >> shift);
    }
 
    gpu constexpr auto setBlock(std::size_t index, block block) noexcept -> void {
#ifndef AESI_UNSAFE
        if(index >= blocksNumber) return;
#endif
        blocks[index] = block;
    }
 
    [[nodiscard]]
    gpu constexpr auto getBlock(std::size_t index) const noexcept -> block {
#ifndef AESI_UNSAFE
        if(index >= blocksNumber) return block();
#endif
        return blocks[index];
    }
 
    [[nodiscard]]
    gpu constexpr auto byteCount() const noexcept -> std::size_t {
        std::size_t lastBlock = blocksNumber - 1;
        for(; lastBlock > 0 && blocks[lastBlock] == 0; --lastBlock)
            ;
 
        for(std::size_t byteN = sizeof(block); byteN-- > 0;) {
            if((blocks[lastBlock] & (0xffU << (byteN * bitsInByte))) >> (byteN * bitsInByte))
                return lastBlock * sizeof(block) + byteN + 1;
        }
        return lastBlock * sizeof(block);
    }
 
    [[nodiscard]]
    gpu constexpr auto bitCount() const noexcept -> std::size_t {
        std::size_t lastBlock = blocksNumber - 1;
        for(; lastBlock > 0 && blocks[lastBlock] == 0; --lastBlock)
            ;
 
        for(std::size_t byteN = sizeof(block); byteN-- > 0;) {
            const auto byte = (blocks[lastBlock] & (0xffU << (byteN * bitsInByte))) >> (byteN * bitsInByte);
            if(!byte) continue;
 
            for(std::size_t bitN = bitsInByte; bitN-- > 0;) {
                if((byte & (0x1u << bitN)) >> bitN)
                    return (lastBlock * sizeof(block) + byteN) * bitsInByte + bitN + 1;
            }
            return ((lastBlock - 1) * sizeof(block) + byteN) * bitsInByte;
        }
        return lastBlock * sizeof(block);
    }
 
    [[nodiscard]]
    gpu constexpr auto isOdd() const noexcept -> bool { return (0x1 & blocks[0]) == 1; }
 
    [[nodiscard]]
    gpu constexpr auto isEven() const noexcept -> bool { return (0x1 & blocks[0]) == 0; }
 
    [[nodiscard]]
    gpu constexpr auto isZero() const noexcept -> bool { return filledBlocksNumber() == 0; }
 
    [[nodiscard]]
    gpu constexpr auto filledBlocksNumber() const noexcept -> std::size_t {
        for(std::size_t i = blocksNumber; i-- > 0;)
            if(blocks[i]) return i + 1;
        return 0;
    }
 
    [[nodiscard]]
    gpu static constexpr auto getBitness() noexcept -> std::size_t { return bitness; }
 
    [[nodiscard]]
    gpu static constexpr auto totalBlocksNumber() noexcept -> std::size_t { return blocksNumber; }
 
    gpu constexpr auto swap(Aeu& other) noexcept -> void {
        Aeu t = other; other.operator=(*this); this->operator=(t);
    }
    /* ----------------------------------------------------------------------- */
 
    
    /* -------------- @name Public arithmetic and number theory. ------------- */
    gpu static constexpr auto divide(const Aeu& number, const Aeu& divisor, Aeu& quotient, Aeu& remainder) noexcept -> void {
        const auto ratio = number.compareTo(divisor);
 
        quotient = Aeu {}; remainder = Aeu {};
 
        if(ratio == Comparison::greater) {
            const auto bitsUsed = number.filledBlocksNumber() * blockBitLength;
            for(std::size_t i = bitsUsed; i-- > 0;) {
                remainder <<= 1u;
                remainder.setBit(0, number.getBit(i));
 
                if(remainder >= divisor) {
                    remainder -= divisor;
                    quotient.setBit(i, true);
                }
            }
        } else if(ratio == Comparison::less)
            remainder = number; else quotient = 1u;
    }
 
    [[nodiscard]]
    gpu static constexpr auto divide(const Aeu& number, const Aeu& divisor) noexcept -> pair<Aeu, Aeu> {
        pair<Aeu, Aeu> results; divide(number, divisor, results.first, results.second); return results;
    }
    
    gpu static constexpr auto gcd(const Aeu& first, const Aeu& second, Aeu& bezoutX, Aeu& bezoutY) noexcept -> Aeu {
        Aeu gcd, quotient, remainder;
 
        const auto ratio = first.compareTo(second);
        if(ratio == Comparison::greater) {
            gcd = second;
            divide(first, second, quotient, remainder);
        } else {
            gcd = first;
            divide(second, first, quotient, remainder);
        }
 
        bezoutX = 0u; bezoutY = 1u;
        for(Aeu tX = 1u, tY = 0u; remainder != 0u; ) {
            Aeu tGcd = gcd; gcd = remainder;
 
            Aeu t = bezoutX; bezoutX = tX - quotient * bezoutX; tX = t;
            t = bezoutY; bezoutY = tY - quotient * bezoutY; tY = t;
 
            divide(tGcd, gcd, quotient, remainder);
        }
 
        if(ratio != Comparison::greater)
            bezoutX.swap(bezoutY);
        
        return gcd;
    }
 
    gpu static constexpr auto gcd(const Aeu& first, const Aeu& second) noexcept -> Aeu {
        Aeu gcd, quotient, remainder;
 
        const auto ratio = first.compareTo(second);
        if(ratio == Comparison::greater) {
            gcd = second;
            divide(first, second, quotient, remainder);
        } else {
            gcd = first;
            divide(second, first, quotient, remainder);
        }
 
        while(remainder != 0u) {
            Aeu tGcd = gcd; gcd = remainder;
            divide(tGcd, gcd, quotient, remainder);
        }
 
        return gcd;
    }
 
    [[nodiscard]]
    gpu static constexpr auto lcm(const Aeu& first, const Aeu& second) noexcept -> Aeu { return first / gcd(first, second) * second; }
    
    template <std::size_t powerBitness = bitness> [[nodiscard]]
    gpu static constexpr auto powm(const Aeu& base, const Aeu<powerBitness>& power, const Aeu& modulo) noexcept -> Aeu {
        if(base == 1u)
            return base;
        if(base == 0u)
            return { 1u };
 
        Aeu output = 1u;
        auto [_, b] = divide(base, modulo);
 
        for(unsigned iteration = 0; power.filledBlocksNumber() * blockBitLength != iteration; iteration++) {
            if(power.getBit(iteration)) {
                const auto [quotient, remainder] = divide(output * b, modulo);
                output = remainder;
            }
 
            const auto [quotient, remainder] = divide(b * b, modulo);
            b = remainder;
        }
 
        return output;
    }
 
    [[nodiscard]]
    gpu static constexpr auto powm(const Aeu& base, const Aeu& power, const Aeu& mod) noexcept -> Aeu {
        Aeu output = 1u;
        auto [_, b] = divide(base, mod);
 
        for(std::size_t iteration = 0; power.filledBlocksNumber() * blockBitLength != iteration; iteration++) {
            if(power.getBit(iteration)) {
                const auto [quotient, remainder] = divide(output * b, mod);
                output = remainder;
            }
 
            const auto [quotient, remainder] = divide(b * b, mod);
            b = remainder;
        }
 
        return output;
    }
 
    [[nodiscard]]
    gpu static constexpr auto power2(std::size_t power) noexcept -> Aeu { Aeu result {}; result.setBit(power, true); return result; }
    
    [[nodiscard]]
    gpu constexpr auto squareRoot() const noexcept -> Aeu {
        Aeu x, y = power2((bitCount() + 1) / 2);
 
        do {
            x = y;
            y = (x + *this / x) >> 1u;
        } while (y < x);
 
        return x;
    }
    /* ----------------------------------------------------------------------- */
 
 
    /* ----------------- @name Public input-output operators. ---------------- */
    template <byte base, bool hexUppercase = false, typename Char> requires (std::is_same_v<Char, char> || (std::is_same_v<Char, wchar_t> && (base == 2 || base == 8 || base == 10 || base == 16)))
    gpu constexpr auto getString(Char* const buffer, std::size_t bufferSize, bool showBase = false) const noexcept -> std::size_t {
        if(bufferSize < 2) return 0;
 
        std::size_t position = 0;
 
        if (showBase && bufferSize > 3) {
            if constexpr (base == 2) {
                if constexpr (std::is_same_v<Char, char>) {
                    buffer[0] = '0'; buffer[1] = 'b';
                } else {
                    buffer[0] = L'0'; buffer[1] = L'b';
                }
                position += 2;
            } else if constexpr (base == 8) {
                if constexpr (std::is_same_v<Char, char>) {
                    buffer[0] = '0'; buffer[1] = 'o';
                } else {
                    buffer[0] = L'0'; buffer[1] = L'o';
                }
                position += 2;
            } else if constexpr (base == 16) {
                if constexpr (std::is_same_v<Char, char>) {
                    buffer[0] = '0'; buffer[1] = 'x';
                } else {
                    buffer[0] = L'0'; buffer[1] = L'x';
                }
                position += 2;
            }
        }
 
        if(isZero()) {
            buffer[position++] = [] { if constexpr (std::is_same_v<Char, char>) { return '0'; } else { return L'0'; }}();
            return position;
        }
 
        if constexpr (base == 16) {
            std::size_t iter = blocks.size() - 1;
            for (; iter > 0 && blocks[iter] == 0; --iter)
                ;
 
            if constexpr (std::is_same_v<Char, char>) {
                if constexpr (hexUppercase) {
                    position += static_cast<std::size_t>(snprintf(buffer + position, bufferSize - position, "%X", blocks[iter]));
                    for (; iter-- > 0;)
                        position += static_cast<std::size_t>(snprintf(buffer + position, bufferSize - position, "%08X", blocks[iter]));
                } else {
                    position += static_cast<std::size_t>(snprintf(buffer + position, bufferSize - position, "%x", blocks[iter]));
                    for (; iter-- > 0;)
                        position += static_cast<std::size_t>(snprintf(buffer + position, bufferSize - position, "%08x", blocks[iter]));
                }
            } else {
                if constexpr (hexUppercase) {
                    position += static_cast<std::size_t>(swprintf(buffer + position, bufferSize - position, L"%X", blocks[iter]));
                    for (; iter-- > 0;)
                        position += static_cast<std::size_t>(swprintf(buffer + position, bufferSize - position, L"%08X", blocks[iter]));
                } else {
                    position += static_cast<std::size_t>(swprintf(buffer + position, bufferSize - position, L"%x", blocks[iter]));
                    for (; iter-- > 0;)
                        position += static_cast<std::size_t>(swprintf(buffer + position, bufferSize - position, L"%08x", blocks[iter]));
                }
            }
        } else {
            const auto startPosition = position;
 
            Aeu copy = *this;
            while (!copy.isZero() && position < bufferSize) {
                auto [quotient, remainder] = divide(copy, base);
                if constexpr (std::is_same_v<Char, char>) {
                    buffer[position++] = static_cast<Char>('0' + remainder.template integralCast<byte>());
                } else {
                    buffer[position++] = static_cast<Char>(L'0' + remainder.template integralCast<byte>());
                }
                copy = quotient;
            }
            const auto digitsTotal = position - startPosition;
            for (std::size_t i = 0; i * 2 < digitsTotal; ++i) {
                Char t = buffer[startPosition + i];
                buffer[startPosition + i] = buffer[startPosition + digitsTotal - 1 - i];
                buffer[startPosition + digitsTotal - 1 - i] = t;
            }
        }
        buffer[position++] = Char {};
        return position;
    }
 
    template <typename Char> requires (std::is_same_v<Char, char> || std::is_same_v<Char, wchar_t>)
    friend constexpr auto operator<<(std::basic_ostream<Char>& os, const Aeu& number) -> std::basic_ostream<Char>& {
        auto flags = os.flags();
 
        const auto base = [] (long baseField, std::basic_ostream<Char>& ss, bool showbase) {
            auto tBase = (baseField == std::ios::hex ? 16u : (baseField == std::ios::oct ? 8u : 10u));
            if(showbase && tBase != 10)
                ss << [&tBase] {
                    if constexpr (std::is_same_v<Char, char>) {
                        return tBase == 8 ? "0o" : "0x";
                    } else {
                        return tBase == 8 ? L"0o" : L"0x";
                    }
                } () << std::noshowbase ;
            return tBase;
        } (flags & std::ios::basefield, os, flags & std::ios::showbase);
 
        if(number.isZero())
            return os << '0';
 
        if(base == 16) {
            std::size_t iter = number.blocks.size() - 1;
            for(; iter > 0 && number.blocks[iter] == 0; --iter)
                ;
 
            os << number.blocks[iter];
            for (; iter-- > 0;) {
                os.fill([] { if constexpr (std::is_same_v<Char, char>) { return '0'; } else { return L'0'; } } ());
                os.width(8); os << std::right << number.blocks[iter];
            }
        } else {
            /* Well, here we use a pre-calculated magic number to ratio the length of numbers in decimal or octal notation according to bitness.
             * * It is 2.95-98 for octal and 3.2 for decimal. */
            constexpr auto bufferSize = static_cast<std::size_t>(static_cast<double>(bitness) / 2.95);
            Char buffer [bufferSize] {}; std::size_t filled = 0;
 
            Aeu copy = number;
            while(!copy.isZero() && filled < bufferSize) {
                const auto [quotient, remainder] = divide(copy, base);
                buffer[filled++] = static_cast<Char>([] { if constexpr (std::is_same_v<Char, char>) { return '0'; } else { return L'0'; } } () + remainder.template integralCast<byte>());
                copy = quotient;
            }
 
            for(; filled > 0; --filled)
                os << buffer[filled - 1];
        }
 
        return os;
    }
 
 
    template <typename Char> requires (std::is_same_v<Char, char> || std::is_same_v<Char, wchar_t>)
    constexpr auto readBinary(std::basic_istream<Char>& is, bool bigEndian = true) -> void {
        blocks = {};
        if(bigEndian) {
            for(auto it = blocks.rbegin(); it != blocks.rend(); ++it)
                if(!is.read(reinterpret_cast<char*>(&*it), sizeof(block))) break;
        } else {
            for(auto& tBlock: blocks)
                if(!is.read(reinterpret_cast<char*>(&tBlock), sizeof(block))) break;
        }
    }
 
 
    template <typename Char> requires (std::is_same_v<Char, char> || std::is_same_v<Char, wchar_t>)
    constexpr auto writeBinary(std::basic_ostream<Char>& os, bool bigEndian = true) const noexcept -> void {
        if(bigEndian) {
            for(auto it = blocks.rbegin(); it != blocks.rend(); ++it)
                if(!os.write(reinterpret_cast<const char*>(&*it), sizeof(block))) break;
        } else {
            for(auto& block: blocks)
                if(!os.write(reinterpret_cast<const char*>(&block), sizeof(block))) break;
        }
    }
    /* ----------------------------------------------------------------------- */
 
 
    /* -------------------- @name Public casting operators. ------------------ */
    template <typename Integral> requires (std::is_integral_v<Integral>) [[nodiscard]]
    gpu constexpr auto integralCast() const noexcept -> Integral {
        const uint64_t value = (static_cast<uint64_t>(blocks[1]) << blockBitLength) | static_cast<uint64_t>(blocks[0]);
        return static_cast<Integral>(value);
    }
 
    template <std::size_t newBitness> requires (newBitness != bitness) [[nodiscard]]
    gpu constexpr auto precisionCast() const noexcept -> Aeu<newBitness> {
        Aeu<newBitness> result {};
 
        const std::size_t blockBoarder = (newBitness > bitness ? Aeu<bitness>::totalBlocksNumber() : Aeu<newBitness>::totalBlocksNumber());
        for(std::size_t blockIdx = 0; blockIdx < blockBoarder; ++blockIdx)
            result.setBlock(blockIdx, getBlock(blockIdx));
 
        return result;
    }
    /* ----------------------------------------------------------------------- */
 
 
#if defined __CUDACC__
    /* ------------------- @name Atomic-like CUDA operators. ----------------- */
    __device__ constexpr auto tryAtomicSet(const Aeu& value) noexcept -> void {
        for(std::size_t i = 0; i < blocksNumber; ++i)
            atomicExch(&blocks[i], value.blocks[i]);
    }
 
    __device__ constexpr auto tryAtomicExchange(const Aeu& value) noexcept -> void {
        for(std::size_t i = 0; i < blocksNumber; ++i)
            atomicExch(&value.blocks[i], atomicExch(&blocks[i], value.blocks[i]));
    }
    /* ----------------------------------------------------------------------- */
#endif
};
 
/* -------------------------------------------- @name Type-definitions  ------------------------------------------- */
using Aeu128 = Aeu<128>;
 
using Aeu256 = Aeu<256>;
 
using Aeu512 = Aeu<512>;
 
using Aeu768 = Aeu<768>;
 
using Aeu1024 = Aeu<1024>;
 
using Aeu1536 = Aeu<1536>;
 
using Aeu2048 = Aeu<2048>;
 
using Aeu3072 = Aeu<3072>;
 
using Aeu4096 = Aeu<4096>;
 
using Aeu6144 = Aeu<6144>;
 
using Aeu8192 = Aeu<8192>;
/* ---------------------------------------------------------------------------------------------------------------- */
 
/* ------------------------------------------ @name Integral conversions  ----------------------------------------- */
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator+(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) + value; }
 
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator-(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) - value; }
 
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator*(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) * value; }
 
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator/(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) / value; }
 
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator%(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) % value; }
 
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator^(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) ^ value; }
 
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator&(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) & value; }
 
template <std::size_t bitness, typename Integral> requires (std::is_integral_v<Integral>)
gpu constexpr auto operator|(Integral number, const Aeu<bitness>& value) noexcept { return Aeu<bitness>(number) | value; }
/* ---------------------------------------------------------------------------------------------------------------- */
 
#endif //AEU_MULTIPRECISION