CP-Algorithms Library

This documentation is automatically generated by competitive-verifier/competitive-verifier

View the Project on GitHub cp-algorithms/cp-algorithms-aux

:heavy_check_mark: cp-algo/math/factorials.hpp

Depends on

Verified with

Code

#ifndef CP_ALGO_MATH_FACTORIALS_HPP
#define CP_ALGO_MATH_FACTORIALS_HPP
#include "cp-algo/util/checkpoint.hpp"
#include "cp-algo/util/bump_alloc.hpp"
#include "cp-algo/util/simd.hpp"
#include "cp-algo/math/combinatorics.hpp"
#include "cp-algo/number_theory/modint.hpp"
#include <ranges>

namespace cp_algo::math {
    template<bool use_bump_alloc = false, int maxn = -1>
    auto facts(auto const& args) {
        static_assert(!use_bump_alloc || maxn > 0, "maxn must be set if use_bump_alloc is true");
        constexpr int max_mod = 1'000'000'000;
        constexpr int accum = 4;
        constexpr int simd_size = 8;
        constexpr int block = 1 << 18;
        constexpr int subblock = block / simd_size;
        using base = std::decay_t<decltype(args[0])>;
        static_assert(modint_type<base>, "Base type must be a modint type");
        using T = std::array<int, 2>;
        using alloc = std::conditional_t<use_bump_alloc,
            bump_alloc<T, 30 * maxn>,
            big_alloc<T>>;
        std::basic_string<T, std::char_traits<T>, alloc> odd_args_per_block[max_mod / subblock];
        std::basic_string<T, std::char_traits<T>, alloc> reg_args_per_block[max_mod / subblock];
        constexpr int limit_reg = max_mod / 64;
        int limit_odd = 0;

        std::vector<base, big_alloc<base>> res(size(args), 1);
        const int mod = base::mod();
        const int imod = -math::inv2(mod);
        for(auto [i, xy]: std::views::zip(args, res) | std::views::enumerate) {
            auto [x, y] = xy;
            int t = x.getr();
            if(t >= mod / 2) {
                t = mod - t - 1;
                y = t % 2 ? 1 : mod-1;
            }
            auto pw = 32ull * (t + 1);
            while(t > limit_reg) {
                limit_odd = std::max(limit_odd, (t - 1) / 2);
                odd_args_per_block[(t - 1) / 2 / subblock].push_back({int(i), (t - 1) / 2});
                t /= 2;
                pw += t;
            }
            reg_args_per_block[t / subblock].push_back({int(i), t});
            y *= pow_fixed<base, 2>(int(pw % (mod - 1)));
        }
        checkpoint("init");
        base bi2x32 = pow_fixed<base, 2>(32).inv();
        auto process = [&](int limit, auto &args_per_block, auto step, auto &&proj) {
            base fact = 1;
            for(int b = 0; b <= limit; b += accum * block) {
                u32x8 cur[accum];
                static std::array<u32x8, subblock> prods[accum];
                for(int z = 0; z < accum; z++) {
                    for(int j = 0; j < simd_size; j++) {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
                        cur[z][j] = uint32_t(b + z * block + j * subblock);
                        cur[z][j] = proj(cur[z][j]);
                        prods[z][0][j] = cur[z][j] + !cur[z][j];
                        prods[z][0][j] = uint32_t(uint64_t(prods[z][0][j]) * bi2x32.getr() % mod);
#pragma GCC diagnostic pop
                    }
                }
                for(int i = 1; i < block / simd_size; i++) {
                    for(int z = 0; z < accum; z++) {
                        cur[z] += step;
                        prods[z][i] = montgomery_mul(prods[z][i - 1], cur[z], mod, imod);
                    }
                }
                checkpoint("inner loop");
                for(int z = 0; z < accum; z++) {
                    for(int j = 0; j < simd_size; j++) {
                        int bl = b + z * block + j * subblock;
                        for(auto [i, x]: args_per_block[bl / subblock]) {
                            res[i] *= fact * prods[z][x - bl][j];
                        }
                        fact *= base(prods[z].back()[j]);
                    }
                }
                checkpoint("mul ans");
            }
        };
        process(limit_reg, reg_args_per_block, 1, std::identity{});
        process(limit_odd, odd_args_per_block, 2, [](uint32_t x) {return 2 * x + 1;});
        auto invs = bulk_invs<base>(res);
        for(auto [i, x]: res | std::views::enumerate) {
            if (args[i] >= mod / 2) {
                x = invs[i];
            }
        }
        checkpoint("inv ans");
        return res;
    }
}
#endif // CP_ALGO_MATH_FACTORIALS_HPP

