numeric.hpp

Go to the documentation of this file.

// Copyright (c) 2026 Judica, Inc.
// Distributed under the MIT software license, see the accompanying
// file COPYING or https://opensource.org/license/mit/.
 
#pragma once
 
#include "purify/common.hpp"
#include "purify/uint.h"
 
namespace purify {
 
namespace detail {
 
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
#if defined(__clang__) || defined(__GNUC__)
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wpedantic"
#endif
#endif
 
// `UInt128` intentionally uses the compiler `__int128` extension when it is
// available; the fallback branch below remains the standards-only path.
class UInt128 {
public:
    [[nodiscard]] static UInt128 from_words(std::uint64_t hi, std::uint64_t lo) {
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
        return UInt128((static_cast<unsigned __int128>(hi) << 64) | lo);
#else
        return UInt128(hi, lo);
#endif
    }
 
    [[nodiscard]] static UInt128 mul_u64(std::uint64_t lhs, std::uint64_t rhs) {
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
        return UInt128(static_cast<unsigned __int128>(lhs) * rhs);
#else
        constexpr std::uint64_t kMask32 = 0xffffffffULL;
 
        std::uint64_t lhs_lo = lhs & kMask32;
        std::uint64_t lhs_hi = lhs >> 32;
        std::uint64_t rhs_lo = rhs & kMask32;
        std::uint64_t rhs_hi = rhs >> 32;
 
        std::uint64_t lo_lo = lhs_lo * rhs_lo;
        std::uint64_t hi_lo = lhs_hi * rhs_lo;
        std::uint64_t lo_hi = lhs_lo * rhs_hi;
        std::uint64_t hi_hi = lhs_hi * rhs_hi;
 
        std::uint64_t cross = (lo_lo >> 32) + (hi_lo & kMask32) + (lo_hi & kMask32);
        return UInt128(hi_hi + (hi_lo >> 32) + (lo_hi >> 32) + (cross >> 32),
                       (lo_lo & kMask32) | (cross << 32));
#endif
    }
 
    [[nodiscard]] UInt128 add_u64(std::uint64_t addend) const {
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
        return UInt128(value_ + addend);
#else
        UInt128 out = *this;
        out.lo_ += addend;
        out.hi_ += (out.lo_ < addend) ? 1U : 0U;
        return out;
#endif
    }
 
    [[nodiscard]] std::pair<UInt128, std::uint32_t> divmod_u32(std::uint32_t divisor) const {
        assert(divisor != 0 && "UInt128::divmod_u32 requires a non-zero divisor");
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
        unsigned __int128 quotient = value_ / divisor;
        return std::make_pair(UInt128(quotient), static_cast<std::uint32_t>(value_ % divisor));
#else
        constexpr std::uint64_t kMask32 = 0xffffffffULL;
 
        std::uint64_t rem = 0;
        auto step = [&](std::uint32_t word) {
            std::uint64_t cur = (rem << 32) | word;
            std::uint32_t qword = static_cast<std::uint32_t>(cur / divisor);
            rem = cur % divisor;
            return qword;
        };
 
        // `step()` mutates `rem`, so sequence each word explicitly instead of
        // relying on constructor-argument evaluation order, which differs across
        // compilers and broke the MSVC fallback decimal conversion path.
        const std::uint32_t q3 = step(static_cast<std::uint32_t>(hi_ >> 32));
        const std::uint32_t q2 = step(static_cast<std::uint32_t>(hi_ & kMask32));
        const std::uint32_t q1 = step(static_cast<std::uint32_t>(lo_ >> 32));
        const std::uint32_t q0 = step(static_cast<std::uint32_t>(lo_ & kMask32));
 
        UInt128 quotient((static_cast<std::uint64_t>(q3) << 32) | q2, (static_cast<std::uint64_t>(q1) << 32) | q0);
        return std::make_pair(quotient, static_cast<std::uint32_t>(rem));
#endif
    }
 
    [[nodiscard]] std::uint64_t low64() const {
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
        return static_cast<std::uint64_t>(value_);
#else
        return lo_;
#endif
    }
 
