Compiler Explorer

Source code

#include <cstdint>
#include <bit>
#include <stdexcept>
#include <string>

#include <immintrin.h>

#if defined(_MSC_VER) && !defined(__clang__)
    #define IS_MSVC 1
#else
    #define IS_MSVC 0
#endif

#if IS_MSVC
    #include <intrin.h>
#endif

#define ASSERT_MSG(cond, msg) { if (!(cond)) throw std::runtime_error("Assertion (" #cond ") failed at line " + std::to_string(__LINE__) + "! Msg: '" + std::string(msg) + "'."); }
#define ASSERT(cond) ASSERT_MSG(cond, "")
#ifdef _DEBUG
    #define DASSERT_MSG(cond, msg) ASSERT_MSG(cond, msg)
#else
    #define DASSERT_MSG(cond, msg)
#endif 
#define DASSERT(cond) DASSERT_MSG(cond, "")

using u16 = uint16_t;
using u32 = uint32_t;
using i64 = int64_t;
using u64 = uint64_t;
using UllPtr = unsigned long long *;

inline int GetExp(double x) {
    return int((std::bit_cast<uint64_t>(x) >> 52) & 0x7FF) - 1023;
}

inline size_t BitSizeWrong(uint64_t x) {
    return x == 0 ? 0 : (GetExp(x) + 1);
}

inline size_t BitSize(u64 x) {
    size_t r = 0;
    if (x >= (u64(1) << 32)) {
        x >>= 32;
        r += 32;
    }
    while (x >= 0x100) {
        x >>= 8;
        r += 8;
    }
    while (x) {
        x >>= 1;
        ++r;
    }
    return r;
}

#if !IS_MSVC
inline u64 UDiv128Asm(u64 h, u64 l, u64 d, u64 * r) {
    u64 q;
    asm (R"(
        .intel_syntax
        mov rdx, %V[h]
        mov rax, %V[l]
        div %V[d]
        mov %V[r], rdx
        mov %V[q], rax
    )"
        : [q] "=r" (q), [r] "=r" (*r)
        : [h] "r" (h), [l] "r" (l), [d] "r" (d)
        : "rax", "rdx"
    );
    return q;
}
#endif

inline std::pair<u64, u64> UDiv128(u64 hi, u64 lo, u64 d) {
    #if IS_MSVC
        u64 r, q = _udiv128(hi, lo, d, &r);
        return std::make_pair(q, r);
    #else
        #if 0
            using u128 = unsigned __int128;
            auto const dnd = (u128(hi) << 64) | lo;
            return std::make_pair(u64(dnd / d), u64(dnd % d));
        #else
            u64 r, q = UDiv128Asm(hi, lo, d, &r);
            return std::make_pair(q, r);
        #endif
    #endif
}

inline std::pair<u64, u64> UMul128(u64 a, u64 b) {
    #if IS_MSVC
        u64 hi, lo = _umul128(a, b, &hi);
        return std::make_pair(hi, lo);
    #else
        using u128 = unsigned __int128;
        auto const x = u128(a) * b;
        return std::make_pair(u64(x >> 64), u64(x));
    #endif
}

inline std::pair<u64, u64> USub128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
    u64 r_hi, r_lo;
    _subborrow_u64(_subborrow_u64(0, a_lo, b_lo, (UllPtr)&r_lo), a_hi, b_hi, (UllPtr)&r_hi);
    return std::make_pair(r_hi, r_lo);
}

inline std::pair<u64, u64> UAdd128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
    u64 r_hi, r_lo;
    _addcarry_u64(_addcarry_u64(0, a_lo, b_lo, (UllPtr)&r_lo), a_hi, b_hi, (UllPtr)&r_hi);
    return std::make_pair(r_hi, r_lo);
}

inline int UCmp128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
    if (a_hi != b_hi)
        return a_hi < b_hi ? -1 : 1;
    return a_lo == b_lo ? 0 : a_lo < b_lo ? -1 : 1;
}

