cp-library
競技プログラミング用 C++ ライブラリ
AND-OR などの論理畳み込み
convolution/boolean_convolution.hpp · GitHub で見る
共通ページ
#include "lib/convolution/boolean_convolution.hpp"
std::vector<unsigned> a{1, 2, 3, 4}, b{5, 6, 7, 8};
auto c = and_or_convolution(a, b); // {7,6,7,0}
*_or_convolution は \(\displaystyle c_k = \bigvee_{i\circ j=k}(a_i \land b_j)\)、*_and_convolution は \(\displaystyle c_k = \bigwedge_{i\circ j=k}(a_i \lor b_j)\) です。各 bit を独立に計算します。
| 添字 \(\circ\) | OR で集約 | AND で集約 | 計算量 |
|---|---|---|---|
| \(i\mathbin{\&}j\) | and_or_convolution(a,b) |
and_and_convolution(a,b) |
\(O(WN\log N)\) |
| \(i\mathbin{\vert}j\) | or_or_convolution(a,b) |
or_and_convolution(a,b) |
\(O(WN\log N)\) |
| \(i\mathbin{\oplus}j\) | xor_or_convolution(a,b) |
xor_and_convolution(a,b) |
\(O(WN\log N)\) |
| \(\gcd(i,j)\) | gcd_or_convolution(a,b) |
gcd_and_convolution(a,b) |
\(O(WN\log\log N)\) |
| \(\operatorname{lcm}(i,j)\) | lcm_or_convolution(a,b) |
lcm_and_convolution(a,b) |
\(O(WN\log\log N)\) |
| \(i\mathbin{\vert}j\)、\(i\mathbin{\&}j = 0\) | subset_or_convolution(a,b) |
subset_and_convolution(a,b) |
\(O(WN\log^2 N)\) |
関数名の前半は添字の操作、後半は係数をまとめる操作です。and_or_convolution なら、全ての組 i,j について c[i & j] |= a[i] & b[j] とする結果になります。AND で集約する版は、各出力先を全 bit が 1 の値で初期化してから &= a[i] | b[j] します。該当する組がない添字には、その初期値が残ります。
係数は bool または 64 bit 以下の整数です。GCD・LCM では添字 0 を使用しません。出力長は GCD・LCM が長い方の長さ、それ以外はそれを 2 の冪に切り上げた長さ \(N\) です。範囲外の LCM は出力しません。\(N \le 2^{30}\) とします。
係数の bit 数を \(W\) として、時間は AND・OR・XOR が \(O(WN \log N)\)、GCD・LCM が \(O(WN \log \log N)\)、Subset が \(O(WN \log^{2} N)\) です。追加空間は Subset が \(O(N \log N)\)、その他が \(O(N)\) です。
依存ライブラリ
- ビット演算畳み込みの modint 用共通実装
- Bitwise Zeta / Möbius Transform
- Divisor Zeta / Möbius Transform
- modint998244353 の GCD 畳み込み
- modint998244353 の共通剰余演算
- Subset Convolution Modular Core
- Walsh–Hadamard Transform
- AND・OR・XOR 畳み込み
- GCD-plus・LCM-plus 畳み込み
- Subset 畳み込み
- modint998244353
- 素数篩
コード
#ifndef BOOLEAN_CONVOLUTION_HPP
#define BOOLEAN_CONVOLUTION_HPP
#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <limits>
#include <type_traits>
#include <vector>
#include "bitwise_convolution.hpp"
#include "gcd_lcm_convolution.hpp"
#include "subset_convolution.hpp"
namespace boolean_convolution_detail {
template <int index, bool dual = false, class T>
std::vector<T> multiply(const std::vector<T> &a, const std::vector<T> &b) {
static_assert(std::is_integral_v<T> && sizeof(T) <= 8 && 0 <= index && index <= 5);
if (a.empty() || b.empty()) return {};
using U = std::make_unsigned_t<std::conditional_t<std::is_same_v<T, bool>, unsigned char, T>>;
std::size_t n = std::max(a.size(), b.size());
if constexpr (index < 3 || index == 5) n = std::bit_ceil(n);
assert(n <= std::size_t(1) << 30);
constexpr U full = std::is_same_v<T, bool> ? U(1) : U(~U(0));
U mask = 0, active = 0;
for (std::size_t i = index == 3 || index == 4; i < a.size() && mask != full; ++i) {
U value = U(a[i]);
mask |= U(dual ? ~value : value) & full;
}
for (std::size_t i = index == 3 || index == 4; i < b.size() && active != mask; ++i) {
U value = U(b[i]);
active |= U(dual ? ~value : value) & mask;
}
std::vector<T> result(n, dual ? T(~U(0)) : T(0));
for (; active; active &= active - 1) {
int bit = std::countr_zero(active);
std::vector<__int128_t> x(n), y(n);
for (std::size_t i = index == 3 || index == 4; i < a.size(); ++i) x[i] = ((U(a[i]) >> bit) & 1) ^ dual;
for (std::size_t i = index == 3 || index == 4; i < b.size(); ++i) y[i] = ((U(b[i]) >> bit) & 1) ^ dual;
std::vector<__int128_t> counts;
if constexpr (index < 3) counts = bitwise_convolution_detail::multiply<index>(std::move(x), std::move(y));
else if constexpr (index == 5) counts = subset_convolution(std::move(x), std::move(y));
else counts = gcd_lcm_convolution_detail::multiply<index == 3>(std::move(x), std::move(y));
for (std::size_t i = 0; i < n; ++i)
if (counts[i]) {
U flag = U(1) << bit;
result[i] = T(dual ? U(result[i]) & ~flag : U(result[i]) | flag);
}
}
if constexpr (index == 3 || index == 4) result[0] = T(0);
return result;
}
}
template <class T>
std::vector<T> and_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<0>(a, b); }
template <class T>
std::vector<T> or_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<1>(a, b); }
template <class T>
std::vector<T> xor_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<2>(a, b); }
template <class T>
std::vector<T> gcd_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<3>(a, b); }
template <class T>
std::vector<T> lcm_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<4>(a, b); }
template <class T>
std::vector<T> and_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<0, true>(a, b); }
template <class T>
std::vector<T> or_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<1, true>(a, b); }
template <class T>
std::vector<T> xor_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<2, true>(a, b); }
template <class T>
std::vector<T> gcd_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<3, true>(a, b); }
template <class T>
std::vector<T> lcm_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<4, true>(a, b); }
template <class T>
std::vector<T> subset_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<5>(a, b); }
template <class T>
std::vector<T> subset_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<5, true>(a, b); }
#endif
#line 1 "convolution/boolean_convolution.hpp"
#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <limits>
#include <type_traits>
#include <vector>
#line 1 "convolution/bitwise_convolution.hpp"
#line 1 "components/bitwise_zeta.hpp"
#line 8 "components/bitwise_zeta.hpp"
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif
namespace bitwise_zeta_detail {
template <bool inverse, bool superset, class T>
void stages(std::vector<T> &a, std::size_t start, std::size_t n, std::size_t half) {
auto step = [&](std::size_t low, std::size_t high) {
T &to = a[superset ? low : high];
const T &from = a[superset ? high : low];
if constexpr (inverse) to -= from;
else to += from;
};
for (; half < n / 2; half *= 4)
for (std::size_t block = start; block < start + n; block += 4 * half)
for (std::size_t j = block; j < block + half; ++j) {
step(j, j + half);
step(j + 2 * half, j + 3 * half);
step(j, j + 2 * half);
step(j + half, j + 3 * half);
}
if (half < n)
for (std::size_t j = start; j < start + half; ++j) step(j, j + half);
}
}
template <bool inverse = false, bool superset = false, class T>
void bitwise_zeta(std::vector<T> &a) {
assert(a.empty() || std::has_single_bit(a.size()));
constexpr std::size_t block = 1024;
if (a.size() <= block) bitwise_zeta_detail::stages<inverse, superset>(a, 0, a.size(), 1);
else {
for (std::size_t start = 0; start < a.size(); start += block)
bitwise_zeta_detail::stages<inverse, superset>(a, start, block, 1);
bitwise_zeta_detail::stages<inverse, superset>(a, 0, a.size(), block);
}
}
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif
#line 1 "components/walsh_hadamard.hpp"
#line 8 "components/walsh_hadamard.hpp"
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif
namespace walsh_hadamard_detail {
template <class T>
void stages(std::vector<T> &a, std::size_t start, std::size_t n, std::size_t half) {
auto step = [&](std::size_t low, std::size_t high) {
T x = a[low], y = a[high];
a[low] = x + y;
a[high] = x - y;
};
for (; half < n / 2; half *= 4)
for (std::size_t block = start; block < start + n; block += 4 * half)
for (std::size_t j = block; j < block + half; ++j) {
step(j, j + half);
step(j + 2 * half, j + 3 * half);
step(j, j + 2 * half);
step(j + half, j + 3 * half);
}
if (half < n)
for (std::size_t j = start; j < start + half; ++j) step(j, j + half);
}
}
template <class T>
void walsh_hadamard(std::vector<T> &a) {
assert(a.empty() || std::has_single_bit(a.size()));
constexpr std::size_t block = 1024;
if (a.size() <= block) walsh_hadamard_detail::stages(a, 0, a.size(), 1);
else {
for (std::size_t start = 0; start < a.size(); start += block)
walsh_hadamard_detail::stages(a, start, block, 1);
walsh_hadamard_detail::stages(a, 0, a.size(), block);
}
}
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif
#line 1 "components/bitwise_convolution_modular.hpp"
#line 1 "components/mint_arithmetic.hpp"
#line 1 "math/mint.hpp"
#line 4 "math/mint.hpp"
#include <istream>
#include <ostream>
#line 7 "math/mint.hpp"
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif
namespace mint_detail {
template <class I>
constexpr bool mint_integral = std::is_integral_v<I>
#if defined(__SIZEOF_INT128__)
|| std::is_same_v<I, __int128> || std::is_same_v<I, unsigned __int128>
#endif
;
struct [[gnu::warn_unused]] Mint {
using mint = Mint;
unsigned v = 0;
[[nodiscard]] static constexpr unsigned mod() { return 998244353; }
[[nodiscard, gnu::always_inline]] static unsigned reduce(unsigned long long x) {
return (x + (unsigned long long)(unsigned(x) * 998244351u) * mod()) >> 32;
}
[[gnu::always_inline]] Mint() = default;
template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
Mint(I x) {
long long r = x % (long long)mod();
if (r < 0) r += mod();
v = reduce((unsigned long long)r * 932051910);
}
[[nodiscard, gnu::always_inline]] static Mint raw(unsigned x) {
Mint a;
a.v = reduce((unsigned long long)x * 932051910);
return a;
}
[[nodiscard, gnu::always_inline]] static Mint word(unsigned x) {
Mint a;
a.v = x;
return a;
}
[[nodiscard, gnu::always_inline]] unsigned val() const {
unsigned r = reduce(v);
return r < mod() ? r : r - mod();
}
friend std::istream &operator>>(std::istream &is, Mint &x) {
long long value;
if (is >> value) x = Mint(value);
return is;
}
friend std::ostream &operator<<(std::ostream &os, const Mint &x) {
return os << x.val();
}
[[gnu::always_inline]] Mint &operator+=(Mint b) {
v += b.v;
if (v >= 2 * mod()) v -= 2 * mod();
return *this;
}
[[gnu::always_inline]] Mint &operator-=(Mint b) {
v -= b.v;
if (v >> 31) v += 2 * mod();
return *this;
}
[[gnu::always_inline]] Mint &operator*=(Mint b) {
v = reduce((unsigned long long)v * b.v);
return *this;
}
Mint &operator/=(Mint b) { return *this *= b.inv(); }
Mint &operator++() { return *this += Mint(1); }
Mint &operator--() { return *this -= Mint(1); }
Mint operator++(int) {
Mint a = *this;
++*this;
return a;
}
Mint operator--(int) {
Mint a = *this;
--*this;
return a;
}
[[nodiscard, gnu::always_inline]] Mint operator+() const { return *this; }
[[nodiscard, gnu::always_inline]] Mint operator-() const { return word(v ? 2 * mod() - v : 0); }
[[nodiscard, gnu::always_inline]] Mint operator+(Mint b) const { return Mint(*this) += b; }
[[nodiscard, gnu::always_inline]] Mint operator-(Mint b) const { return Mint(*this) -= b; }
[[nodiscard, gnu::always_inline]] Mint operator*(Mint b) const { return Mint(*this) *= b; }
[[nodiscard]] Mint operator/(Mint b) const { return Mint(*this) /= b; }
[[nodiscard, gnu::always_inline]] bool operator==(Mint b) const { return (v < mod() ? v : v - mod()) == (b.v < mod() ? b.v : b.v - mod()); }
[[nodiscard, gnu::always_inline]] bool operator!=(Mint b) const { return !(*this == b); }
template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
[[nodiscard]] friend Mint operator+(I a, Mint b) { return Mint(a) + b; }
template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
[[nodiscard]] friend Mint operator-(I a, Mint b) { return Mint(a) - b; }
template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
[[nodiscard]] friend Mint operator*(I a, Mint b) { return Mint(a) * b; }
template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
[[nodiscard]] friend Mint operator/(I a, Mint b) { return Mint(a) / b; }
template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
[[nodiscard]] friend bool operator==(I a, Mint b) { return Mint(a) == b; }
template <class I, std::enable_if_t<mint_integral<I>, int> = 0>
[[nodiscard]] friend bool operator!=(I a, Mint b) { return Mint(a) != b; }
[[nodiscard]] Mint pow(long long n) const {
assert(n >= 0);
if (!n) return Mint(1);
Mint r(1), a = *this;
for (; n > 1; n >>= 1, a *= a)
if (n & 1) r *= a;
return r * a;
}
[[nodiscard, gnu::noinline]] Mint inv() const {
assert(*this != Mint(0));
Mint a2 = *this * *this, a3 = a2 * *this;
Mint a7 = a3 * a3 * *this;
Mint a31 = a7;
for (int i = 0; i < 2; ++i) a31 *= a31;
a31 *= a3;
Mint r = a7;
for (int i = 0; i < 3; ++i) r *= r;
r *= a3;
r *= r;
for (int j = 0; j < 4; ++j) {
for (int i = 0; i < 5; ++i) r *= r;
r *= a31;
}
for (int i = 0; i < 3; ++i) r *= r;
return r * a7;
}
};
} // namespace mint_detail
using mint = mint_detail::Mint;
using fps_mint = mint;
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif
#line 5 "components/mint_arithmetic.hpp"
#if defined(__x86_64__) || defined(__i386__)
#include <immintrin.h>
#endif
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif
namespace mint_detail {
struct Ops {
using T = mint;
typedef unsigned int U __attribute__((may_alias));
static constexpr U mod = 998244353, twice = 2 * mod, one = 301989884;
static U reduce(unsigned long long x) { return T::reduce(x); }
static U mul(U a, U b) { return reduce((unsigned long long)a * b); }
static U add(U a, U b) {
U c = a + b;
return c < twice ? c : c - twice;
}
static U sub(U a, U b) {
U c = a - b;
return c < twice ? c : c + twice;
}
static U in(U a) { return mul(a, 932051910); }
static U out(U a) {
a = reduce(a);
return a < mod ? a : a - mod;
}
static U quotient(U r) { return ((unsigned long long)r << 32) / mod; }
static U fixed_mul(U a, U r, U q) {
return a * r - U(((unsigned long long)a * q) >> 32) * mod;
}
#if defined(__x86_64__) || defined(__i386__)
static bool vectorized() {
#if defined(__AVX2__)
return true;
#else
static const bool available = __builtin_cpu_supports("avx2");
return available;
#endif
}
using I8 = __m256i;
__attribute__((target("avx2"))) static I8 load8(const U *a) {
return _mm256_loadu_si256(reinterpret_cast<const I8 *>(a));
}
__attribute__((target("avx2"))) static void store8(U *a, I8 x) {
_mm256_storeu_si256(reinterpret_cast<I8 *>(a), x);
}
__attribute__((target("avx2"), always_inline)) static inline I8 reduce8(I8 x, I8 y) {
I8 r = _mm256_set1_epi32(998244351), p = _mm256_set1_epi32(mod);
x = _mm256_add_epi64(x, _mm256_mul_epu32(_mm256_mul_epu32(x, r), p));
y = _mm256_add_epi64(y, _mm256_mul_epu32(_mm256_mul_epu32(y, r), p));
return _mm256_blend_epi32(_mm256_srli_epi64(x, 32), y, 0xaa);
}
__attribute__((target("avx2"))) static I8 mul8(I8 a, I8 b) {
return reduce8(_mm256_mul_epu32(a, b), _mm256_mul_epu32(_mm256_srli_epi64(a, 32), _mm256_srli_epi64(b, 32)));
}
__attribute__((target("avx2"))) static I8 fixed_mul8(I8 a, I8 r, I8 q) {
I8 x = _mm256_srli_epi64(_mm256_mul_epu32(a, q), 32);
I8 y = _mm256_mul_epu32(_mm256_srli_epi64(a, 32), q);
I8 t = _mm256_blend_epi32(x, y, 0xaa);
return _mm256_sub_epi32(_mm256_mullo_epi32(a, r),
_mm256_mullo_epi32(t, _mm256_set1_epi32(mod)));
}
__attribute__((target("avx2"))) static I8 fixed_even8(I8 a, I8 r, I8 q) {
I8 t = _mm256_srli_epi64(_mm256_mul_epu32(a, q), 32);
return _mm256_sub_epi64(_mm256_mul_epu32(a, r),
_mm256_mul_epu32(t, _mm256_set1_epi32(mod)));
}
__attribute__((target("avx2"))) static I8 add8(I8 a, I8 b) {
I8 x = _mm256_add_epi32(a, b);
return _mm256_min_epu32(x, _mm256_sub_epi32(x, _mm256_set1_epi32(twice)));
}
__attribute__((target("avx2"))) static I8 sub8(I8 a, I8 b) {
I8 x = _mm256_sub_epi32(a, b);
return _mm256_min_epu32(x, _mm256_add_epi32(x, _mm256_set1_epi32(twice)));
}
#endif
};
}
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif
#line 9 "components/bitwise_convolution_modular.hpp"
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif
namespace bitwise_convolution_detail {
struct Modular : mint_detail::Ops {
static_assert(sizeof(mint) == sizeof(U) && std::is_trivially_copyable_v<mint>);
#if defined(__x86_64__) || defined(__i386__)
__attribute__((target("avx2"))) static I8 load8(const mint *a) {
return mint_detail::Ops::load8(reinterpret_cast<const U *>(a));
}
__attribute__((target("avx2"))) static void store8(mint *a, I8 value) {
mint_detail::Ops::store8(reinterpret_cast<U *>(a), value);
}
template <int kind, bool inverse>
__attribute__((target("avx2"))) static void pair(I8 &x, I8 &y) {
if constexpr (kind == 2) {
I8 old = x;
x = add8(x, y), y = sub8(old, y);
} else if constexpr (kind == 0) x = inverse ? sub8(x, y) : add8(x, y);
else y = inverse ? sub8(y, x) : add8(y, x);
}
template <int kind, bool inverse, int mask>
__attribute__((target("avx2"))) static I8 lanes(I8 x, I8 y) {
if constexpr (kind == 2) return _mm256_blend_epi32(add8(x, y), sub8(y, x), mask);
else return _mm256_blend_epi32(x, inverse ? sub8(x, y) : add8(x, y), mask);
}
template <int kind, bool inverse>
__attribute__((target("avx2"))) static void layers(mint *a, std::size_t start, std::size_t n, std::size_t half) {
for (; half < n / 4; half *= 8)
for (std::size_t block = start; block < start + n; block += 8 * half)
for (std::size_t i = block; i < block + half; i += 8) {
I8 x[8];
#pragma GCC unroll 8
for (int j = 0; j < 8; ++j) x[j] = load8(a + i + j * half);
#pragma GCC unroll 3
for (int h = 1; h <= 4; h *= 2)
#pragma GCC unroll 8
for (int j = 0; j < 8; ++j) if (!(j & h)) pair<kind, inverse>(x[j], x[j + h]);
#pragma GCC unroll 8
for (int j = 0; j < 8; ++j) store8(a + i + j * half, x[j]);
}
for (; half < n / 2; half *= 4)
for (std::size_t block = start; block < start + n; block += 4 * half)
for (std::size_t i = block; i < block + half; i += 8) {
I8 x = load8(a + i), y = load8(a + i + half), z = load8(a + i + 2 * half), w = load8(a + i + 3 * half);
pair<kind, inverse>(x, y), pair<kind, inverse>(z, w);
pair<kind, inverse>(x, z), pair<kind, inverse>(y, w);
store8(a + i, x), store8(a + i + half, y), store8(a + i + 2 * half, z), store8(a + i + 3 * half, w);
}
if (half < n)
for (std::size_t i = start; i < start + half; i += 8) {
I8 x = load8(a + i), y = load8(a + i + half);
pair<kind, inverse>(x, y);
store8(a + i, x), store8(a + i + half, y);
}
}
template <int kind, bool inverse>
__attribute__((target("avx2"))) static void transform(mint *a, std::size_t n) {
constexpr std::size_t block = 1024;
for (std::size_t start = 0; start < n; start += block) {
std::size_t size = std::min(block, n - start);
for (std::size_t i = start; i < start + size; i += 8) {
I8 x = load8(a + i);
x = lanes<kind, inverse, kind == 0 ? 0x55 : 0xaa>(x, _mm256_shuffle_epi32(x, 0xb1));
x = lanes<kind, inverse, kind == 0 ? 0x33 : 0xcc>(x, _mm256_shuffle_epi32(x, 0x4e));
x = lanes<kind, inverse, kind == 0 ? 0x0f : 0xf0>(x, _mm256_permute2x128_si256(x, x, 1));
store8(a + i, x);
}
layers<kind, inverse>(a, start, size, 8);
}
if (n > block) layers<kind, inverse>(a, 0, n, block);
}
template <int kind>
__attribute__((target("avx2"))) static void multiply(std::vector<mint> &a, std::vector<mint> &b) {
std::size_t n = a.size();
auto *x = a.data(), *y = b.data();
transform<kind, false>(x, n), transform<kind, false>(y, n);
I8 factor{};
if constexpr (kind == 2) factor = _mm256_set1_epi32(mint(n).inv().v);
for (std::size_t i = 0; i < n; i += 8) {
I8 value = mul8(load8(x + i), load8(y + i));
if constexpr (kind == 2) value = mul8(value, factor);
store8(x + i, value);
}
transform<kind, true>(x, n);
}
#endif
};
}
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif
#line 12 "convolution/bitwise_convolution.hpp"
#include <cstdint>
#line 14 "convolution/bitwise_convolution.hpp"
#include <utility>
#line 16 "convolution/bitwise_convolution.hpp"
namespace bitwise_convolution_detail {
template <int kind, class T>
std::vector<T> multiply(std::vector<T> a, std::vector<T> b) {
if (a.empty() || b.empty()) return {};
std::size_t n = std::bit_ceil(std::max(a.size(), b.size()));
a.resize(n, T(0));
b.resize(n, T(0));
#if (defined(__x86_64__) || defined(__i386__)) && !defined(BITWISE_CONVOLUTION_DISABLE_AVX2)
if constexpr (std::is_same_v<T, mint>) {
if (n >= 8 && Modular::vectorized()) {
Modular::multiply<kind>(a, b);
return a;
}
}
#endif
if constexpr (kind == 2) {
walsh_hadamard(a);
walsh_hadamard(b);
} else {
bitwise_zeta<false, kind == 0>(a);
bitwise_zeta<false, kind == 0>(b);
}
for (std::size_t i = 0; i < n; ++i) a[i] *= b[i];
if constexpr (kind == 2) {
walsh_hadamard(a);
if constexpr (std::numeric_limits<T>::is_integer)
for (T &x : a) x /= std::int64_t(n);
else {
T inverse = T(1) / T(n);
for (T &x : a) x *= inverse;
}
} else bitwise_zeta<true, kind == 0>(a);
return a;
}
}
template <class T>
std::vector<T> and_convolution(std::vector<T> a, std::vector<T> b) {
return bitwise_convolution_detail::multiply<0>(std::move(a), std::move(b));
}
template <class T>
std::vector<T> or_convolution(std::vector<T> a, std::vector<T> b) {
return bitwise_convolution_detail::multiply<1>(std::move(a), std::move(b));
}
template <class T>
std::vector<T> xor_convolution(std::vector<T> a, std::vector<T> b) {
return bitwise_convolution_detail::multiply<2>(std::move(a), std::move(b));
}
#line 1 "convolution/gcd_lcm_convolution.hpp"
#line 1 "components/gcd_convolution_mint.hpp"
#line 1 "components/divisor_zeta.hpp"
#line 1 "math/prime_sieve.hpp"
#line 5 "math/prime_sieve.hpp"
#include <array>
#line 7 "math/prime_sieve.hpp"
#include <cmath>
#line 10 "math/prime_sieve.hpp"
#include <cstring>
#line 12 "math/prime_sieve.hpp"
inline std::vector<std::size_t> prime_sieve(std::size_t n);
template<class F> void for_each_prime(std::size_t n, F emit) {
if (n <= 2) return;
emit(std::size_t(2));
if (n <= 65537) {
const std::size_t size = n / 2;
std::vector<std::uint8_t> prime(size, 1);
prime[0] = 0;
for (std::size_t p = 3; p <= (n - 1) / p; p += 2) {
if (!prime[p / 2]) continue;
for (std::size_t j = p * p / 2; j < size; j += p) prime[j] = 0;
}
std::size_t i = 0;
#if defined(__BYTE_ORDER__) && __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
for (; i + 8 <= size; i += 8) {
std::uint64_t word;
std::memcpy(&word, prime.data() + i, sizeof(word));
while (word) {
emit(2 * (i + (std::countr_zero(word) >> 3)) + 1);
word &= word - 1;
}
}
#endif
for (; i < size; ++i) if (prime[i]) emit(2 * i + 1);
return;
}
emit(3), emit(5);
constexpr std::array<unsigned, 8> residues{1, 7, 11, 13, 17, 19, 23, 29};
constexpr auto index = [&] {
std::array<unsigned, 30> index{};
for (unsigned i = 0; i < 8; ++i) index[residues[i]] = i;
return index;
}();
constexpr std::size_t period = 7 * 11 * 13, block_size = period * 32;
static const auto pattern = [&] {
std::array<std::uint8_t, period> pattern{};
for (std::size_t i = 0; i < period; ++i)
for (unsigned d = 0; d < 8; ++d) {
const auto x = 30 * i + residues[d];
pattern[i] |= std::uint8_t(bool(x % 7 && x % 11 && x % 13) << d);
}
return pattern;
}();
std::size_t root = std::sqrt(static_cast<long double>(n - 1));
while (root > (n - 1) / root) --root;
while (root + 1 <= (n - 1) / (root + 1)) ++root;
const auto primes = prime_sieve(root + 1);
struct State {
std::size_t p;
std::array<std::size_t, 8> next;
std::array<std::uint8_t, 8> mask;
};
std::vector<State> state;
state.reserve(primes.size());
for (const auto p : primes) {
if (p < 17) continue;
State current;
current.p = p;
for (unsigned d = 0; d < 8; ++d) {
const std::size_t k = p + (residues[d] + 30 - p % 30) % 30;
const __uint128_t x = __uint128_t(p) * k;
current.next[d] = x / 30;
current.mask[d] = std::uint8_t(~(1u << index[x % 30]));
}
state.push_back(current);
}
const std::size_t size = (n - 1) / 30 + 1;
std::vector<std::uint8_t> prime(block_size);
constexpr auto offsets = [&] {
std::array<unsigned, 64> offsets{};
for (unsigned i = 0; i < 64; ++i) offsets[i] = 30 * (i / 8) + residues[i % 8];
return offsets;
}();
for (std::size_t base = 0; base < size;) {
const std::size_t length = std::min(block_size, size - base);
for (std::size_t i = 0; i < length; i += period)
std::copy_n(pattern.data(), std::min(period, length - i), prime.data() + i);
if (base == 0) {
prime[0] &= ~1u;
prime[0] |= (1u << 1) | (1u << 2) | (1u << 3);
}
for (auto ¤t : state) {
const auto p = current.p;
for (unsigned d = 0; d < 8; ++d) {
std::size_t j = current.next[d];
const auto mask = current.mask[d];
for (; j + 3 * p < length; j += 4 * p) {
prime[j] &= mask;
prime[j + p] &= mask;
prime[j + 2 * p] &= mask;
prime[j + 3 * p] &= mask;
}
for (; j < length; j += p) prime[j] &= mask;
current.next[d] = j - length;
}
}
if (base + length == size) {
const auto remaining = n - 30 * (size - 1);
for (unsigned d = 0; d < 8; ++d) if (residues[d] >= remaining) prime[length - 1] &= ~(1u << d);
}
std::fill(prime.begin() + length, prime.begin() + ((length + 7) / 8 * 8), 0);
for (std::size_t i = 0; i < length; i += 8) {
std::uint64_t word = 0;
std::memcpy(&word, prime.data() + i, sizeof(word));
if constexpr (std::endian::native == std::endian::big) word = std::byteswap(word);
else if constexpr (std::endian::native != std::endian::little) {
word = 0;
for (unsigned d = 0; d < 8; ++d) word |= std::uint64_t(prime[i + d]) << (8 * d);
}
while (word) {
const auto x = 30 * (base + i) + offsets[std::countr_zero(word)];
emit(x);
word &= word - 1;
}
}
base += length;
}
}
inline std::vector<std::size_t> prime_sieve(std::size_t n) {
std::vector<std::size_t> result;
if (n >= 1024) result.reserve(std::size_t(double(n) / (std::log(double(n)) - 1.1)));
for_each_prime(n, [&](std::size_t p) { result.push_back(p); });
return result;
}
#line 7 "components/divisor_zeta.hpp"
namespace divisor_zeta_detail {
template <bool inverse, bool multiples, class F>
void for_each(std::size_t n, const std::vector<std::size_t> &ps, std::size_t count, F update) {
if (!n) return;
for (std::size_t j = 0; j < count; ++j) {
const std::size_t p = ps[j], end = (n - 1) / p;
auto apply = [&](std::size_t i) { update(multiples ? i : i * p, multiples ? i * p : i); };
if constexpr (inverse == multiples) {
std::size_t i = 1;
for (; i + 3 <= end; i += 4) apply(i), apply(i + 1), apply(i + 2), apply(i + 3);
for (; i <= end; ++i) apply(i);
} else {
std::size_t i = end;
for (; i >= 4; i -= 4) apply(i), apply(i - 1), apply(i - 2), apply(i - 3);
for (; i; --i) apply(i);
}
}
}
}
template <bool inverse = false, bool multiples = false, class T>
void divisor_zeta(std::vector<T> &a, const std::vector<std::size_t> &ps) {
divisor_zeta_detail::for_each<inverse, multiples>(a.size(), ps, ps.size(), [&](std::size_t to, std::size_t from) {
if constexpr (inverse) a[to] -= a[from];
else a[to] += a[from];
});
}
template <bool inverse = false, bool multiples = false, class T>
void divisor_zeta(std::vector<T> &a) { divisor_zeta<inverse, multiples>(a, prime_sieve(a.size())); }
#line 10 "components/gcd_convolution_mint.hpp"
#if defined(__x86_64__) || defined(__i386__)
#include <immintrin.h>
#endif
namespace mint_detail { struct Mint; }
namespace gcd_convolution_mint_detail {
#if defined(__x86_64__) || defined(__i386__)
template <bool inverse, class T>
__attribute__((target("avx2")))
void prime_pass(T *a, T *b, std::size_t n, const std::vector<std::size_t> &ps, std::size_t count) {
const __m256i mod = _mm256_set1_epi32(2 * T::mod());
for (std::size_t j = 0; j < count; ++j) {
const std::size_t p = ps[j], end = (n - 1) / p;
auto scalar = [&](std::size_t i) {
if constexpr (inverse) a[i] -= a[i * p];
else a[i] += a[i * p], b[i] += b[i * p];
};
auto gather = [&](const T *s) __attribute__((target("avx2"))) {
return _mm256_setr_epi32(s[0].v, s[p].v, s[2 * p].v, s[3 * p].v,
s[4 * p].v, s[5 * p].v, s[6 * p].v, s[7 * p].v);
};
auto vector = [&](std::size_t i) __attribute__((target("avx2"))) {
__m256i x = _mm256_loadu_si256(reinterpret_cast<const __m256i *>(a + i)), y = gather(a + i * p);
if constexpr (inverse) {
x = _mm256_sub_epi32(x, y);
x = _mm256_min_epu32(x, _mm256_add_epi32(x, mod));
} else {
x = _mm256_add_epi32(x, y);
x = _mm256_min_epu32(x, _mm256_sub_epi32(x, mod));
}
_mm256_storeu_si256(reinterpret_cast<__m256i *>(a + i), x);
if constexpr (!inverse) {
x = _mm256_add_epi32(_mm256_loadu_si256(reinterpret_cast<const __m256i *>(b + i)), gather(b + i * p));
x = _mm256_min_epu32(x, _mm256_sub_epi32(x, mod));
_mm256_storeu_si256(reinterpret_cast<__m256i *>(b + i), x);
}
};
if constexpr (inverse) {
std::size_t i = 1;
for (; i <= end && i < 8; ++i) scalar(i);
for (; i + 7 <= end; i += 8) vector(i);
for (; i <= end; ++i) scalar(i);
} else {
std::size_t i = end;
for (; i >= 15; i -= 8) vector(i - 7);
for (; i; --i) scalar(i);
}
}
}
#endif
template <bool inverse, class T>
void transform(std::vector<T> &a, std::vector<T> &b, const std::vector<std::size_t> &ps, std::size_t count) {
#if defined(__x86_64__) || defined(__i386__)
if (__builtin_cpu_supports("avx2")) {
prime_pass<inverse>(a.data(), b.data(), a.size(), ps, count);
return;
}
#endif
divisor_zeta_detail::for_each<inverse, true>(a.size(), ps, count, [&](std::size_t to, std::size_t from) {
if constexpr (inverse) a[to] -= a[from];
else a[to] += a[from], b[to] += b[from];
});
}
template <bool inverse, class T>
void large_multiples(std::vector<T> &a, std::vector<T> &b, const std::vector<std::size_t> &ps,
std::size_t first) {
if (first == ps.size()) return;
constexpr std::uint64_t twice = 2 * T::mod();
const std::size_t n = a.size() - 1;
for (std::size_t i = 1; i <= n / ps[first]; ++i) {
const auto end = std::upper_bound(ps.begin() + first, ps.end(), n / i);
std::uint64_t x = inverse ? 0 : a[i].v, y = inverse ? a[i].v : b[i].v;
for (auto p = ps.begin() + first; p != end; ++p) {
x += a[i * *p].v;
if constexpr (!inverse) y += b[i * *p].v;
}
if constexpr (inverse) a[i] = T::word(y % twice) - T::word(x % twice);
else {
a[i] = T::word(x % twice), b[i] = T::word(y % twice);
}
}
}
template <class T>
void multiply(std::vector<T> &a, std::vector<T> &b, const std::vector<std::size_t> &ps) {
const std::size_t n = a.size() - 1;
std::size_t cutoff = std::sqrt(n);
while (cutoff + 1 <= n / (cutoff + 1)) ++cutoff;
while (cutoff && cutoff > n / cutoff) --cutoff;
const std::size_t first = std::upper_bound(ps.begin(), ps.end(), cutoff) - ps.begin();
transform<false>(a, b, ps, first);
large_multiples<false>(a, b, ps, first);
for (std::size_t i = 1; i <= n; ++i) a[i] *= b[i];
transform<true>(a, b, ps, first);
large_multiples<true>(a, b, ps, first);
}
}
#line 12 "convolution/gcd_lcm_convolution.hpp"
namespace gcd_lcm_convolution_detail {
template <bool gcd, class T>
std::vector<T> multiply(std::vector<T> a, std::vector<T> b) {
if (a.empty() || b.empty()) return {};
assert(a[0] == T(0) && b[0] == T(0));
std::size_t n = std::max(a.size(), b.size());
a.resize(n, T(0));
b.resize(n, T(0));
auto ps = prime_sieve(n);
if constexpr (gcd && std::is_same_v<T, mint_detail::Mint>) {
if (n >= 4096 && n - 1 <= UINT32_MAX) {
gcd_convolution_mint_detail::multiply(a, b, ps);
return a;
}
}
divisor_zeta_detail::for_each<false, gcd>(n, ps, ps.size(), [&](std::size_t to, std::size_t from) {
a[to] += a[from], b[to] += b[from];
});
for (std::size_t i = 1; i < n; ++i) a[i] *= b[i];
divisor_zeta<true, gcd>(a, ps);
return a;
}
}
template <class T>
std::vector<T> gcd_convolution(std::vector<T> a, std::vector<T> b) {
return gcd_lcm_convolution_detail::multiply<true>(std::move(a), std::move(b));
}
template <class T>
std::vector<T> lcm_convolution(std::vector<T> a, std::vector<T> b) {
return gcd_lcm_convolution_detail::multiply<false>(std::move(a), std::move(b));
}
#line 1 "convolution/subset_convolution.hpp"
#line 1 "components/subset_convolution_modular.hpp"
#line 11 "components/subset_convolution_modular.hpp"
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
#include <immintrin.h>
#endif
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif
namespace subset_convolution_detail {
template <class T, class = void>
struct SmallModint : std::false_type {};
template <class T>
struct SmallModint<T, std::void_t<decltype(std::integral_constant<std::uint32_t, T::mod()>{}),
decltype(std::declval<const T &>().val()), decltype(T::raw(std::uint32_t{}))>>
: std::bool_constant<(T::mod() > 1 && T::mod() < (std::uint32_t(1) << 30) && (T::mod() & 1))
&& std::is_convertible_v<decltype(std::declval<const T &>().val()), std::uint32_t>
&& std::is_same_v<decltype(T::raw(std::uint32_t{})), T>> {};
inline bool avx2() {
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__)) && !defined(SUBSET_CONVOLUTION_DISABLE_AVX2)
static const bool enabled = __builtin_cpu_supports("avx2");
return enabled;
#else
return false;
#endif
}
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
template <std::uint32_t mod>
struct Modular {
using V = std::uint32_t;
static V add(V a, V b) { V c = a + b; return c < mod ? c : c - mod; }
static V sub(V a, V b) { return a >= b ? a - b : a + mod - b; }
static constexpr int width = 24;
__attribute__((target("avx2"))) static void transform(V *a, int n) {
const __m256i p = _mm256_set1_epi32(mod);
if (n >= 128) {
int q = n / 4;
for (int i = 0; i < q; ++i)
for (int rank = 0; rank < width; rank += 8) {
auto *p0 = reinterpret_cast<__m256i *>(a + std::size_t(i) * width + rank);
auto *p1 = reinterpret_cast<__m256i *>(a + std::size_t(i + q) * width + rank);
auto *p2 = reinterpret_cast<__m256i *>(a + std::size_t(i + q * 2) * width + rank);
auto *p3 = reinterpret_cast<__m256i *>(a + std::size_t(i + q * 3) * width + rank);
__m256i x0 = _mm256_loadu_si256(p0), x1 = _mm256_loadu_si256(p1), x2 = _mm256_loadu_si256(p2), x3 = _mm256_loadu_si256(p3);
__m256i u = _mm256_add_epi32(x0, x2), v = _mm256_add_epi32(_mm256_sub_epi32(x0, x2), p);
x0 = _mm256_min_epu32(u, _mm256_sub_epi32(u, p));
x2 = _mm256_min_epu32(v, _mm256_sub_epi32(v, p));
u = _mm256_add_epi32(x1, x3);
v = _mm256_add_epi32(_mm256_sub_epi32(x1, x3), p);
x1 = _mm256_min_epu32(u, _mm256_sub_epi32(u, p));
x3 = _mm256_min_epu32(v, _mm256_sub_epi32(v, p));
u = _mm256_add_epi32(x0, x1);
v = _mm256_add_epi32(_mm256_sub_epi32(x0, x1), p);
x0 = _mm256_min_epu32(u, _mm256_sub_epi32(u, p));
x1 = _mm256_min_epu32(v, _mm256_sub_epi32(v, p));
u = _mm256_add_epi32(x2, x3);
v = _mm256_add_epi32(_mm256_sub_epi32(x2, x3), p);
x2 = _mm256_min_epu32(u, _mm256_sub_epi32(u, p));
x3 = _mm256_min_epu32(v, _mm256_sub_epi32(v, p));
_mm256_storeu_si256(p0, x0);
_mm256_storeu_si256(p1, x1);
_mm256_storeu_si256(p2, x2);
_mm256_storeu_si256(p3, x3);
}
for (int i = 0; i < 4; ++i) transform(a + std::size_t(i) * q * width, q);
return;
}
for (int half = 1; half < n; half *= 2)
for (int start = 0; start < n; start += 2 * half)
for (int j = 0; j < half; ++j)
for (int rank = 0; rank < width; rank += 8) {
auto *left = reinterpret_cast<__m256i *>(a + std::size_t(start + j) * width + rank);
auto *right = reinterpret_cast<__m256i *>(a + std::size_t(start + half + j) * width + rank);
__m256i x = _mm256_loadu_si256(left), y = _mm256_loadu_si256(right);
__m256i u = _mm256_add_epi32(x, y), v = _mm256_add_epi32(_mm256_sub_epi32(x, y), p);
_mm256_storeu_si256(left, _mm256_min_epu32(u, _mm256_sub_epi32(u, p)));
_mm256_storeu_si256(right, _mm256_min_epu32(v, _mm256_sub_epi32(v, p)));
}
}
static constexpr V inverse = [] {
V x = 1;
for (int i = 0; i < 5; ++i) x *= 2 - mod * x;
return V(0) - x;
}();
__attribute__((target("avx2"))) static __m256i reduce(__m256i x) {
const __m256i p = _mm256_set1_epi64x(mod), inv = _mm256_set1_epi64x(inverse);
__m256i q = _mm256_mul_epu32(x, inv);
x = _mm256_srli_epi64(_mm256_add_epi64(x, _mm256_mul_epu32(q, p)), 32);
x = _mm256_min_epu32(x, _mm256_sub_epi64(x, p));
return _mm256_min_epu32(x, _mm256_sub_epi64(x, p));
}
__attribute__((target("avx2"))) static void product(V *a, const V *b, int n) {
constexpr V r2 = V((~std::uint64_t(0) % mod + 1) % mod);
const __m256i p = _mm256_set1_epi64x(mod), scale = _mm256_set1_epi64x(r2);
for (int pos = 0; pos < n; pos += 4) {
__m256i x[20], y[20], out[20];
for (int rank = 0; rank < 20; rank += 4) {
for (int side = 0; side < 2; ++side) {
const V *v = (side ? b : a) + std::size_t(pos) * width + rank;
__m128i a0 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(v));
__m128i a1 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(v + width));
__m128i a2 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(v + width * 2));
__m128i a3 = _mm_loadu_si128(reinterpret_cast<const __m128i *>(v + width * 3));
__m128i u = _mm_unpacklo_epi32(a0, a1), v0 = _mm_unpackhi_epi32(a0, a1), w = _mm_unpacklo_epi32(a2, a3), z = _mm_unpackhi_epi32(a2, a3);
auto *dest = side ? y : x;
dest[rank] = _mm256_cvtepu32_epi64(_mm_unpacklo_epi64(u, w));
dest[rank + 1] = _mm256_cvtepu32_epi64(_mm_unpackhi_epi64(u, w));
dest[rank + 2] = _mm256_cvtepu32_epi64(_mm_unpacklo_epi64(v0, z));
dest[rank + 3] = _mm256_cvtepu32_epi64(_mm_unpackhi_epi64(v0, z));
}
}
out[0] = _mm256_setzero_si256();
#pragma GCC unroll 20
for (int rank = 1; rank < 20; ++rank) {
__m256i value = _mm256_setzero_si256();
for (int first = 0; first < rank; first += 8) {
__m256i sum = _mm256_setzero_si256();
for (int i = first; i < std::min(rank, first + 8); ++i)
sum = _mm256_add_epi64(sum, _mm256_mul_epu32(x[i], y[rank - i - 1]));
value = _mm256_add_epi64(value, reduce(sum));
value = _mm256_min_epu32(value, _mm256_sub_epi64(value, p));
}
out[rank] = reduce(_mm256_mul_epu32(value, scale));
}
for (int rank = 0; rank < 20; rank += 4) {
__m128i z[4];
for (int j = 0; j < 4; ++j) z[j] = _mm256_castsi256_si128(_mm256_permutevar8x32_epi32(out[rank + j], _mm256_setr_epi32(0, 2, 4, 6, 0, 2, 4, 6)));
__m128i u = _mm_unpacklo_epi32(z[0], z[1]), v = _mm_unpackhi_epi32(z[0], z[1]), w = _mm_unpacklo_epi32(z[2], z[3]), zz = _mm_unpackhi_epi32(z[2], z[3]);
V *dest = a + std::size_t(pos) * width + rank;
_mm_storeu_si128(reinterpret_cast<__m128i *>(dest), _mm_unpacklo_epi64(u, w));
_mm_storeu_si128(reinterpret_cast<__m128i *>(dest + width), _mm_unpackhi_epi64(u, w));
_mm_storeu_si128(reinterpret_cast<__m128i *>(dest + width * 2), _mm_unpacklo_epi64(v, zz));
_mm_storeu_si128(reinterpret_cast<__m128i *>(dest + width * 3), _mm_unpackhi_epi64(v, zz));
}
}
}
template <class T>
__attribute__((target("avx2"))) static std::vector<T> convolution(const std::vector<T> &a, const std::vector<T> &b, int n) {
int levels = std::bit_width(unsigned(n)) - 1;
int bottom = n >> (std::min(levels - 3, 2) + 1);
std::vector<V> x(std::size_t(bottom) * width), y(x.size()), counts(n), result(n);
std::vector<V> aa(n), bb(n);
for (int i = 0; i < int(a.size()); ++i) aa[i] = a[i].val();
for (int i = 0; i < int(b.size()); ++i) bb[i] = b[i].val();
for (int i = 1; i < n; ++i) counts[i] = std::popcount(unsigned(i)) - 1;
for (int top = 0; top < n / 2; top += bottom) {
std::fill(x.begin(), x.end(), 0);
std::fill(y.begin(), y.end(), 0);
for (int start = 0; start < n; start += 2 * bottom) {
bool negative = std::popcount(unsigned((start >> 1) & top)) & 1;
for (int mask = start == 0; mask < 2 * bottom; ++mask) {
std::size_t index = std::size_t(mask >> 1) * width + counts[start + mask];
if (negative) x[index] = sub(x[index], aa[start + mask]), y[index] = sub(y[index], bb[start + mask]);
else x[index] = add(x[index], aa[start + mask]), y[index] = add(y[index], bb[start + mask]);
}
}
transform(x.data(), bottom);
transform(y.data(), bottom);
product(x.data(), y.data(), bottom);
transform(x.data(), bottom);
for (int start = 0; start < n; start += 2 * bottom) {
bool negative = std::popcount(unsigned((start >> 1) & top)) & 1;
for (int mask = start == 0; mask < 2 * bottom; ++mask) {
V value = x[std::size_t(mask >> 1) * width + counts[start + mask]];
result[start + mask] = negative ? sub(result[start + mask], value) : add(result[start + mask], value);
}
}
}
V scale = 1;
for (int i = 1; i < levels; ++i) scale = std::uint64_t(scale) * ((mod + 1) / 2) % mod;
std::vector<T> out(n, T(0));
out[0] = T::raw(std::uint64_t(aa[0]) * bb[0] % mod);
for (int i = 1; i < n; ++i) out[i] = T::raw((std::uint64_t(result[i]) * scale + std::uint64_t(aa[i]) * bb[0] + std::uint64_t(aa[0]) * bb[i]) % mod);
return out;
}
};
#endif
}
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif
#line 9 "convolution/subset_convolution.hpp"
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif
template <class T>
std::vector<T> subset_convolution(const std::vector<T> &a, const std::vector<T> &b) {
if (a.empty() || b.empty()) return {};
std::size_t n = std::bit_ceil(std::max(a.size(), b.size()));
#if (defined(__x86_64__) || defined(__i386__)) && (defined(__GNUC__) || defined(__clang__))
if constexpr (subset_convolution_detail::SmallModint<T>::value)
if (n >= 256 && n <= (std::size_t(1) << 20) && subset_convolution_detail::avx2())
return subset_convolution_detail::Modular<T::mod()>::convolution(a, b, int(n));
#endif
int levels = std::bit_width(n) - 1, width = levels + 1;
std::vector<T> x(n * width, T(0)), y(n * width, T(0)), result(n, T(0));
for (std::size_t i = 0; i < a.size(); ++i) x[i * width + std::popcount(i)] = a[i];
for (std::size_t i = 0; i < b.size(); ++i) y[i * width + std::popcount(i)] = b[i];
auto zeta = [&](std::size_t first, std::size_t length, std::size_t half) {
for (; half < length; half *= 2)
for (std::size_t start = first; start < first + length; start += 2 * half) {
int low = std::popcount(start);
for (std::size_t j = start; j < start + half; ++j)
for (int k = low, end = std::popcount(j); k <= end; ++k) {
x[(j + half) * width + k] += x[j * width + k];
y[(j + half) * width + k] += y[j * width + k];
}
}
};
std::size_t block = std::min<std::size_t>(n, 256);
for (std::size_t start = 0; start < n; start += block) zeta(start, block, 1);
zeta(0, n, block);
for (std::size_t mask = 0; mask < n; ++mask) {
int count = std::popcount(mask);
for (int k = std::min(levels, 2 * count); k >= count; --k) {
T sum(0);
for (int i = k - count; i <= count; ++i)
sum += x[mask * width + i] * y[mask * width + k - i];
x[mask * width + k] = sum;
}
}
auto mobius = [&](std::size_t first, std::size_t length, std::size_t half) {
for (; half < length; half *= 2)
for (std::size_t start = first; start < first + length; start += 2 * half) {
int low = std::popcount(start) + 1, high = levels - std::countr_zero(half);
for (std::size_t j = start; j < start + half; ++j) {
int count = std::popcount(j - start);
for (int k = low + count; k <= high + count; ++k)
x[(j + half) * width + k] -= x[j * width + k];
}
}
};
for (std::size_t start = 0; start < n; start += block) mobius(start, block, 1);
mobius(0, n, block);
for (std::size_t mask = 0; mask < n; ++mask) result[mask] = x[mask * width + std::popcount(mask)];
return result;
}
#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif
#line 15 "convolution/boolean_convolution.hpp"
namespace boolean_convolution_detail {
template <int index, bool dual = false, class T>
std::vector<T> multiply(const std::vector<T> &a, const std::vector<T> &b) {
static_assert(std::is_integral_v<T> && sizeof(T) <= 8 && 0 <= index && index <= 5);
if (a.empty() || b.empty()) return {};
using U = std::make_unsigned_t<std::conditional_t<std::is_same_v<T, bool>, unsigned char, T>>;
std::size_t n = std::max(a.size(), b.size());
if constexpr (index < 3 || index == 5) n = std::bit_ceil(n);
assert(n <= std::size_t(1) << 30);
constexpr U full = std::is_same_v<T, bool> ? U(1) : U(~U(0));
U mask = 0, active = 0;
for (std::size_t i = index == 3 || index == 4; i < a.size() && mask != full; ++i) {
U value = U(a[i]);
mask |= U(dual ? ~value : value) & full;
}
for (std::size_t i = index == 3 || index == 4; i < b.size() && active != mask; ++i) {
U value = U(b[i]);
active |= U(dual ? ~value : value) & mask;
}
std::vector<T> result(n, dual ? T(~U(0)) : T(0));
for (; active; active &= active - 1) {
int bit = std::countr_zero(active);
std::vector<__int128_t> x(n), y(n);
for (std::size_t i = index == 3 || index == 4; i < a.size(); ++i) x[i] = ((U(a[i]) >> bit) & 1) ^ dual;
for (std::size_t i = index == 3 || index == 4; i < b.size(); ++i) y[i] = ((U(b[i]) >> bit) & 1) ^ dual;
std::vector<__int128_t> counts;
if constexpr (index < 3) counts = bitwise_convolution_detail::multiply<index>(std::move(x), std::move(y));
else if constexpr (index == 5) counts = subset_convolution(std::move(x), std::move(y));
else counts = gcd_lcm_convolution_detail::multiply<index == 3>(std::move(x), std::move(y));
for (std::size_t i = 0; i < n; ++i)
if (counts[i]) {
U flag = U(1) << bit;
result[i] = T(dual ? U(result[i]) & ~flag : U(result[i]) | flag);
}
}
if constexpr (index == 3 || index == 4) result[0] = T(0);
return result;
}
}
template <class T>
std::vector<T> and_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<0>(a, b); }
template <class T>
std::vector<T> or_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<1>(a, b); }
template <class T>
std::vector<T> xor_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<2>(a, b); }
template <class T>
std::vector<T> gcd_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<3>(a, b); }
template <class T>
std::vector<T> lcm_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<4>(a, b); }
template <class T>
std::vector<T> and_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<0, true>(a, b); }
template <class T>
std::vector<T> or_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<1, true>(a, b); }
template <class T>
std::vector<T> xor_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<2, true>(a, b); }
template <class T>
std::vector<T> gcd_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<3, true>(a, b); }
template <class T>
std::vector<T> lcm_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<4, true>(a, b); }
template <class T>
std::vector<T> subset_or_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<5>(a, b); }
template <class T>
std::vector<T> subset_and_convolution(const std::vector<T> &a, const std::vector<T> &b) { return boolean_convolution_detail::multiply<5, true>(a, b); }