    [[nodiscard]] std::uint64_t high64() const {
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
        return static_cast<std::uint64_t>(value_ >> 64);
#else
        return hi_;
#endif
    }
 
private:
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
    explicit UInt128(unsigned __int128 value) : value_(value) {}
    unsigned __int128 value_ = 0;
#else
    UInt128(std::uint64_t hi, std::uint64_t lo) : lo_(lo), hi_(hi) {}
    std::uint64_t lo_ = 0;
    std::uint64_t hi_ = 0;
#endif
};
 
#if defined(__SIZEOF_INT128__) && !defined(_MSC_VER)
#if defined(__clang__) || defined(__GNUC__)
#pragma GCC diagnostic pop
#endif
#endif
 
[[nodiscard]] inline std::size_t bit_length_u64(std::uint64_t value) {
    return value == 0 ? 0U : 64U - static_cast<std::size_t>(std::countl_zero(value));
}
 
template <std::size_t Words>
struct BigUIntCBridge {
    static constexpr bool kAvailable = false;
};
 
#define PURIFY_DEFINE_BIGUINT_C_BRIDGE(words, bits) \
template <> \
struct BigUIntCBridge<words> { \
    static constexpr bool kAvailable = true; \
    static void set_zero(std::uint64_t* out) { purify_u##bits##_set_zero(out); } \
    static void set_u64(std::uint64_t* out, std::uint64_t value) { purify_u##bits##_set_u64(out, value); } \
    static void from_bytes_be(std::uint64_t* out, const unsigned char* data, std::size_t size) { \
        purify_u##bits##_from_bytes_be(out, data, size); \
    } \
    static bool is_zero(const std::uint64_t* value) { return purify_u##bits##_is_zero(value) != 0; } \
    static int compare(const std::uint64_t* lhs, const std::uint64_t* rhs) { return purify_u##bits##_compare(lhs, rhs); } \
    static bool try_add_small(std::uint64_t* value, std::uint32_t addend) { \
        return purify_u##bits##_try_add_small(value, addend) != 0; \
    } \
    static bool try_mul_small(std::uint64_t* value, std::uint32_t factor) { \
        return purify_u##bits##_try_mul_small(value, factor) != 0; \
    } \
    static bool try_add(std::uint64_t* value, const std::uint64_t* addend) { \
        return purify_u##bits##_try_add(value, addend) != 0; \
    } \
    static bool try_sub(std::uint64_t* value, const std::uint64_t* subtrahend) { \
        return purify_u##bits##_try_sub(value, subtrahend) != 0; \
    } \
    static std::size_t bit_length(const std::uint64_t* value) { return purify_u##bits##_bit_length(value); } \
    static bool bit(const std::uint64_t* value, std::size_t index) { return purify_u##bits##_bit(value, index) != 0; } \
    static bool try_set_bit(std::uint64_t* value, std::size_t index) { \
        return purify_u##bits##_try_set_bit(value, index) != 0; \
    } \
    static void shifted_left(std::uint64_t* out, const std::uint64_t* value, std::size_t bits_value) { \
        purify_u##bits##_shifted_left(out, value, bits_value); \
    } \
    static void shifted_right(std::uint64_t* out, const std::uint64_t* value, std::size_t bits_value) { \
        purify_u##bits##_shifted_right(out, value, bits_value); \
    } \
    static void shift_right_one(std::uint64_t* value) { purify_u##bits##_shift_right_one(value); } \
    static void mask_bits(std::uint64_t* value, std::size_t bits_value) { purify_u##bits##_mask_bits(value, bits_value); } \
    static std::uint32_t divmod_small(std::uint64_t* value, std::uint32_t divisor) { \
        return purify_u##bits##_divmod_small(value, divisor); \
    } \
    static void to_bytes_be(unsigned char* out, const std::uint64_t* value) { purify_u##bits##_to_bytes_be(out, value); } \
}
 
PURIFY_DEFINE_BIGUINT_C_BRIDGE(4, 256);
PURIFY_DEFINE_BIGUINT_C_BRIDGE(5, 320);
PURIFY_DEFINE_BIGUINT_C_BRIDGE(8, 512);
 
#undef PURIFY_DEFINE_BIGUINT_C_BRIDGE
 
struct FieldElementAccess;
 
}  // namespace detail
 
template <std::size_t Words>
struct BigUInt {
    std::array<std::uint64_t, Words> limbs{};
 
    static BigUInt zero() {
        BigUInt out;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::set_zero(out.limbs.data());
        }
        return out;
    }
 
    static BigUInt one() {
        BigUInt out;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::set_u64(out.limbs.data(), 1);
        } else {
            out.limbs[0] = 1;
        }
        return out;
    }
 
    static BigUInt from_u64(std::uint64_t value) {
        BigUInt out;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::set_u64(out.limbs.data(), value);
        } else {
            out.limbs[0] = value;
        }
        return out;
    }
 
    static BigUInt from_bytes_be(const unsigned char* data, std::size_t size) {
        BigUInt out;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::from_bytes_be(out.limbs.data(), data, size);
        } else {
            for (std::size_t i = 0; i < size; ++i) {
                out.mul_small(256);
                out.add_small(data[i]);
            }
        }
        return out;
    }
 
    static Result<BigUInt> try_from_hex(std::string_view hex) {
        BigUInt out;
        if (hex.size() >= 2 && hex[0] == '0' && (hex[1] == 'x' || hex[1] == 'X')) {
            hex.remove_prefix(2);
        }
        for (char ch : hex) {
            if (std::isspace(static_cast<unsigned char>(ch)) != 0) {
                continue;
            }
            unsigned value;
            if (ch >= '0' && ch <= '9') {
                value = static_cast<unsigned>(ch - '0');
            } else if (ch >= 'a' && ch <= 'f') {
                value = static_cast<unsigned>(10 + ch - 'a');
            } else if (ch >= 'A' && ch <= 'F') {
                value = static_cast<unsigned>(10 + ch - 'A');
            } else {
                return unexpected_error(ErrorCode::InvalidHex, "BigUInt::try_from_hex:invalid_digit");
            }
            if (!out.try_mul_small(16) || !out.try_add_small(value)) {
                return unexpected_error(ErrorCode::Overflow, "BigUInt::try_from_hex:overflow");
            }
        }
        return out;
    }
 
    static BigUInt from_hex(std::string_view hex) {
        Result<BigUInt> out = try_from_hex(hex);
        assert(out.has_value() && "BigUInt::from_hex() requires valid in-range input");
        return std::move(*out);
    }
 
    bool is_zero() const {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            return detail::BigUIntCBridge<Words>::is_zero(limbs.data());
        } else {
            for (std::uint64_t limb : limbs) {
                if (limb != 0) {
                    return false;
                }
            }
            return true;
        }
    }
 
    int compare(const BigUInt& other) const {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            return detail::BigUIntCBridge<Words>::compare(limbs.data(), other.limbs.data());
        } else {
            for (std::size_t i = Words; i != 0; --i) {
                const std::size_t idx = i - 1;
                if (limbs[idx] < other.limbs[idx]) {
                    return -1;
                }
                if (limbs[idx] > other.limbs[idx]) {
                    return 1;
                }
            }
            return 0;
        }
    }
 
    bool operator==(const BigUInt& other) const {
        return limbs == other.limbs;
    }
 
    bool operator!=(const BigUInt& other) const {
        return !(*this == other);
    }
 
    bool operator<(const BigUInt& other) const {
        return compare(other) < 0;
    }
 
    bool operator>=(const BigUInt& other) const {
        return compare(other) >= 0;
    }
 
    bool try_add_small(std::uint32_t value) {
        BigUInt out = *this;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            if (!detail::BigUIntCBridge<Words>::try_add_small(out.limbs.data(), value)) {
                return false;
            }
        } else {
            std::uint64_t carry = value;
            for (std::size_t i = 0; i < Words && carry != 0; ++i) {
                std::uint64_t sum = out.limbs[i] + carry;
                carry = (sum < out.limbs[i]) ? 1U : 0U;
                out.limbs[i] = sum;
            }
            if (carry != 0) {
                return false;
            }
        }
        *this = out;
        return true;
    }
 
    void add_small(std::uint32_t value) {
        bool ok = try_add_small(value);
        assert(ok && "BigUInt::add_small() overflow");
        (void)ok;
    }
 
    bool try_mul_small(std::uint32_t value) {
        BigUInt out = *this;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            if (!detail::BigUIntCBridge<Words>::try_mul_small(out.limbs.data(), value)) {
                return false;
            }
        } else {
            std::uint64_t carry = 0;
            for (std::size_t i = 0; i < Words; ++i) {
                detail::UInt128 accum = detail::UInt128::mul_u64(out.limbs[i], value);
                accum = accum.add_u64(carry);
                out.limbs[i] = accum.low64();
                carry = accum.high64();
            }
            if (carry != 0) {
                return false;
            }
        }
        *this = out;
        return true;
    }
 