std::pair<u64, size_t> BarrettRS64(u64 n) {
    // https://www.nayuki.io/page/barrett-reduction-algorithm
    ASSERT_MSG(n >= 3 && (n & (n - 1)) != 0, "n " + std::to_string(n))
    size_t const nbits = BitSize(n);
    // 2^s = q * n + r;   2^s = (2^64 + q0) * n + r;  2^s - n * 2^64 = q0 * n + r
    u64 const dnd_hi = (nbits >= 64 ? 0ULL : (u64(1) << nbits)) - n;
    auto const q0 = UDiv128(dnd_hi, 0, n).first;
    return std::make_pair(q0, nbits);
}

template <bool Use128 = true, bool Adjust = true>
std::pair<u64, u64> BarrettDivMod64(u64 x_hi, u64 x_lo, u64 n, u64 r, size_t s) {
    // ((x_hi * 2^64 + x_lo) * (2^64 + r)) >> (64 + s)
    
    DASSERT(x_hi < n);
    
    #if !IS_MSVC
    if constexpr(Use128) {
        using u128 = unsigned __int128;
        u128 const xf = (u128(x_hi) << 64) | x_lo;
        u64 q = u64((u128(x_hi) * r + xf + u64((u128(x_lo) * r) >> 64)) >> s);
        if (s < 64) {
            u64 t = x_lo - q * n;
            if constexpr(Adjust) {
                u64 const mask = ~u64(i64(t - n) >> 63);
                q += mask & 1;
                t -= mask & n;
            }
            return std::make_pair(q, t);
        } else {
            u128 t = xf - u128(q) * n;
            return t < n ? std::make_pair(q, u64(t)) : std::make_pair(q + 1, u64(t) - n);
        }
    } else
    #endif
    {
        auto const w1a = UMul128(x_lo, r).first;
        auto const [w2b, w1b] = UMul128(x_hi, r);
        auto const w2c = x_hi, w1c = x_lo;
        u64 w1, w2 = _addcarry_u64(0, w1a, w1b, (UllPtr)&w1);
        w2 += _addcarry_u64(0, w1, w1c, (UllPtr)&w1);
        w2 += w2b + w2c;
        
        if (s < 64) {
            u64 q = (w2 << (64 - s)) | (w1 >> s);
            u64 t = x_lo - q * n;
            if constexpr(Adjust) {
                u64 const mask = ~u64(i64(t - n) >> 63);
                q += mask & 1;
                t -= mask & n;
            }
            return std::make_pair(q, t);
        } else {
            u64 const q = w2;
            auto const [b_hi, b_lo] = UMul128(q, n);
            auto const [t_hi, t_lo] = USub128(x_hi, x_lo, b_hi, b_lo);
            return t_hi != 0 || t_lo >= n ? std::make_pair(q + 1, t_lo - n) : std::make_pair(q, t_lo);
        }
    }
}

#include <random>
#include <iomanip>
#include <iostream>
#include <chrono>

void TestBarrett() {
    std::mt19937_64 rng{123}; //{std::random_device{}()};
    for (size_t i = 0; i < (1 << 11); ++i) {
        size_t const nbits = rng() % 63 + 2;
        u64 n = 0;
        do {
            n = (u64(1) << (nbits - 1)) + rng() % (u64(1) << (nbits - 1));
        } while (!(n >= 3 && (n & (n - 1)) != 0));
        auto const [br, bs] = BarrettRS64(n);
        for (size_t j = 0; j < (1 << 6); ++j) {
            u64 const hi = rng() % n, lo = rng();
            auto const [ref_q, ref_r] = UDiv128(hi, lo, n);
            u64 bar_q = 0, bar_r = 0;
            for (size_t k = 0; k < 2; ++k) {
                bar_q = 0; bar_r = 0;
                if (k == 0)
                    std::tie(bar_q, bar_r) = BarrettDivMod64<true>(hi, lo, n, br, bs);
                else
                    std::tie(bar_q, bar_r) = BarrettDivMod64<false>(hi, lo, n, br, bs);
                ASSERT_MSG(bar_q == ref_q && bar_r == ref_r, "i " + std::to_string(i) + ", j " + std::to_string(j) + ", k " + std::to_string(k) +
                    ", nbits " + std::to_string(nbits) + ", n " + std::to_string(n) + ", bar_q " + std::to_string(bar_q) +
                    ", ref_q " + std::to_string(ref_q) + ", bar_r " + std::to_string(bar_r) + ", ref_r " + std::to_string(ref_r));
            }
        }
    }
}