#line 1 "cp-algo/math/factorials.hpp"


#line 1 "cp-algo/util/checkpoint.hpp"


#include <iostream>
#include <chrono>
#include <string>
#include <map>
namespace cp_algo {
    std::map<std::string, double> checkpoints;
    template<bool final = false>
    void checkpoint([[maybe_unused]] std::string const& msg = "") {
#ifdef CP_ALGO_CHECKPOINT
        static double last = 0;
        double now = (double)clock() / CLOCKS_PER_SEC;
        double delta = now - last;
        last = now;
        if(msg.size() && !final) {
            checkpoints[msg] += delta;
        }
        if(final) {
            for(auto const& [key, value] : checkpoints) {
                std::cerr << key << ": " << value * 1000 << " ms\n";
            }
            std::cerr << "Total: " << now * 1000 << " ms\n";
        }
#endif
    }
}

#line 1 "cp-algo/util/bump_alloc.hpp"


#include <cstddef>
#line 1 "cp-algo/util/big_alloc.hpp"



#line 6 "cp-algo/util/big_alloc.hpp"

// Single macro to detect POSIX platforms (Linux, Unix, macOS)
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
#  define CP_ALGO_USE_MMAP 1
#  include <sys/mman.h>
#else
#  define CP_ALGO_USE_MMAP 0
#endif

namespace cp_algo {
    template <typename T, std::size_t Align = 32>
    class big_alloc {
        static_assert( Align >= alignof(void*), "Align must be at least pointer-size");
        static_assert(std::popcount(Align) == 1, "Align must be a power of two");
    public:
        using value_type = T;
        template <class U> struct rebind { using other = big_alloc<U, Align>; };
        constexpr bool operator==(const big_alloc&) const = default;
        constexpr bool operator!=(const big_alloc&) const = default;

        big_alloc() noexcept = default;
        template <typename U, std::size_t A>
        big_alloc(const big_alloc<U, A>&) noexcept {}

        [[nodiscard]] T* allocate(std::size_t n) {
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align = std::max<std::size_t>(alignof(T),  Align);
#if CP_ALGO_USE_MMAP
            if (padded >= MEGABYTE) {
                void* raw = mmap(nullptr, padded,
                                PROT_READ | PROT_WRITE,
                                MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
                madvise(raw, padded, MADV_HUGEPAGE);
                madvise(raw, padded, MADV_POPULATE_WRITE);
                return static_cast<T*>(raw);
            }
#endif
            return static_cast<T*>(::operator new(padded, std::align_val_t(align)));
        }

        void deallocate(T* p, std::size_t n) noexcept {
            if (!p) return;
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align  = std::max<std::size_t>(alignof(T),  Align);
    #if CP_ALGO_USE_MMAP
            if (padded >= MEGABYTE) { munmap(p, padded); return; }
    #endif
            ::operator delete(p, padded, std::align_val_t(align));
        }

    private:
        static constexpr std::size_t MEGABYTE = 1 << 20;
        static constexpr std::size_t round_up(std::size_t x) noexcept {
            return (x + Align - 1) / Align * Align;
        }
    };
}

#line 5 "cp-algo/util/bump_alloc.hpp"
namespace cp_algo {
    template<class T, size_t max_len>
    struct bump_alloc {
        static char* buf;
        static size_t buf_ind;
        using value_type = T;
        template <class U> struct rebind { using other = bump_alloc<U, max_len>; };
        constexpr bool operator==(const bump_alloc&) const = default;
        constexpr bool operator!=(const bump_alloc&) const = default;
        bump_alloc() = default;
        template<class U> bump_alloc(const U&) {}
        T* allocate(size_t n) {
            buf_ind -= n * sizeof(T);
            buf_ind &= 0 - alignof(T);
            return (T*)(buf + buf_ind);
        }
        void deallocate(T*, size_t) {}
    };
    template<class T, size_t max_len>
    char* bump_alloc<T, max_len>::buf = big_alloc<char>().allocate(max_len * sizeof(T));
    template<class T, size_t max_len>
    size_t bump_alloc<T, max_len>::buf_ind = max_len * sizeof(T);
}

#line 1 "cp-algo/util/simd.hpp"


#include <experimental/simd>
#include <cstdint>
#line 6 "cp-algo/util/simd.hpp"
#include <memory>
namespace cp_algo {
    template<typename T, size_t len>
    using simd [[gnu::vector_size(len * sizeof(T))]] = T;
    using i64x4 = simd<int64_t, 4>;
    using u64x4 = simd<uint64_t, 4>;
    using u32x8 = simd<uint32_t, 8>;
    using i32x4 = simd<int32_t, 4>;
    using u32x4 = simd<uint32_t, 4>;
    using i16x4 = simd<int16_t, 4>;
    using u8x32 = simd<uint8_t, 32>;
    using dx4 = simd<double, 4>;