    void mul_small(std::uint32_t value) {
        bool ok = try_mul_small(value);
        assert(ok && "BigUInt::mul_small() overflow");
        (void)ok;
    }
 
    bool try_add_assign(const BigUInt& other) {
        BigUInt out = *this;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            if (!detail::BigUIntCBridge<Words>::try_add(out.limbs.data(), other.limbs.data())) {
                return false;
            }
        } else {
            std::uint64_t carry = 0;
            for (std::size_t i = 0; i < Words; ++i) {
                std::uint64_t sum = out.limbs[i] + other.limbs[i];
                std::uint64_t carry1 = (sum < out.limbs[i]) ? 1U : 0U;
                std::uint64_t next = sum + carry;
                std::uint64_t carry2 = (next < sum) ? 1U : 0U;
                out.limbs[i] = next;
                carry = carry1 | carry2;
            }
            if (carry != 0) {
                return false;
            }
        }
        *this = out;
        return true;
    }
 
    void add_assign(const BigUInt& other) {
        bool ok = try_add_assign(other);
        assert(ok && "BigUInt::add_assign() overflow");
        (void)ok;
    }
 
    bool try_sub_assign(const BigUInt& other) {
        BigUInt out = *this;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            if (!detail::BigUIntCBridge<Words>::try_sub(out.limbs.data(), other.limbs.data())) {
                return false;
            }
        } else {
            std::uint64_t borrow = 0;
            for (std::size_t i = 0; i < Words; ++i) {
                std::uint64_t rhs = other.limbs[i] + borrow;
                std::uint64_t rhs_overflow = (rhs < other.limbs[i]) ? 1U : 0U;
                std::uint64_t next = out.limbs[i] - rhs;
                std::uint64_t needs_borrow = (out.limbs[i] < rhs) ? 1U : 0U;
                out.limbs[i] = next;
                borrow = rhs_overflow | needs_borrow;
            }
            if (borrow != 0) {
                return false;
            }
        }
        *this = out;
        return true;
    }
 
