add bulk_invs and pow_fixed helpers

Author: adamant-pwn

cp-algo/math/combinatorics.hpp

@@ -33,17 +33,41 @@ namespace cp_algo::math {
         }
         return F[n];
     }
-    template<typename T>
-    T small_inv(auto n) {
-        static std::vector<T> F(maxn);
+    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) {
-            for(int i = 1; i < maxn; i++) {
-                F[i] = rfact<T>(i) * fact<T>(i - 1);
-            }
             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 F[n];
+        return prec_low[n & 0xFFFF] * prec_high[n >> 16];
+    }
+    template<typename T>
+    std::vector<T> bulk_invs(auto const& args) {
+        std::vector<T> res(size(args), args[0]);
+        for(size_t i = 1; i < size(args); i++) {
+            res[i] = res[i - 1] * args[i];
+        }
+        auto all_invs = T(1) / res.back();
+        for(size_t i = 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 = builk_invs<T>(std::views::iota(1) | std::views::take(maxn));
+        return F[n - 1];
     }
     template<typename T>
     T binom_large(T n, auto r) {

cp-algo/math/factorials.hpp

@@ -3,7 +3,7 @@
 #include "cp-algo/util/checkpoint.hpp"
 #include "cp-algo/util/bump_alloc.hpp"
 #include "cp-algo/util/simd.hpp"
-#include "cp-algo/math/common.hpp"
+#include "cp-algo/math/combinatorics.hpp"
 #include "cp-algo/number_theory/modint.hpp"
 #include <ranges>
 
@@ -37,38 +37,37 @@ namespace cp_algo::math {
                 t = mod - t - 1;
                 y = t % 2 ? 1 : mod-1;
             }
-            int pw = 0;
+            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 *= bpow(base(2), pw);
+            y *= pow_fixed<base, 2>(int(pw % (mod - 1)));
         }
         checkpoint("init");
-        uint32_t b2x32 = (1ULL << 32) % mod;
+        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];
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
-                        cur[z][j] = uint32_t(uint64_t(cur[z][j]) * b2x32 % mod);
+                        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;
-                        cur[z] = cur[z] >= mod ? cur[z] - mod : cur[z];
                         prods[z][i] = montgomery_mul(prods[z][i - 1], cur[z], mod, imod);
                     }
                 }
@@ -85,12 +84,12 @@ namespace cp_algo::math {
                 checkpoint("mul ans");
             }
         };
-        uint32_t b2x33 = (1ULL << 33) % mod;
-        process(limit_reg, reg_args_per_block, b2x32, std::identity{});
-        process(limit_odd, odd_args_per_block, b2x33, [](uint32_t x) {return 2 * x + 1;});
+        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 = x.inv();
+                x = invs[i];
             }
         }
         checkpoint("inv ans");