    [[gnu::target("avx2")]] inline dx4 abs(dx4 a) {
    return a < 0 ? -a : a;
    }

    // https://stackoverflow.com/a/77376595
    // works for ints in (-2^51, 2^51)
    static constexpr dx4 magic = dx4() + (3ULL << 51);
    [[gnu::target("avx2")]] inline i64x4 lround(dx4 x) {
        return i64x4(x + magic) - i64x4(magic);
    }
    [[gnu::target("avx2")]] inline dx4 to_double(i64x4 x) {
        return dx4(x + i64x4(magic)) - magic;
    }

    [[gnu::target("avx2")]] inline dx4 round(dx4 a) {
        return dx4{
            std::nearbyint(a[0]),
            std::nearbyint(a[1]),
            std::nearbyint(a[2]),
            std::nearbyint(a[3])
        };
    }

    [[gnu::target("avx2")]] inline u64x4 low32(u64x4 x) {
        return x & uint32_t(-1);
    }
    [[gnu::target("avx2")]] inline auto swap_bytes(auto x) {
        return decltype(x)(__builtin_shufflevector(u32x8(x), u32x8(x), 1, 0, 3, 2, 5, 4, 7, 6));
    }
    [[gnu::target("avx2")]] inline u64x4 montgomery_reduce(u64x4 x, uint32_t mod, uint32_t imod) {
        auto x_ninv = u64x4(_mm256_mul_epu32(__m256i(x), __m256i() + imod));
        x += u64x4(_mm256_mul_epu32(__m256i(x_ninv), __m256i() + mod));
        return swap_bytes(x);
    }

    [[gnu::target("avx2")]] inline u64x4 montgomery_mul(u64x4 x, u64x4 y, uint32_t mod, uint32_t imod) {
        return montgomery_reduce(u64x4(_mm256_mul_epu32(__m256i(x), __m256i(y))), mod, imod);
    }
    [[gnu::target("avx2")]] inline u32x8 montgomery_mul(u32x8 x, u32x8 y, uint32_t mod, uint32_t imod) {
        return u32x8(montgomery_mul(u64x4(x), u64x4(y), mod, imod)) |
               u32x8(swap_bytes(montgomery_mul(u64x4(swap_bytes(x)), u64x4(swap_bytes(y)), mod, imod)));
    }
    [[gnu::target("avx2")]] inline dx4 rotate_right(dx4 x) {
        static constexpr u64x4 shuffler = {3, 0, 1, 2};
        return __builtin_shuffle(x, shuffler);
    }

    template<std::size_t Align = 32>
    [[gnu::target("avx2")]] inline bool is_aligned(const auto* p) noexcept {
        return (reinterpret_cast<std::uintptr_t>(p) % Align) == 0;
    }

    template<class Target>
    [[gnu::target("avx2")]] inline Target& vector_cast(auto &&p) {
        return *reinterpret_cast<Target*>(std::assume_aligned<alignof(Target)>(&p));
    }
}

#line 1 "cp-algo/math/combinatorics.hpp"


#line 1 "cp-algo/math/common.hpp"


#include <functional>
#line 5 "cp-algo/math/common.hpp"
#include <cassert>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
    const int maxn = CP_ALGO_MAXN;
#else
    const int maxn = 1 << 19;
#endif
    const int magic = 64; // threshold for sizes to run the naive algo

    auto bpow(auto const& x, auto n, auto const& one, auto op) {
        if(n == 0) {
            return one;
        } else {
            auto t = bpow(x, n / 2, one, op);
            t = op(t, t);
            if(n % 2) {
                t = op(t, x);
            }
            return t;
        }
    }
    auto bpow(auto x, auto n, auto ans) {
        return bpow(x, n, ans, std::multiplies{});
    }
    template<typename T>
    T bpow(T const& x, auto n) {
        return bpow(x, n, T(1));
    }
    inline constexpr auto inv2(auto x) {
        assert(x % 2);
        std::make_unsigned_t<decltype(x)> y = 1;
        while(y * x != 1) {
            y *= 2 - x * y;
        }
        return y;
    }
}

