PROBLEM LINK:
Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4
Author: raysh07
Tester: iceknight1093
Editorialist: iceknight1093
DIFFICULTY:
TBD
PREREQUISITES:
NTT
PROBLEM:
You’re given a permutation P.
For each i = 0, 1, \ldots, N count the number of pairs of arrays (A, B) such that |A| = i and merge(A, B) = P, where merge(A, B) is the result of greedily merging A and B (i.e. always take the smaller element among the first elements of A and B.)
N \le 2\cdot 10^5
EXPLANATION:
From the easy version, we know the solution is as follows:
- Break the permutation P into blocks, where a block is formed by taking the maximum remaining element and everything after it.
- Then, if the block sizes are \{s_1, \ldots, s_k\}, the answer for |A|=x is the number of subsets of that multiset whose sum equals x.
In the easy version, we could use a knapsack DP to do the second part; but with higher constraints that’s no longer viable.
Instead of doing DP, we will use generating functions.
Specifically, to a block of size s_i we’ll associate the generating function (1 + x^{s_i}), where the 1 denotes not picking the block (equivalently, sending the block to array B) and x^{s_i} denotes picking the block for array A.
Consider the product of all these generating functions, i.e.
p(x) = \prod_{i=1}^k \left(1 + x^{s_i}\right)
For any integer m, the coefficient of x^m in p(x) will tell us exactly the number of ways of obtaining m as the sum of some subset of the x_i, due to how polynomial multiplication works.
Thus, our goal is to compute the first N+1 coefficients of p(x).
(Because \sum s_i = N, coefficients of higher degree terms than N in p(x) are 0 anyway.)
Or, putting it simply, we simply want to find the polynomial p(x) quickly.
Since p(x) is the product of several polynomials, we can simply use a fast polynomial multiplication algorithm (here, number theoretic transform) to do this quickly; but since we can have \mathcal{O}(N) polynomials we need to be careful about how we do it.
There are two ways:
- Start with all the polynomials (1+x^{s_i}).
Repeatedly extract the two lowest-degree polynomials, multiply them, and then insert their product into the set.
This is essentially Huffman encoding, and since \sum s_i =N the overall complexity is \mathcal{O}(N\log^2 N).
(The extra log is because multiplying two polynomials of sizes x and y is done in \mathcal{O}((x+y)\log(x+y)).) - Alternately, you can do a divide-and conquer.
Split the polynomials into two halves of approximately equal size, then recursively find the product of each half and do a final multiplication to obtain their product.
With k polynomials having total degree N, this is \mathcal{O}(N\log N\log k).
As long as NTT is used for the convolution, either method will work.
TIME COMPLEXITY:
\mathcal{O}(N \log^2 N) per testcase.
CODE:
Editorialist's code (C++)
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include "bits/stdc++.h"
#include <cassert>
using namespace std;
using ll = long long int;
mt19937_64 RNG(chrono::high_resolution_clock::now().time_since_epoch().count());
// https://judge.yosupo.jp/submission/69895
namespace ntt {
template <class T, class F = multiplies<T>>
T power(T a, long long n, F op = multiplies<T>(), T e = {1}) {
// assert(n >= 0);
T res = e;
while (n) {
if (n & 1) res = op(res, a);
if (n >>= 1) a = op(a, a);
}
return res;
}
constexpr int mod = int(1e9) + 7;
constexpr int nttmod = 998'244'353;
template <std::uint32_t P>
struct ModInt32 {
public:
using i32 = std::int32_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
using u64 = std::uint64_t;
using m32 = ModInt32;
using internal_value_type = u32;
private:
u32 v;
static constexpr u32 get_r() {
u32 iv = P;
for (u32 i = 0; i != 4; ++i) iv *= 2U - P * iv;
return -iv;
}
static constexpr u32 r = get_r(), r2 = -u64(P) % P;
static_assert((P & 1) == 1);
static_assert(-r * P == 1);
static_assert(P < (1 << 30));
static constexpr u32 pow_mod(u32 x, u64 y) {
u32 res = 1;
for (; y != 0; y >>= 1, x = u64(x) * x % P)
if (y & 1) res = u64(res) * x % P;
return res;
}
static constexpr u32 reduce(u64 x) {
return (x + u64(u32(x) * r) * P) >> 32;
}
static constexpr u32 norm(u32 x) { return x - (P & -(x >= P)); }
public:
static constexpr u32 get_pr() {
u32 tmp[32] = {}, cnt = 0;
const u64 phi = P - 1;
u64 m = phi;
for (u64 i = 2; i * i <= m; ++i)
if (m % i == 0) {
tmp[cnt++] = i;
while (m % i == 0) m /= i;
}
if (m != 1) tmp[cnt++] = m;
for (u64 res = 2; res != P; ++res) {
bool flag = true;
for (u32 i = 0; i != cnt && flag; ++i)
flag &= pow_mod(res, phi / tmp[i]) != 1;
if (flag) return res;
}
return 0;
}
constexpr ModInt32() : v(0){};
~ModInt32() = default;
constexpr ModInt32(u32 _v) : v(reduce(u64(_v) * r2)) {}
constexpr ModInt32(i32 _v) : v(reduce(u64(_v % P + P) * r2)) {}
constexpr ModInt32(u64 _v) : v(reduce((_v % P) * r2)) {}
constexpr ModInt32(i64 _v) : v(reduce(u64(_v % P + P) * r2)) {}
constexpr ModInt32(const m32& rhs) : v(rhs.v) {}
constexpr u32 get() const { return norm(reduce(v)); }
explicit constexpr operator u32() const { return get(); }
explicit constexpr operator i32() const { return i32(get()); }
constexpr m32& operator=(const m32& rhs) { return v = rhs.v, *this; }
constexpr m32 operator-() const {
m32 res;
return res.v = (P << 1 & -(v != 0)) - v, res;
}
constexpr m32 inv() const { return pow(P - 2); }
constexpr m32& operator+=(const m32& rhs) {
return v += rhs.v - (P << 1), v += P << 1 & -(v >> 31), *this;
}
constexpr m32& operator-=(const m32& rhs) {
return v -= rhs.v, v += P << 1 & -(v >> 31), *this;
}
constexpr m32& operator*=(const m32& rhs) {
return v = reduce(u64(v) * rhs.v), *this;
}
constexpr m32& operator/=(const m32& rhs) {
return this->operator*=(rhs.inv());
}
friend m32 operator+(const m32& lhs, const m32& rhs) {
return m32(lhs) += rhs;
}
friend m32 operator-(const m32& lhs, const m32& rhs) {
return m32(lhs) -= rhs;
}
friend m32 operator*(const m32& lhs, const m32& rhs) {
return m32(lhs) *= rhs;
}
friend m32 operator/(const m32& lhs, const m32& rhs) {
return m32(lhs) /= rhs;
}
friend bool operator==(const m32& lhs, const m32& rhs) {
return norm(lhs.v) == norm(rhs.v);
}
friend bool operator!=(const m32& lhs, const m32& rhs) {
return norm(lhs.v) != norm(rhs.v);
}
friend std::istream& operator>>(std::istream& is, m32& rhs) {
return is >> rhs.v, rhs.v = reduce(u64(rhs.v) * r2), is;
}
friend std::ostream& operator<<(std::ostream& os, const m32& rhs) {
return os << rhs.get();
}
constexpr m32 pow(i64 y) const {
// assumes P is a prime
i64 rem = y % (P - 1);
if (y > 0 && rem == 0)
y = P - 1;
else
y = rem;
m32 res(1), x(*this);
for (; y != 0; y >>= 1, x *= x)
if (y & 1) res *= x;
return res;
}
};
using mint = ModInt32<nttmod>;
void ntt(vector<mint>& a, bool inverse) {
static array<mint, 30> dw{}, idw{};
if (dw[0] == 0) {
mint root = 2;
while (power(root, (nttmod - 1) / 2) == 1) root += 1;
for (int i = 0; i < 30; ++i)
dw[i] = -power(root, (nttmod - 1) >> (i + 2)),
idw[i] = 1 / dw[i];
}
int n = (int)a.size();
assert((n & (n - 1)) == 0);
if (not inverse) {
for (int m = n; m >>= 1;) {
mint w = 1;
for (int s = 0, k = 0; s < n; s += 2 * m) {
for (int i = s, j = s + m; i < s + m; ++i, ++j) {
auto x = a[i], y = a[j] * w;
a[i] = x + y;
a[j] = x - y;
}
w *= dw[__builtin_ctz(++k)];
}
}
} else {
for (int m = 1; m < n; m *= 2) {
mint w = 1;
for (int s = 0, k = 0; s < n; s += 2 * m) {
for (int i = s, j = s + m; i < s + m; ++i, ++j) {
auto x = a[i], y = a[j];
a[i] = x + y;
a[j] = (x - y) * w;
}
w *= idw[__builtin_ctz(++k)];
}
}
auto inv = 1 / mint(n);
for (auto&& e : a) e *= inv;
}
}
vector<mint> operator*(vector<mint> l, vector<mint> r) {
if (l.empty() or r.empty()) return {};
int n = (int)l.size(), m = (int)r.size(),
sz = 1 << __lg(2 * (n + m - 1) - 1);
if (min(n, m) < 30) {
vector<mint> res(n + m - 1);
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) res[i + j] += (l[i] * r[j]);
return {begin(res), end(res)};
}
bool eq = l == r;
l.resize(sz), ntt(l, false);
if (eq)
r = l;
else
r.resize(sz), ntt(r, false);
for (int i = 0; i < sz; ++i) l[i] *= r[i];
ntt(l, true), l.resize(n + m - 1);
return l;
}
vector<mint>& operator*=(vector<mint>& l, vector<mint> r) {
if (l.empty() or r.empty()) {
l.clear();
return l;
}
int n = (int)l.size(), m = (int)r.size(),
sz = 1 << __lg(2 * (n + m - 1) - 1);
if (min(n, m) < 30) {
vector<mint> res(n + m - 1);
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) res[i + j] += (l[i] * r[j]);
l = {begin(res), end(res)};
return l;
}
bool eq = l == r;
l.resize(sz), ntt(l, false);
if (eq)
r = l;
else
r.resize(sz), ntt(r, false);
for (int i = 0; i < sz; ++i) l[i] *= r[i];
ntt(l, true), l.resize(n + m - 1);
return l;
}
} // namespace ntt
int main()
{
ios::sync_with_stdio(false); cin.tie(0);
using mint = ntt::mint;
using poly = vector<mint>;
int t; cin >> t;
while (t--) {
int n; cin >> n;
vector p(n, 0);
for (int &x : p) cin >> x;
vector<poly> v;
int pmx = 0, sz = 0;
for (int x : p) {
if (x < pmx) ++sz;
else {
if (sz) {
poly q(sz+1, 0);
q[0] = q[sz] = 1;
v.push_back(q);
}
pmx = x;
sz = 1;
}
}
poly q(sz+1, 0);
q[0] = q[sz] = 1;
v.push_back(q);
auto mul = [&] (const auto &self, int L, int R) -> poly {
if (L+1 == R) return v[L];
int mid = (L+R)/2;
return self(self, L, mid) * self(self, mid, R);
};
auto ans = mul(mul, 0, size(v));
for (auto x : ans) cout << x << ' ';
cout << '\n';
}
}