better reading for bitpack

Author: adamant-pwn

cp-algo/structures/bit_array.hpp

@@ -14,6 +14,9 @@ namespace cp_algo::structures {
         void set(size_t x) {
             data[x / width] |= 1ULL << (x % width);
         }
+        void reset(size_t x) {
+            data[x / width] &= ~(1ULL << (x % width));
+        }
         void flip(size_t x) {
             data[x / width] ^= 1ULL << (x % width);
         }

cp-algo/structures/bitpack.hpp

@@ -1,28 +1,24 @@
 #ifndef CP_ALGO_STRUCTURES_BITPACK_HPP
 #define CP_ALGO_STRUCTURES_BITPACK_HPP
 #include "../structures/bit_array.hpp"
+#include "../util/simd.hpp"
 #include <cstdint>
 #include <cstddef>
 #include <string>
 #include <array>
 namespace cp_algo::structures {
-    template<size_t n, typename Int = uint64_t>
-    struct bitpack: bit_array<n, Int> {
-        using Base = bit_array<n, Int>;
+    template<size_t n>
+    struct bitpack: bit_array<n, uint64_t> {
+        using Base = bit_array<n, uint64_t>;
         using Base::width, Base::blocks, Base::data;
+        using Base::set, Base::reset;
         auto operator <=> (bitpack const& t) const = default;
 
         bitpack() {}
-        bitpack(std::string bits) {
-            size_t rem = size(bits) % width;
-            if(rem) {
-                bits += std::string(width - rem, '0');
-            }
-            for(size_t i = 0, pos = 0; pos < size(bits); i++, pos += width) {
-                for(size_t j = width; j; j--) {
-                    data[i] *= 2;
-                    data[i] ^= bits[pos + j - 1] == '1';
-                }
+        bitpack(std::string &bits) {
+            bits.resize((size(bits) + width - 1) / width * width);
+            for(size_t i = 0; i < blocks; i++) {
+                data[i] = read_bits64(bits.data() + i * width);
             }
         }
 
@@ -42,7 +38,7 @@ namespace cp_algo::structures {
         std::string to_string() const {
             std::string res(blocks * width, '0');
             for(size_t i = 0, pos = 0; i < blocks; i++, pos += width) {
-                Int block = data[i];
+                auto block = data[i];
                 for(size_t j = 0; j < width; j++) {
                     res[pos + j] = '0' + block % 2;
                     block /= 2;

cp-algo/util/bit.hpp

@@ -1,6 +1,6 @@
 #ifndef CP_ALGO_UTIL_BIT_HPP
 #define CP_ALGO_UTIL_BIT_HPP
-#include <immintrin.h>
+#include "../util/simd.hpp"
 #include <cstdint>
 #include <array>
 #include <bit>
@@ -25,5 +25,12 @@ namespace cp_algo {
             callback.template operator()<1ULL << fl>();
         }
     }
+
+    [[gnu::target("avx2"), gnu::always_inline]] inline uint32_t read_bits(char const* p) {
+        return _mm256_movemask_epi8(__m256i(vector_cast<u8x32 const>(p[0]) + (127 - '0')));
+    }
+    [[gnu::always_inline]] inline uint64_t read_bits64(char const* p) {
+        return read_bits(p) | (uint64_t(read_bits(p + 32)) << 32);
+    }
 }
 #endif // CP_ALGO_UTIL_BIT_HPP

cp-algo/util/simd.hpp

@@ -7,12 +7,13 @@
 namespace cp_algo {
     template<typename T, size_t len>
     using simd [[gnu::vector_size(len * sizeof(T))]] = T;
-    using u32x8 = simd<uint32_t, 8>;
     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::always_inline]] inline dx4 abs(dx4 a) {
@@ -44,23 +45,14 @@ namespace cp_algo {
     [[gnu::always_inline]] inline auto swap_bytes(auto x) {
         return decltype(x)(__builtin_shufflevector(u32x8(x), u32x8(x), 1, 0, 3, 2, 5, 4, 7, 6));
     }
-    [[gnu::always_inline]] inline u64x4 montgomery_reduce(u64x4 x, uint32_t mod, uint32_t imod) {
-#ifdef __AVX2__
+    [[gnu::target("avx2"), gnu::always_inline]] 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));
-#else
-        auto x_ninv = x * imod;
-        x += low32(x_ninv) * mod;
-#endif
         return swap_bytes(x);
     }
 
-    [[gnu::always_inline]] inline u64x4 montgomery_mul(u64x4 x, u64x4 y, uint32_t mod, uint32_t imod) {
-#ifdef __AVX2__
+    [[gnu::target("avx2"), gnu::always_inline]] 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);
-#else
-        return montgomery_reduce(low32(x) * low32(y), mod, imod);
-#endif
     }
     [[gnu::always_inline]] inline u32x8 montgomery_mul(u32x8 x, u32x8 y, uint32_t mod, uint32_t imod) {
         return u32x8(montgomery_mul(u64x4(x), u64x4(y), mod, imod)) |

verify/structures/bitpack/prod_mod_2.test.cpp

@@ -40,11 +40,12 @@ void solve() {
         b[i] = row;
     }
     cp_algo::checkpoint("read");
-    for(int j = 0; j < m; j += 64) {
-        for(int z = 0; z < 64 / K; z++) {
+    const int width = bitpack<maxn>::width;
+    for(int j = 0; j < m; j += width) {
+        for(int z = 0; z < width / K; z++) {
             process_precalc(j / K + z);
             for(int i = 0; i < n; i++) {
-                c[i] ^= precalc[uint8_t(a[i].word(j / 64) >> K * z)];
+                c[i] ^= precalc[uint8_t(a[i].word(j / width) >> K * z)];
             }
         }
     }