#line 5 "cp-algo/math/combinatorics.hpp"
#include <ranges>
namespace cp_algo::math {
    // fact/rfact/small_inv are caching
    // Beware of usage with dynamic mod
    template<typename T>
    T fact(auto n) {
        static std::vector<T> F(maxn);
        static bool init = false;
        if(!init) {
            F[0] = T(1);
            for(int i = 1; i < maxn; i++) {
                F[i] = F[i - 1] * T(i);
            }
            init = true;
        }
        return F[n];
    }
    // Only works for modint types
    template<typename T>
    T rfact(auto n) {
        static std::vector<T> F(maxn);
        static bool init = false;
        if(!init) {
            int t = (int)std::min<int64_t>(T::mod(), maxn) - 1;
            F[t] = T(1) / fact<T>(t);
            for(int i = t - 1; i >= 0; i--) {
                F[i] = F[i + 1] * T(i + 1);
            }
            init = true;
        }
        return F[n];
    }
    template<typename T, int base>
    T pow_fixed(int n) {
        static std::vector<T> prec_low(1 << 16);
        static std::vector<T> prec_high(1 << 16);
        static bool init = false;
        if(!init) {
            init = true;
            prec_low[0] = prec_high[0] = T(1);
            T step_low = T(base);
            T step_high = bpow(T(base), 1 << 16);
            for(int i = 1; i < (1 << 16); i++) {
                prec_low[i] = prec_low[i - 1] * step_low;
                prec_high[i] = prec_high[i - 1] * step_high;
            }
        }
        return prec_low[n & 0xFFFF] * prec_high[n >> 16];
    }
    template<typename T>
    std::vector<T> bulk_invs(auto const& args) {
        std::vector<T> res(std::size(args), args[0]);
        for(size_t i = 1; i < std::size(args); i++) {
            res[i] = res[i - 1] * args[i];
        }
        auto all_invs = T(1) / res.back();
        for(size_t i = std::size(args) - 1; i > 0; i--) {
            res[i] = all_invs * res[i - 1];
            all_invs *= args[i];
        }
        res[0] = all_invs;
        return res;
    }
    template<typename T>
    T small_inv(auto n) {
        static auto F = bulk_invs<T>(std::views::iota(1, maxn));
        return F[n - 1];
    }
    template<typename T>
    T binom_large(T n, auto r) {
        assert(r < maxn);
        T ans = 1;
        for(decltype(r) i = 0; i < r; i++) {
            ans = ans * T(n - i) * small_inv<T>(i + 1);
        }
        return ans;
    }
    template<typename T>
    T binom(auto n, auto r) {
        if(r < 0 || r > n) {
            return T(0);
        } else if(n >= maxn) {
            return binom_large(T(n), r);
        } else {
            return fact<T>(n) * rfact<T>(r) * rfact<T>(n - r);
        }
    }
}

#line 1 "cp-algo/number_theory/modint.hpp"


#line 6 "cp-algo/number_theory/modint.hpp"
namespace cp_algo::math {