    void sub_assign(const BigUInt& other) {
        bool ok = try_sub_assign(other);
        assert(ok && "BigUInt::sub_assign() underflow");
        (void)ok;
    }
 
    std::size_t bit_length() const {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            return detail::BigUIntCBridge<Words>::bit_length(limbs.data());
        } else {
            for (std::size_t i = Words; i != 0; --i) {
                const std::size_t idx = i - 1;
                if (limbs[idx] != 0) {
                    return idx * 64 + detail::bit_length_u64(limbs[idx]);
                }
            }
            return 0;
        }
    }
 
    bool bit(std::size_t index) const {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            return detail::BigUIntCBridge<Words>::bit(limbs.data(), index);
        } else {
            std::size_t word = index / 64;
            std::size_t shift = index % 64;
            if (word >= Words) {
                return false;
            }
            return ((limbs[word] >> shift) & 1U) != 0;
        }
    }
 
    bool try_set_bit(std::size_t index) {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            return detail::BigUIntCBridge<Words>::try_set_bit(limbs.data(), index);
        } else {
            std::size_t word = index / 64;
            std::size_t shift = index % 64;
            if (word >= Words) {
                return false;
            }
            limbs[word] |= (static_cast<std::uint64_t>(1) << shift);
            return true;
        }
    }
 
    void set_bit(std::size_t index) {
        bool ok = try_set_bit(index);
        assert(ok && "BigUInt::set_bit() index out of range");
        (void)ok;
    }
 
    BigUInt shifted_left(std::size_t bits) const {
        BigUInt out;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::shifted_left(out.limbs.data(), limbs.data(), bits);
        } else {
            std::size_t word_shift = bits / 64;
            std::size_t bit_shift = bits % 64;
            for (std::size_t i = Words; i != 0; --i) {
                const std::size_t idx = i - 1;
                if (idx < word_shift) {
                    continue;
                }
                std::size_t src = idx - word_shift;
                out.limbs[idx] |= limbs[src] << bit_shift;
                if (bit_shift != 0 && src > 0) {
                    out.limbs[idx] |= limbs[src - 1] >> (64 - bit_shift);
                }
            }
        }
        return out;
    }
 
    BigUInt shifted_right(std::size_t bits) const {
        BigUInt out;
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::shifted_right(out.limbs.data(), limbs.data(), bits);
        } else {
            std::size_t word_shift = bits / 64;
            std::size_t bit_shift = bits % 64;
            for (std::size_t i = 0; i < Words; ++i) {
                std::size_t src = i + word_shift;
                if (src >= Words) {
                    break;
                }
                out.limbs[i] |= limbs[src] >> bit_shift;
                if (bit_shift != 0 && src + 1 < Words) {
                    out.limbs[i] |= limbs[src + 1] << (64 - bit_shift);
                }
            }
        }
        return out;
    }
 
    void shift_right_one() {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::shift_right_one(limbs.data());
        } else {
            for (std::size_t i = 0; i < Words; ++i) {
                std::uint64_t next = (i + 1 < Words) ? limbs[i + 1] : 0;
                limbs[i] = (limbs[i] >> 1) | (next << 63);
            }
        }
    }
 
    void mask_bits(std::size_t bits) {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::mask_bits(limbs.data(), bits);
        } else {
            std::size_t full_words = bits / 64;
            std::size_t extra_bits = bits % 64;
            for (std::size_t i = full_words + (extra_bits != 0 ? 1 : 0); i < Words; ++i) {
                limbs[i] = 0;
            }
            if (extra_bits != 0 && full_words < Words) {
                std::uint64_t mask =
                    (extra_bits == 64) ? ~static_cast<std::uint64_t>(0) : ((static_cast<std::uint64_t>(1) << extra_bits) - 1);
                limbs[full_words] &= mask;
            }
        }
    }
 
    std::uint32_t divmod_small(std::uint32_t divisor) {
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            return detail::BigUIntCBridge<Words>::divmod_small(limbs.data(), divisor);
        } else {
            std::uint64_t rem = 0;
            for (std::size_t i = Words; i != 0; --i) {
                const std::size_t idx = i - 1;
                auto [quotient, next_rem] = detail::UInt128::from_words(rem, limbs[idx]).divmod_u32(divisor);
                assert(quotient.high64() == 0 && "BigUInt::divmod_small() quotient must fit in 64 bits");
                limbs[idx] = quotient.low64();
                rem = next_rem;
            }
            return static_cast<std::uint32_t>(rem);
        }
    }
 
    std::array<unsigned char, Words * 8> to_bytes_be() const {
        std::array<unsigned char, Words * 8> out{};
        if constexpr (detail::BigUIntCBridge<Words>::kAvailable) {
            detail::BigUIntCBridge<Words>::to_bytes_be(out.data(), limbs.data());
        } else {
            for (std::size_t i = 0; i < Words; ++i) {
                std::uint64_t limb = limbs[i];
                for (std::size_t j = 0; j < 8; ++j) {
                    out[out.size() - 1 - (i * 8 + j)] = static_cast<unsigned char>(limb & 0xffU);
                    limb >>= 8;
                }
            }
        }
        return out;
    }
 
    std::string to_hex() const {
        auto bytes = to_bytes_be();
        std::ostringstream out;
        bool started = false;
        out << std::hex << std::nouppercase;
        for (unsigned char byte : bytes) {
            if (!started) {
                if (byte == 0) {
                    continue;
                }
                out << static_cast<unsigned>(byte);
                started = true;
            } else {
                out << std::setw(2) << std::setfill('0') << static_cast<unsigned>(byte);
            }
        }
        return started ? out.str() : "0";
    }
 
    std::string to_decimal() const {
        if (is_zero()) {
            return "0";
        }
        BigUInt copy = *this;
        std::vector<std::uint32_t> parts;
        while (!copy.is_zero()) {
            parts.push_back(copy.divmod_small(1000000000U));
        }
        std::ostringstream out;
        out << parts.back();
        for (std::size_t i = parts.size(); i > 1; --i) {
            out << std::setw(9) << std::setfill('0') << parts[i - 2];
        }
        return out.str();
    }
};
 