class bi_uint
{
public:
    using u64_t = std::uint64_t;
    constexpr static std::size_t u64_bits = 8 * sizeof(u64_t);

//little-endian
    std::vector<u64_t> data;

static auto constexpr DefPrep = [](auto n){
        return std::make_pair(false, false);
    };
    static auto constexpr DefDivMod = [](auto dnd_hi, auto dnd_lo, auto dsr, auto const & prep){
        return UDiv128(dnd_hi, dnd_lo, dsr);
    };
    
    //User should guarantee data.size()>0 and val>0
    template <typename PrepT = decltype(DefPrep), typename DivModT = decltype(DefDivMod)>
    void self_div(const u64_t val, PrepT const & Prep = DefPrep, DivModT const & DivMod = DefDivMod)
    {
        auto it = data.rbegin();

if(data.size() == 1) {
            *it /= val;
            return;
        }    
        
        u64_t rem_hi = 0, rem_lo = 0;
        if(*it < val) {
            rem_lo = *it++;
            //data.pop_back();
        }
        
        auto const prep = Prep(val);
        
        u64_t r = rem_lo % val;
        u64_t q = 0;
        while(it != data.rend()) {
            rem_hi = r;
            rem_lo = *it;
            std::tie(q, r) = DivMod(rem_hi, rem_lo, val, prep);
            //*it++ = static_cast<u64_t>(q);
            it++;
            auto volatile out = static_cast<u64_t>(q);
        }
    }
};

void TestSpeed() {
    auto Time = []{
        static auto const gtb = std::chrono::high_resolution_clock::now();
        return std::chrono::duration_cast<std::chrono::duration<double>>(
            std::chrono::high_resolution_clock::now() - gtb).count();
    };
    std::mt19937_64 rng{123};
    std::vector<u64> limbs, divisors;
    for (size_t i = 0; i < (1 << 17); ++i)
        limbs.push_back(rng());
    for (size_t i = 0; i < (1 << 8); ++i) {
        size_t const nbits = rng() % 63 + 2;
        u64 n = 0;
        do {
            n = (u64(1) << (nbits - 1)) + rng() % (u64(1) << (nbits - 1));
        } while (!(n >= 3 && (n & (n - 1)) != 0));
        divisors.push_back(n);
    }
    std::cout << std::fixed << std::setprecision(3);
    double div_time = 0;
    {
        bi_uint x;
        x.data = limbs;
        auto const tim = Time();
        for (auto dsr: divisors)
            x.self_div(dsr);
        div_time = Time() - tim;
        std::cout << "Divide time " << div_time << " sec" << std::endl;
    }
    {
        bi_uint x;
        x.data = limbs;
        for (size_t i = 0; i < 2; ++i) {
            if (IS_MSVC && i == 0)
                continue;
            auto const tim = Time();
            for (auto dsr: divisors)
                x.self_div(dsr, [](auto n){ return BarrettRS64(n); },
                    [i](auto dnd_hi, auto dnd_lo, auto dsr, auto const & prep){
                        return i == 0 ? BarrettDivMod64<true>(dnd_hi, dnd_lo, dsr, prep.first, prep.second) :
                            BarrettDivMod64<false>(dnd_hi, dnd_lo, dsr, prep.first, prep.second);
                    });
            double const bar_time = Time() - tim;
            std::cout << "Barrett" << (i == 0 ? "128" : "64 ") << " time " << bar_time << " sec, boost " << div_time / bar_time << std::endl;
        }
    }
}

int main() {
    try {
        TestBarrett();
        TestSpeed();
        return 0;
    } catch (std::exception const & ex) {
        std::cout << "Exception: " << ex.what() << std::endl;
        return -1;
    }
}