    template<typename modint, typename _Int>
    struct modint_base {
        using Int = _Int;
        using UInt = std::make_unsigned_t<Int>;
        static constexpr size_t bits = sizeof(Int) * 8;
        using Int2 = std::conditional_t<bits <= 32, int64_t, __int128_t>;
        using UInt2 = std::conditional_t<bits <= 32, uint64_t, __uint128_t>;
        constexpr static Int mod() {
            return modint::mod();
        }
        constexpr static Int remod() {
            return modint::remod();
        }
        constexpr static UInt2 modmod() {
            return UInt2(mod()) * mod();
        }
        constexpr modint_base() = default;
        constexpr modint_base(Int2 rr) {
            to_modint().setr(UInt((rr + modmod()) % mod()));
        }
        modint inv() const {
            return bpow(to_modint(), mod() - 2);
        }
        modint operator - () const {
            modint neg;
            neg.r = std::min(-r, remod() - r);
            return neg;
        }
        modint& operator /= (const modint &t) {
            return to_modint() *= t.inv();
        }
        modint& operator *= (const modint &t) {
            r = UInt(UInt2(r) * t.r % mod());
            return to_modint();
        }
        modint& operator += (const modint &t) {
            r += t.r; r = std::min(r, r - remod());
            return to_modint();
        }
        modint& operator -= (const modint &t) {
            r -= t.r; r = std::min(r, r + remod());
            return to_modint();
        }
        modint operator + (const modint &t) const {return modint(to_modint()) += t;}
        modint operator - (const modint &t) const {return modint(to_modint()) -= t;}
        modint operator * (const modint &t) const {return modint(to_modint()) *= t;}
        modint operator / (const modint &t) const {return modint(to_modint()) /= t;}
        // Why <=> doesn't work?..
        auto operator == (const modint &t) const {return to_modint().getr() == t.getr();}
        auto operator != (const modint &t) const {return to_modint().getr() != t.getr();}
        auto operator <= (const modint &t) const {return to_modint().getr() <= t.getr();}
        auto operator >= (const modint &t) const {return to_modint().getr() >= t.getr();}
        auto operator < (const modint &t) const {return to_modint().getr() < t.getr();}
        auto operator > (const modint &t) const {return to_modint().getr() > t.getr();}
        Int rem() const {
            UInt R = to_modint().getr();
            return R - (R > (UInt)mod() / 2) * mod();
        }
        constexpr void setr(UInt rr) {
            r = rr;
        }
        constexpr UInt getr() const {
            return r;
        }

        // Only use these if you really know what you're doing!
        static UInt modmod8() {return UInt(8 * modmod());}
        void add_unsafe(UInt t) {r += t;}
        void pseudonormalize() {r = std::min(r, r - modmod8());}
        modint const& normalize() {
            if(r >= (UInt)mod()) {
                r %= mod();
            }
            return to_modint();
        }
        void setr_direct(UInt rr) {r = rr;}
        UInt getr_direct() const {return r;}
    protected:
        UInt r;
    private:
        constexpr modint& to_modint() {return static_cast<modint&>(*this);}
        constexpr modint const& to_modint() const {return static_cast<modint const&>(*this);}
    };
    template<typename modint>
    concept modint_type = std::is_base_of_v<modint_base<modint, typename modint::Int>, modint>;
    template<modint_type modint>
    decltype(std::cin)& operator >> (decltype(std::cin) &in, modint &x) {
        typename modint::UInt r;
        auto &res = in >> r;
        x.setr(r);
        return res;
    }
    template<modint_type modint>
    decltype(std::cout)& operator << (decltype(std::cout) &out, modint const& x) {
        return out << x.getr();
    }

    template<auto m>
    struct modint: modint_base<modint<m>, decltype(m)> {
        using Base = modint_base<modint<m>, decltype(m)>;
        using Base::Base;
        static constexpr Base::Int mod() {return m;}
        static constexpr Base::UInt remod() {return m;}
        auto getr() const {return Base::r;}
    };

    template<typename Int = int>
    struct dynamic_modint: modint_base<dynamic_modint<Int>, Int> {
        using Base = modint_base<dynamic_modint<Int>, Int>;
        using Base::Base;

        static Base::UInt m_reduce(Base::UInt2 ab) {
            if(mod() % 2 == 0) [[unlikely]] {
                return typename Base::UInt(ab % mod());
            } else {
                typename Base::UInt2 m = typename Base::UInt(ab) * imod();
                return typename Base::UInt((ab + m * mod()) >> Base::bits);
            }
        }
        static Base::UInt m_transform(Base::UInt a) {
            if(mod() % 2 == 0) [[unlikely]] {
                return a;
            } else {
                return m_reduce(a * pw128());
            }
        }
        dynamic_modint& operator *= (const dynamic_modint &t) {
            Base::r = m_reduce(typename Base::UInt2(Base::r) * t.r);
            return *this;
        }
        void setr(Base::UInt rr) {
            Base::r = m_transform(rr);
        }
        Base::UInt getr() const {
            typename Base::UInt res = m_reduce(Base::r);
            return std::min(res, res - mod());
        }
        static Int mod() {return m;}
        static Int remod() {return 2 * m;}
        static Base::UInt imod() {return im;}
        static Base::UInt2 pw128() {return r2;}
        static void switch_mod(Int nm) {
            m = nm;
            im = m % 2 ? inv2(-m) : 0;
            r2 = static_cast<Base::UInt>(static_cast<Base::UInt2>(-1) % m + 1);
        }