template <std::size_t OutWords, std::size_t InWords>
BigUInt<OutWords> widen(const BigUInt<InWords>& value) {
    static_assert(OutWords >= InWords, "Cannot narrow with widen");
    BigUInt<OutWords> out;
    if constexpr (OutWords == 5 && InWords == 4) {
        purify_u320_widen_u256(out.limbs.data(), value.limbs.data());
    } else if constexpr (OutWords == 8 && InWords == 4) {
        purify_u512_widen_u256(out.limbs.data(), value.limbs.data());
    } else {
        for (std::size_t i = 0; i < InWords; ++i) {
            out.limbs[i] = value.limbs[i];
        }
    }
    return out;
}
 
template <std::size_t OutWords, std::size_t InWords>
Result<BigUInt<OutWords>> try_narrow(const BigUInt<InWords>& value) {
    static_assert(OutWords <= InWords, "Cannot widen with narrow");
    BigUInt<OutWords> out;
    if constexpr (OutWords == 4 && InWords == 5) {
        if (!purify_u256_try_narrow_u320(out.limbs.data(), value.limbs.data())) {
            return unexpected_error(ErrorCode::NarrowingOverflow, "try_narrow:high_bits_set");
        }
    } else if constexpr (OutWords == 4 && InWords == 8) {
        if (!purify_u256_try_narrow_u512(out.limbs.data(), value.limbs.data())) {
            return unexpected_error(ErrorCode::NarrowingOverflow, "try_narrow:high_bits_set");
        }
    } else {
        for (std::size_t i = 0; i < OutWords; ++i) {
            out.limbs[i] = value.limbs[i];
        }
        for (std::size_t i = OutWords; i < InWords; ++i) {
            if (value.limbs[i] != 0) {
                return unexpected_error(ErrorCode::NarrowingOverflow, "try_narrow:high_bits_set");
            }
        }
    }
    return out;
}
 
