PROBLEM LINK:
Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4
Author: iceknight1093
Tester: sushil2006
Editorialist: iceknight1093
DIFFICULTY:
Medium
PREREQUISITES:
Dynamic programming
PROBLEM:
For a permutation P, define f(P) to be the minimum number of operations of this type needed to sort it:
- Choose L, R such that P[L\ldots R] is not sorted.
- Then, either move P_L to the end of P, or P_R to the start of P.
Given N, K, and X, count the number of permutations of length N for which f(P) = K and P_1 = X.
EXPLANATION:
First, let’s understand how to compute f(P).
To do this, it’s useful to look at the elements that aren’t moved by any operation.
Since each operation moves an element to either the front or the end, the elements that haven’t been moved by an operation yet will form a contiguous subarray.
This means, in the end when the permutation is sorted, the not-moved elements will be of the form [x, x+1, \ldots, y] for some x and y.
Note that this means that the elements x, x+1, x+2, \ldots, y formed a subsequence in the initial permutation.
Let’s look at the longest subsequence of this form, i.e. consisting of consecutive integers placed in ascending order.
We can then always sort P while ensuring that these elements are not operated on, while everything else is operated on.
How?
We have the subsequence x, x+1, \ldots, y
Consider the value y+1.
Since this subsequence was maximal, y+1 cannot appear after y - it must be before y.
This means the suffix starting with the element y+1 is not sorted, since it contains y after y+1.
Choose this suffix, and then place y+1 at the end.Note that after this is done, the suffixes starting at each of the values y+2, y+3, \ldots, N will all not be sorted either (since y+1 will appear after them), so by choosing each of these suffixes in order they can all be placed at the end.
Similarly, x-1 must appear after x, so using the prefix ending at it we can place it at the beginning; following which all of x-2, x-3, \ldots, 2, 1 can be placed at the beginning in order.
This sorts P.
So, quite simply, if the longest increasing subsequence in P consisting of consecutive values has length L, we have f(P) = N-L.
Using this characterization, we can compute the answer.
Let’s ignore the condition of P_1 = X for now.
We want f(P) = K, which means the count of permutations whose longest increasing consecutive subsequence has length N-K.
This can be computed using dynamic programming.
Define dp(i, j, l) to be the number of permutations of [1, i] for which P_j = i, the increasing consecutive subsequence ending at value i has length l, and every increasing consecutive subsequence has length \leq N-K.
Let’s try to place value i+1 into this permutation.
- If we place i+1 before value i, i.e. at an index \lt j, the old subsequence breaks and a new one starts at i+1.
So, for each 0 \leq x \lt j, we add dp(i, j, l) to dp(i+1, x+1, 1). - If we place i+1 after i, the subsequence extends to include it.
For each x \geq j, we add dp(i, j, l) to dp(i+1, x+1, l+1).
Note that this is only allowed if l \lt N-K since we cannot exceed N-K.
In the end, the sum of dp(N, i, l) for all (i, l) gives the number of permutations whose longest increasing consecutive subsequence has length \leq N-K.
To obtain the count of it being exactly N-K, subtract out the count of permutations for which the answer is \leq N-K-1, which can be obtained by running basically the same DP again.
Our DP has \mathcal{O}(N^3) states and we do \mathcal{O}(N) updates from each, for \mathcal{O}(N^4) overall.
It’s possible to improve this to \mathcal{O}(N^3) using prefix sums, but that optimization wasn’t necessary here.
Now, we bring in P_1 = X.
First, let’s solve the case of X = 1, i.e. the smallest element is always the leftmost.
This is not too hard: you can either modify the DP transitions slightly (to not allow any element \gt 1 to be placed at the very beginning), or simply note that counting configurations with P_1 = 1 is equivalent to counting configurations with P_N = N (consider what happens if you compute the DP by iterating in reverse order and store the position of the minimum rather than the maximum)
So, simply summing up dp(N, N, l) across all l \leq N-K with our original DP definition will give us the required value (of course, don’t forget to subtract out the count for answer \leq N-K-1 afterwards.)
Next, what about X \gt 1?
Here, observe that since the very first element is X, the “chain” is broken at X-1.
That is, the sets of elements [1, X-1] and [X, N] are essentially completely independent of each other in terms of forming increasing consecutive subsequences.
This allows us to do the following:
- Fix the set of positions that contain elements \leq X-1.
This can be done in \binom{N-1}{X-1} ways. - Then, use the existing DP table to compute the number of valid configurations of these elements with answer \leq N-K, which is just the sum of dp(X-1, i, l) across all (i, l).
- This leaves the range [X, N].
However, this is equivalent to solving for [1, N-X+1] with the first element being 1, and we already saw how to do that above.
So, after computing all the dp(i, j, l) values for 1 \leq i \leq N, 1 \leq j \leq i, 1 \leq l \leq \min(i, N-K) just once, we’ve done some simple combinatorics to combine them into the overall answer.
TIME COMPLEXITY:
\mathcal{O}(N^4) or \mathcal{O}(N^3) per testcase.
CODE:
Editorialist's code (C++)
// #include <bits/allocator.h>
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include "bits/stdc++.h"
using namespace std;
using ll = long long int;
mt19937_64 RNG(chrono::high_resolution_clock::now().time_since_epoch().count());
/**
* Integers modulo p, where p is a prime
* Source: Aeren (modified from tourist?)
* Modmul for 64-bit mod from kactl:ModMulLL
* Works with p < 7.2e18 with x87 80-bit long double, and p < 2^52 ~ 4.5e12 with 64-bit
*/
template<typename T>
struct Z_p{
using Type = typename decay<decltype(T::value)>::type;
static vector<Type> MOD_INV;
constexpr Z_p(): value(){ }
template<typename U> Z_p(const U &x){ value = normalize(x); }
template<typename U> static Type normalize(const U &x){
Type v;
if(-mod() <= x && x < mod()) v = static_cast<Type>(x);
else v = static_cast<Type>(x % mod());
if(v < 0) v += mod();
return v;
}
const Type& operator()() const{ return value; }
template<typename U> explicit operator U() const{ return static_cast<U>(value); }
constexpr static Type mod(){ return T::value; }
Z_p &operator+=(const Z_p &otr){ if((value += otr.value) >= mod()) value -= mod(); return *this; }
Z_p &operator-=(const Z_p &otr){ if((value -= otr.value) < 0) value += mod(); return *this; }
template<typename U> Z_p &operator+=(const U &otr){ return *this += Z_p(otr); }
template<typename U> Z_p &operator-=(const U &otr){ return *this -= Z_p(otr); }
Z_p &operator++(){ return *this += 1; }
Z_p &operator--(){ return *this -= 1; }
Z_p operator++(int){ Z_p result(*this); *this += 1; return result; }
Z_p operator--(int){ Z_p result(*this); *this -= 1; return result; }
Z_p operator-() const{ return Z_p(-value); }
template<typename U = T>
typename enable_if<is_same<typename Z_p<U>::Type, int>::value, Z_p>::type &operator*=(const Z_p& rhs){
#ifdef _WIN32
uint64_t x = static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value);
uint32_t xh = static_cast<uint32_t>(x >> 32), xl = static_cast<uint32_t>(x), d, m;
asm(
"divl %4; \n\t"
: "=a" (d), "=d" (m)
: "d" (xh), "a" (xl), "r" (mod())
);
value = m;
#else
value = normalize(static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value));
#endif
return *this;
}
template<typename U = T>
typename enable_if<is_same<typename Z_p<U>::Type, int64_t>::value, Z_p>::type &operator*=(const Z_p &rhs){
uint64_t ret = static_cast<uint64_t>(value) * static_cast<uint64_t>(rhs.value) - static_cast<uint64_t>(mod()) * static_cast<uint64_t>(1.L / static_cast<uint64_t>(mod()) * static_cast<uint64_t>(value) * static_cast<uint64_t>(rhs.value));
value = normalize(static_cast<int64_t>(ret + static_cast<uint64_t>(mod()) * (ret < 0) - static_cast<uint64_t>(mod()) * (ret >= static_cast<uint64_t>(mod()))));
return *this;
}
template<typename U = T>
typename enable_if<!is_integral<typename Z_p<U>::Type>::value, Z_p>::type &operator*=(const Z_p &rhs){
value = normalize(value * rhs.value);
return *this;
}
template<typename U>
Z_p &operator^=(U e){
if(e < 0) *this = 1 / *this, e = -e;
Z_p res = 1;
for(; e; *this *= *this, e >>= 1) if(e & 1) res *= *this;
return *this = res;
}
template<typename U>
Z_p operator^(U e) const{
return Z_p(*this) ^= e;
}
Z_p &operator/=(const Z_p &otr){
Type a = otr.value, m = mod(), u = 0, v = 1;
if(a < (int)MOD_INV.size()) return *this *= MOD_INV[a];
while(a){
Type t = m / a;
m -= t * a; swap(a, m);
u -= t * v; swap(u, v);
}
assert(m == 1);
return *this *= u;
}
template<typename U> friend const Z_p<U> &abs(const Z_p<U> &v){ return v; }
Type value;
};
template<typename T> bool operator==(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value == rhs.value; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(const Z_p<T>& lhs, U rhs){ return lhs == Z_p<T>(rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) == rhs; }
template<typename T> bool operator!=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return !(lhs == rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(const Z_p<T> &lhs, U rhs){ return !(lhs == rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(U lhs, const Z_p<T> &rhs){ return !(lhs == rhs); }
template<typename T> bool operator<(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value < rhs.value; }
template<typename T> bool operator>(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value > rhs.value; }
template<typename T> bool operator<=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value <= rhs.value; }
template<typename T> bool operator>=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value >= rhs.value; }
template<typename T> Z_p<T> operator+(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(const Z_p<T> &lhs, U rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T> Z_p<T> operator-(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T> Z_p<T> operator*(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T> Z_p<T> operator/(const Z_p<T> &lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(const Z_p<T>& lhs, U rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(U lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T> istream &operator>>(istream &in, Z_p<T> &number){
typename common_type<typename Z_p<T>::Type, int64_t>::type x;
in >> x;
number.value = Z_p<T>::normalize(x);
return in;
}
template<typename T> ostream &operator<<(ostream &out, const Z_p<T> &number){ return out << number(); }
// constexpr int mod = 1e9 + 7; // 1000000007
constexpr int mod = (119 << 23) + 1; // 998244353
// constexpr int mod = 1e9 + 9; // 1000000009
using Zp = Z_p<integral_constant<decay<decltype(mod)>::type, mod>>;
const int N = 105;
Zp dp[N][N][N], C[N][N];
// dp[i][j][k] = permutations of 1...i with p[j] = i, last segment length = k
int main()
{
ios::sync_with_stdio(false); cin.tie(0);
for (int n = 0; n < N; ++n) {
C[n][0] = 1;
for (int r = 1; r <= n; ++r)
C[n][r] = C[n-1][r] + C[n-1][r-1];
}
auto calc = [&] (int n, int k, int x) {
// length n, max chain <= k, p[1] = x
if (k <= 0) return Zp(0);
// Reset
for (int i = 1; i <= n; ++i) for (int j = 1; j <= i; ++j) for (int l = 1; l <= i and l <= k; ++l)
dp[i][j][l] = 0;
dp[1][1][1] = 1;
for (int i = 1; i < n; ++i) {
for (int j = 1; j <= i; ++j) {
for (int l = 1; l <= min(i, k); ++l) {
// Place i+1
// If placed before j, l becomes 1
// If placed after j, l increases by 1
for (int p = 0; p <= i; ++p) {
if (p < j) dp[i+1][p+1][1] += dp[i][j][l];
else {
if (l+1 <= k) dp[i+1][p+1][l+1] += dp[i][j][l];
}
}
}
}
}
// [1, x-1]
Zp ans = C[n-1][x-1], val = 0;
for (int i = 1; i <= x-1; ++i)
for (int j = 1; j <= min(x-1, k); ++j)
val += dp[x-1][i][j];
if (x == 1) val = 1;
ans *= val;
// [x, n]
// length n-x+1, but smallest element first
// equivalent to largest element last, can just use existing dp
int len = n-x+1;
val = 0;
for (int i = 1; i <= min(len, k); ++i)
val += dp[len][len][i];
ans *= val;
return ans;
};
int t; cin >> t;
while (t--) {
int n, k, x; cin >> n >> k >> x;
k = n-k;
cout << calc(n, k, x) - calc(n, k-1, x) << '\n';
}
}