        // Wrapper for temp switching
        auto static with_mod(Int tmp, auto callback) {
            struct scoped {
                Int prev = mod();
                ~scoped() {switch_mod(prev);}
            } _;
            switch_mod(tmp);
            return callback();
        }
    private:
        static thread_local Int m;
        static thread_local Base::UInt im, r2;
    };
    template<typename Int>
    Int thread_local dynamic_modint<Int>::m = 1;
    template<typename Int>
    dynamic_modint<Int>::Base::UInt thread_local dynamic_modint<Int>::im = -1;
    template<typename Int>
    dynamic_modint<Int>::Base::UInt thread_local dynamic_modint<Int>::r2 = 0;
}

#line 9 "cp-algo/math/factorials.hpp"

namespace cp_algo::math {
    template<bool use_bump_alloc = false, int maxn = -1>
    auto facts(auto const& args) {
        static_assert(!use_bump_alloc || maxn > 0, "maxn must be set if use_bump_alloc is true");
        constexpr int max_mod = 1'000'000'000;
        constexpr int accum = 4;
        constexpr int simd_size = 8;
        constexpr int block = 1 << 18;
        constexpr int subblock = block / simd_size;
        using base = std::decay_t<decltype(args[0])>;
        static_assert(modint_type<base>, "Base type must be a modint type");
        using T = std::array<int, 2>;
        using alloc = std::conditional_t<use_bump_alloc,
            bump_alloc<T, 30 * maxn>,
            big_alloc<T>>;
        std::basic_string<T, std::char_traits<T>, alloc> odd_args_per_block[max_mod / subblock];
        std::basic_string<T, std::char_traits<T>, alloc> reg_args_per_block[max_mod / subblock];
        constexpr int limit_reg = max_mod / 64;
        int limit_odd = 0;

        std::vector<base, big_alloc<base>> res(size(args), 1);
        const int mod = base::mod();
        const int imod = -math::inv2(mod);
        for(auto [i, xy]: std::views::zip(args, res) | std::views::enumerate) {
            auto [x, y] = xy;
            int t = x.getr();
            if(t >= mod / 2) {
                t = mod - t - 1;
                y = t % 2 ? 1 : mod-1;
            }
            auto pw = 32ull * (t + 1);
            while(t > limit_reg) {
                limit_odd = std::max(limit_odd, (t - 1) / 2);
                odd_args_per_block[(t - 1) / 2 / subblock].push_back({int(i), (t - 1) / 2});
                t /= 2;
                pw += t;
            }
            reg_args_per_block[t / subblock].push_back({int(i), t});
            y *= pow_fixed<base, 2>(int(pw % (mod - 1)));
        }
        checkpoint("init");
        base bi2x32 = pow_fixed<base, 2>(32).inv();
        auto process = [&](int limit, auto &args_per_block, auto step, auto &&proj) {
            base fact = 1;
            for(int b = 0; b <= limit; b += accum * block) {
                u32x8 cur[accum];
                static std::array<u32x8, subblock> prods[accum];
                for(int z = 0; z < accum; z++) {
                    for(int j = 0; j < simd_size; j++) {
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
                        cur[z][j] = uint32_t(b + z * block + j * subblock);
                        cur[z][j] = proj(cur[z][j]);
                        prods[z][0][j] = cur[z][j] + !cur[z][j];
                        prods[z][0][j] = uint32_t(uint64_t(prods[z][0][j]) * bi2x32.getr() % mod);
#pragma GCC diagnostic pop
                    }
                }
                for(int i = 1; i < block / simd_size; i++) {
                    for(int z = 0; z < accum; z++) {
                        cur[z] += step;
                        prods[z][i] = montgomery_mul(prods[z][i - 1], cur[z], mod, imod);
                    }
                }
                checkpoint("inner loop");
                for(int z = 0; z < accum; z++) {
                    for(int j = 0; j < simd_size; j++) {
                        int bl = b + z * block + j * subblock;
                        for(auto [i, x]: args_per_block[bl / subblock]) {
                            res[i] *= fact * prods[z][x - bl][j];
                        }
                        fact *= base(prods[z].back()[j]);
                    }
                }
                checkpoint("mul ans");
            }
        };
        process(limit_reg, reg_args_per_block, 1, std::identity{});
        process(limit_odd, odd_args_per_block, 2, [](uint32_t x) {return 2 * x + 1;});
        auto invs = bulk_invs<base>(res);
        for(auto [i, x]: res | std::views::enumerate) {
            if (args[i] >= mod / 2) {
                x = invs[i];
            }
        }
        checkpoint("inv ans");
        return res;
    }
}
Back to top page