template <std::size_t OutWords, std::size_t InWords>
BigUInt<OutWords> narrow(const BigUInt<InWords>& value) {
    Result<BigUInt<OutWords>> out = try_narrow<OutWords>(value);
    assert(out.has_value() && "narrow() requires the source value to fit");
    return std::move(*out);
}
 
template <std::size_t Words>
Result<std::pair<BigUInt<Words>, BigUInt<Words>>> try_divmod_same(const BigUInt<Words>& numerator, const BigUInt<Words>& denominator) {
    BigUInt<Words> quotient;
    BigUInt<Words> remainder = numerator;
    if constexpr (Words == 8) {
        if (!purify_u512_try_divmod_same(quotient.limbs.data(), remainder.limbs.data(),
                                         numerator.limbs.data(), denominator.limbs.data())) {
            return unexpected_error(ErrorCode::DivisionByZero, "try_divmod_same:zero_denominator");
        }
    } else {
        if (denominator.is_zero()) {
            return unexpected_error(ErrorCode::DivisionByZero, "try_divmod_same:zero_denominator");
        }
        std::size_t n_bits = remainder.bit_length();
        std::size_t d_bits = denominator.bit_length();
        if (n_bits < d_bits) {
            return std::make_pair(quotient, remainder);
        }
        std::size_t shift = n_bits - d_bits;
        BigUInt<Words> shifted = denominator.shifted_left(shift);
        for (std::size_t i = shift + 1; i != 0; --i) {
            const std::size_t idx = i - 1;
            if (remainder.compare(shifted) >= 0) {
                bool sub_ok = remainder.try_sub_assign(shifted);
                bool bit_ok = quotient.try_set_bit(idx);
                assert(sub_ok && "divmod_same() subtraction should stay in range");
                assert(bit_ok && "divmod_same() quotient bit index should stay in range");
                if (!sub_ok || !bit_ok) {
                    return unexpected_error(ErrorCode::InternalMismatch, "try_divmod_same:internal_step");
                }
            }
            if (idx != 0) {
                shifted.shift_right_one();
            }
        }
    }
    return std::make_pair(quotient, remainder);
}
 
