include/fastmath.hpp - main - Ant010ff/ffpp

include/fastmath.hpp
160 строк6 KB
Ant010ff
Version 2.11.3.
08 авг 2026, 16:44
08 авг 2026, 16:44fd187d3
/*
This header is part of a Functional Flow Processing Primitives (FFPP) library, version 2.11.3.

Official repository: https://gitlab.com/ant010ff/ffpp

Licensed under the MIT License <http://opensource.org/licenses/MIT>.
SPDX-License-Identifier: MIT
Copyright (c) 2021 - 2026 Anton Nasonov <ant010fff @ gmail . com>.
*/

#ifndef FFPP_FAST_MATH_HPP
#define FFPP_FAST_MATH_HPP

namespace FFPP::math {

    FFPP_ATTR_INLINE
    inline size_t MulHi(size_t x, size_t y) noexcept {
        if constexpr(sizeof(size_t) > 4) {
            #if defined(_MSC_VER) && defined(_M_X64)
                return __umulh(x, y);
            #elif (defined(__GNUC__) || defined(__clang__)) && (defined(__x86_64__) || defined(__aarch64__))
                return static_cast<uint64_t>(((__uint128_t)(x) * y) >> 64);
            #else
                uint64_t const
                    x0 = static_cast<uint64_t>(x) & uint64_t(0xFFFFFFFF)
                    , x1 = static_cast<uint64_t>(x) >> 32
                    , y0 = static_cast<uint64_t>(y) & uint64_t(0xFFFFFFFF)
                    , y1 = static_cast<uint64_t>(y) >> 32
                    , w0 = x0 * y0
                    , t = x1 * y0 + (w0 >> 32)
                    , w1 = (t & uint64_t(0xFFFFFFFF)) + (x0 * y1)
                    , w2 = (t >> 32) + (w1 >> 32)
                ;
                return static_cast<size_t>(x1 * y1 + w2);
            #endif
        }
        else {
            return static_cast<size_t>((static_cast<uint64_t>(x) * y) >> 32);
        }
    }

    FFPP_ATTR_INLINE
    inline size_t Sqrt(size_t nSquare) noexcept {
        if(nSquare < 3) return nSquare == 0 ? 0 : 1;
        size_t constexpr 
            c_nHalf = std::numeric_limits<size_t>::digits / 2,
            c_nThreshold = size_t(1) << c_nHalf;
        size_t nLower = 1, nUpper = nSquare, nSqrt = 0;
        if(nSquare < c_nThreshold) while(nLower <= nUpper) {
            size_t const nMid = nLower + ((nUpper - nLower) >> 1);
            if(nMid * nMid <= nSquare) {
                nSqrt = nMid;
                nLower = nMid + 1;
            } 
            else {
                nUpper = nMid - 1;
            }
        }
        else while(nLower <= nUpper) {
            size_t const nMid = nLower + ((nUpper - nLower) >> 1);
            if(MulHi(nMid, nMid) == 0) {
                if(nMid * nMid <= nSquare) {
                    nSqrt = nMid;
                    nLower = nMid + 1;
                } 
                else {
                    nUpper = nMid - 1;
                }
            } 
            else {
                nUpper = nMid - 1;
            }
        }            
        return nSqrt;
    }

    //FastMod: implements fast modulo calculation for a run-time constant divisor based on the Barret reduction algorithm.
    class FastMod {
    public:
        FastMod() { }

        FastMod(size_t nDivisor) 
        : nMultiplier_([nDivisor = nDivisor < 2 ? 1 : nDivisor] () {
            size_t constexpr nHalf = size_t(1) << (std::numeric_limits<size_t>::digits - 1);
            if(nDivisor >= nHalf) Except<std::range_error>::Throw("Divisor (", nDivisor, ") is out of supported range.");
            size_t const 
                q0 = nHalf / nDivisor
                , r0 = nHalf % nDivisor
                , nMultiplier = q0 << 1
            ;
            if(r0 >= ((nDivisor + 1) >> 1)) return nMultiplier + 1;
            else return nMultiplier;
        } ()) 
        , nDivisor_(nDivisor < 2 ? 1 : nDivisor)
        , mskDivisor_(std::has_single_bit(nDivisor) ? (nDivisor - 1) : 0)
        { }

		FastMod(FastMod const& that) = default;		
        FastMod(FastMod&& that) = default;

		FastMod& operator = (FastMod const& that) = default;
        FastMod& operator = (FastMod&& that) = default;

        FFPP_ATTR_INLINE
        size_t operator () (size_t nDividend) const noexcept {
            assert(nDivisor_ != 0);
            if(mskDivisor_ != 0) [[likely]] {
                return nDividend & mskDivisor_;
            }
            else if(nDivisor_ == 1) [[unlikely]] {
                return 0;
            }
            else {
                size_t nMod = nDividend - static_cast<size_t>(MulHi(nMultiplier_, nDividend)) * nDivisor_;
                return (nMod >= nDivisor_) ? (nMod - nDivisor_) : nMod; //expecting conditional move here
            }
        }

        FFPP_ATTR_INLINE
        size_t Divisor() const noexcept { return nDivisor_; }

        FFPP_ATTR_INLINE
        size_t Mask() const noexcept { return mskDivisor_; }

        #if FFPP_ENABLE_DEBUG
            static void Assertions() {
                auto constexpr nMax = ((size_t(1) << (std::numeric_limits<size_t>::digits - 1)) - 1);
                assert(FastMod { 7 } (0) == 0);
                assert(FastMod { 7 } (6) == 6);
                assert(FastMod { 7 } (7) == 0);
                assert(FastMod { 7 } (8) == 1);
                assert(FastMod { 7 } (49) == 0);
                assert(FastMod { 7 } (50) == 1);
                assert(FastMod { 7 } (56) == 0);
                assert(FastMod { 7 } (57) == 1);
                assert(FastMod { nMax } (nMax) == 0);
                assert(FastMod { nMax - 1 } (nMax - 3) == nMax - 3);
            }

            static void RunStressTest(size_t nSeed = 12345, size_t nIterations = 10000000) {
                for(size_t i = 0; i < nIterations; ++i) {
                    size_t const 
                        nDivisor = prng::SplitMix(nSeed) & ((size_t(1) << (std::numeric_limits<size_t>::digits - 1)) - 1)
                        , nDividend = prng::SplitMix(nSeed);
                    FastMod const fmod(nDivisor);
                    assert(fmod(nDividend) == nDividend % fmod.Divisor());
                }
            }
        #endif

    private:
        size_t 
            nMultiplier_ = 0,
            nDivisor_ = 1, 
            mskDivisor_ = 0;
    };//FastMod

}//FFPP::math

#endif//FFPP_FAST_MATH_HPP