template <std::size_t Words>
std::pair<BigUInt<Words>, BigUInt<Words>> divmod_same(const BigUInt<Words>& numerator, const BigUInt<Words>& denominator) {
    Result<std::pair<BigUInt<Words>, BigUInt<Words>>> out = try_divmod_same(numerator, denominator);
    assert(out.has_value() && "divmod_same() requires a non-zero denominator");
    return std::move(*out);
}
 
template <std::size_t LeftWords, std::size_t RightWords>
BigUInt<LeftWords + RightWords> multiply(const BigUInt<LeftWords>& lhs, const BigUInt<RightWords>& rhs) {
    BigUInt<LeftWords + RightWords> out;
    if constexpr (LeftWords == 4 && RightWords == 4) {
        purify_u512_multiply_u256(out.limbs.data(), lhs.limbs.data(), rhs.limbs.data());
    } else {
        for (std::size_t i = 0; i < LeftWords; ++i) {
            std::uint64_t carry = 0;
            for (std::size_t j = 0; j < RightWords; ++j) {
                detail::UInt128 accum = detail::UInt128::mul_u64(lhs.limbs[i], rhs.limbs[j]);
                accum = accum.add_u64(out.limbs[i + j]);
                accum = accum.add_u64(carry);
                out.limbs[i + j] = accum.low64();
                carry = accum.high64();
            }
            out.limbs[i + RightWords] += carry;
        }
    }
    return out;
}
 
using UInt256 = BigUInt<4>;
using UInt320 = BigUInt<5>;
using UInt512 = BigUInt<8>;
 
class FieldElement;
 
const UInt256& prime_p();
 
class FieldElement {
public:
    FieldElement();
 
    static FieldElement zero();
 
    static FieldElement one();
 
    static FieldElement from_u64(std::uint64_t value);
 
    static FieldElement from_int(std::int64_t value);
 
    static Result<FieldElement> try_from_bytes32(const std::array<unsigned char, 32>& bytes);
 
    static FieldElement from_bytes32(const std::array<unsigned char, 32>& bytes);
 
    static Result<FieldElement> try_from_uint256(const UInt256& value);
 
    static FieldElement from_uint256(const UInt256& value);
 
    UInt256 to_uint256() const;
 
    std::array<unsigned char, 32> to_bytes_be() const;
 
    std::array<unsigned char, 32> to_bytes_le() const;
 
    std::string to_hex() const;
 
    std::string to_decimal() const;
 
    bool is_zero() const;
 
    bool is_one() const;
 
    bool is_odd() const;
 
    bool is_square() const;
 
    FieldElement negate() const;
 
    void conditional_assign(const FieldElement& other, bool flag);
 
    FieldElement inverse_consttime() const;
 
    FieldElement inverse() const;
 
    std::optional<FieldElement> sqrt() const;
 
    FieldElement pow(const UInt256& exponent) const;
 
    friend bool operator==(const FieldElement& lhs, const FieldElement& rhs);
 
    friend bool operator!=(const FieldElement& lhs, const FieldElement& rhs);
 
    friend FieldElement operator+(const FieldElement& lhs, const FieldElement& rhs);
 
    friend FieldElement operator-(const FieldElement& lhs, const FieldElement& rhs);
 
    friend FieldElement operator*(const FieldElement& lhs, const FieldElement& rhs);
 
private:
    friend struct detail::FieldElementAccess;
    explicit FieldElement(const purify_scalar& raw) : value_(raw) {}
 
    purify_scalar value_{};
};
 
namespace detail {
 
struct FieldElementAccess {
    static const purify_scalar& raw(const FieldElement& value) noexcept {
        return value.value_;
    }
 
    static FieldElement from_raw(const purify_scalar& raw) noexcept {
        return FieldElement(raw);
    }
};
 
}  // namespace detail
 
FieldElement square(const FieldElement& value);
 
int legendre_symbol(const FieldElement& value);
 
}  // namespace purify