# MATPERM - Editorial

Author: Vansh Kaul
Tester: Riley Borgard
Editorialist: Istvan Nagy

Medium

# PREREQUISITES:

permutation, multi source shortest path, modular multiplicative inverse, difference array

# PROBLEM:

There is a grid of N\times M size. K of the grid cells has a cannon placed on them (A grid cell can have at most one cannon). Each cannon can throw a ball to a cell which is at a distance at most equal to the strength of the cannon (manhattan distance is to be considered). Cannon can also throw a ball at the same cell in which it’s present. The strength of all the cannons is 0 initially and it increases by 1 every second. At 0^{th} second, cannons start throwing the balls following the given rules:

• At each second exactly one cannon throws a ball to a cell which is at a distance at most X from the cannon where X is the current strength of the cannon.
• Cannon can’t throw a ball at a cell that has a ball thrown at it previously or outside the grid.
• When a ball is thrown by a cannon at the i^{th} second, number i will get printed on the cell the ball lands on.

At the end of N*M-1 seconds. The grid will be representing a permutation of numbers from 0 to N*M-1 . (Two permutations P_1 and P_2 are said to be different if at least one cell has different numbers printed on it in P_1 and P_2).

A permutation is considered to be good if it can be formed at the end of N*M - 1 seconds using the given method.

Your task is to find the fraction of total permutations that are good.

It can be proved that the fraction can be represented in the form P/Q where P and Q are co-prime integers. Print the required fraction as P\ Q^{-1} \quad mod (2500000001) where Q^{-1} denotes the modular multiplicative inverse of Q \quad mod (2500000001).

# EXPLANATION

Let’s try to construct the permutation one-by-one. On 0-th second we have to put 0 in a cell with a cannon. On 1-st second we have to put 1 in a cell within distance 1 from some cannon, but not in already used cell. Generally, on i-th second we have to put i in a cell within distance i from some cannon, but not in already used cells. It seems that we have to remember which cells are used, but let’s take a closer look. For any j < i it is true that if cell is within distance j from some cannon, it is also within distance i. Therefore, on i-th second all used cells are within distance i independent of our previous choices. Let dist2[i] be the number of cells within distance i from some cannon. Then the number of good permutations is:

\prod_{i=0}^{N*M-1} (dist2[i] - i)

Let dist[i] be the number of cells on distance exactly i from closest cannon. Then dist2[i] is just an array of prefix sums of dist. For a subtask we can find all distances using multi-source BFS, but for the full points we’ll have to do something smarter.

Observation: A cell distance from a cannon cannot exceeds N+M - 2.

Consequence: From the N+M-2-th second any cell is in reachable distance from the cannons.

The number of good permutations (using the Consequence above):

\prod_{i=0}^{N+M-3}({dist2[i]-i}) \prod_{i=N+M-2}^{N*M-1}(N*M-i).

The total permutation is simple permutation of N*M cells:

\prod_{i=0}^{N*M-1}(N*M-i)

So after the N+M-3-th second the numerator and the denominator will be equal so we can get rid off that part of the product, to get the final form of the result:

\prod_{i=0}^{N+M-3}\frac{({dist2[i]-i})}{(N*M-i)}

So now we have “just” to calculate the dist array.

First step:

Calculate the dist only the columns which contains cannons. We can do this with the multiple source dijkstra algorithm it has a O(E*log(V)) time complexity. Currently the V = N*K and E < 4*V. In the current case the distances are in a relatively small range ([0, N+M-2]) so we can replace the priority queue with buckets to speed up the algorithm. More details here. The time complexity of this step: O(N*K+M).

Second step:

Now let just concentrate on two points in the same row and two adjacent cannon column: p_1 = (x_1, y), p_2 = (x_2,y) and look at how the closest cannon distances between these two points can looks like. Let d_1 the closest cannon distance from p_1 and similarly d_2 the closest cannon distance for p_2. We can make few observation about the intermediate points distance: [d_1,...,d_2]

1. Two neighbor value difference is 0 or 1.
2. if d_1 < d_2 then the adjacent cell distance to the p_2 direction is d_1+1.
3. if d_1 = d_2 and exists at least a point between these two then the adjacent cell distance is d_1+1.

With these observation we can split the [d_1,...,d_2] array two parts:

[d_1,...,d_c], [d_{c+1},...,d_2] where the first interval values are strictly increasing by 1, and the second interval values are strictly decreasing by one. We can calculate this two interval in constant time, so generate all of this interval time complexity is: O(N*K+M)

Third step:

We can use difference array technique to calculate effectively the dist array from this monotone arrays. Create a new auxiliary array dist3 and when we have an a [d1,...,d_c] monotone increasing array, we can do the following trick:

• increase the dist3[d_1] value by one
• decrease the dist3[d_c+1] value by one

After we add all the arrays from the second step to the dist3, we can calculate the dist array :

dist[i] = \sum_{j=0}^idist3[j]

This step time complexity is:

1. add the arrays to dist3 : O(N*K+M)
2. calculate dist from dist3 : O(N+M)

# TIME COMPLEXITY:

O(N*K+M+K*log(K)) per test case

The first part overall:

\sum_t O(N*k)\leq O(N)*\sum_tk\leq 1000*O(N)

The third part overall:

\sum_t O(K*log(K))\leq 10* \sum_t O(k)

SOLUTIONS:

Setter's Solution
#pragma GCC optimize("Ofast")
#pragma GCC target("avx,avx2,fma")
#include<bits/stdc++.h>
// using namespace std;
#include <algorithm>
#include <chrono>
#include <iostream>
#include<vector>
using namespace std;
using namespace std::chrono;

#define ll long long
#define int long long
#define vll vector<ll>
#define fast  ios_base::sync_with_stdio(false);  cin.tie(NULL); cout.tie(NULL);
#define mod 2500000001

ll power(ll x, ll y, ll p) {
ll res = 1; x = x % p; if (x == 0) return 0; while (y > 0) { if (y & 1) res = (res * x) % p; y = y >> 1; x = (x * x) % p;} return res;
}
ll mmul (ll a, ll b) {
return ( (a % mod) * (b % mod) ) % mod;
}
ll mdiv(ll a, ll b) {
return ((a % mod) * power(b, mod - 2, mod) ) % mod;
}

void swap(ll &a, ll &b) {
ll temp = a;
a = b;
b = temp;
}

ll distance(pair<ll, ll> p1, pair<ll, ll> p2) {
return abs(p1.first - p2.first) + abs(p1.second - p2.second);
}
ll solve(ll n, ll m, ll k, vll x, vll y) {
vector<pair<ll, ll>> xy;
for (int i = 0; i < k; i++) {
xy.push_back({x[i] - 1, y[i] - 1});
}
sort(xy.begin(), xy.end());
vector<ll> dist(n + m + 1);
vector<pair<ll, ll>> st;
for (int yy = 0; yy < n; yy++) {
st.clear();
for (int i = 0; i < k; i++) {
if (st.size() > 0 && distance(xy[i], {xy[i].first, yy}) > distance(st.back(), {xy[i].first, yy}))
continue;
while (st.size() > 0 && distance(xy[i], {st.back().first, yy}) < distance(st.back(), {st.back().first, yy}))
st.pop_back();
st.push_back(xy[i]);
}
for (int i = 1; i < (ll)st.size(); i++) {
ll panelty1 = abs(st[i - 1].second - yy);
ll panelty2 = abs(st[i].second - yy);
//To avoid recounting
dist[abs(st[i].second - yy)]--;
dist[abs(st[i].second - yy) + 1]++;

ll d = abs(st[i].first - st[i - 1].first) + 1;

if (panelty1 > panelty2)
swap(panelty1, panelty2);

//Bringing two penalties to equal numbers
dist[panelty1]++;
dist[panelty2]--;

d -= (panelty2 - panelty1);
dist[panelty2] += 2;
dist[panelty2 + (d / 2)]--;
dist[panelty2 + (d + 1) / 2]--;
}
if (st.size() > 0) {
ll la = abs(st.back().second - yy);
ll fi = abs(st[0].second - yy);
dist[fi + 1]++;
dist[fi + st[0].first + 1]--;
dist[la]++;
dist[la + m - st.back().first]--;
}

}
// cout << dist[0] << " ";
for (int i = 1; i < n + m + 1; i++) {
dist[i] += dist[i - 1];
// cout << dist[i] << " ";
}
vector<ll> pre = dist;
for (int i = 1; i < (ll)pre.size(); i++)
pre[i] += pre[i - 1];

ll total = n * m;
ll num = 1, den = 1;
for (ll i = 0; i < min(n * m, n + m); i++) {
num = mmul(num, pre[i] - i);
den = mmul(den, total - i);
}
return mdiv(num, den);
}
signed main() {

fast;
//	auto start = high_resolution_clock::now();
#ifndef ONLINE_JUDGE
freopen("input1.txt", "r", stdin);
freopen("output1.txt", "w", stdout);
#endif
fast;
ll t;
cin >> t;
while (t--) {
ll N, M, K;
cin >> N >> M >> K;

vector<ll> x(K), y(K);

for (int i = 0; i < K; i++) {
cin >> y[i] >> x[i];
}
cout << solve(N, M, K, x, y) << '\n';
}}

Tester's Solution

#include <bits/stdc++.h>

#define ll long long
#define sz(x) ((int) (x).size())
#define all(x) (x).begin(), (x).end()
#define vi vector<int>
#define pii pair<int, int>
#define rep(i, a, b) for(int i = (a); i < (b); i++)
using namespace std;
template<typename T>
using minpq = priority_queue<T, vector<T>, greater<T>>;

template<long long m>
struct mod {
long long x;
mod() : x(0) {}
mod(long long xx) : x(xx) {
if(abs(x) >= m) x %= m;
if(x < 0) x += m;
}
mod operator+(const mod &a) const {
mod n;
n.x = x + a.x;
if(n.x >= m) n.x -= m;
return n;
}
mod operator-(const mod &a) const {
mod n;
n.x = x - a.x;
if(n.x < 0) n.x += m;
return n;
}
mod operator*(const mod &a) const {
return x * a.x;
}
mod operator+=(const mod &a) {
x += a.x;
if(x >= m) x -= m;
return *this;
}
mod operator-=(const mod &a) {
x -= a.x;
if(x < 0) x += m;
return *this;
}
mod operator*=(const mod &a) {
x = (x * a.x) % m;
return *this;
}
mod pow(long long b) const {
mod ans = 1;
mod a = *this;
while(b > 0) {
if(b & 1) ans *= a;
a *= a;
b /= 2;
}
return ans;
}
mod inv() const {
return pow(m - 2);
}
mod operator/(const mod &a) const {
return (*this) * a.inv();
}
mod operator/=(const mod &a) {
return (*this) *= a.inv();
}
bool operator==(const mod &o) const {
return x == o.x;
}
bool operator!=(const mod &o) const {
return x != o.x;
}
long long operator()() const {
return x;
}
};

template<ll m>
istream &operator>>(istream &is, mod<m> &n) {
long long x;
is >> x;
n = x;
return is;
}

template<ll m>
ostream &operator<<(ostream &os, const mod<m> &n) {
return os << n();
}

const ll M = 2500000001LL;
using base = mod<M>;

void solve() {
int n, m, k;
cin >> n >> m >> k;
vector<pii> vp;
rep(i, 0, k) {
int x, y;
cin >> x >> y;
x--; y--;
vp.push_back({x, y});
}
sort(all(vp));
vi st;
auto dist = [&](int x1, int y1, int x2, int y2) {
return abs(x1 - x2) + abs(y1 - y2);
};
vector<ll> d(n + m + 5, 0);
rep(c, 0, m) {
st.clear();
rep(i, 0, k) {
int x, y; tie(x, y) = vp[i];
if(st.empty()) {
st.push_back(i);
}else {
int x2, y2; tie(x2, y2) = vp[st.back()];
if(dist(x2, y2, x, c) > abs(y - c)) {
while(!st.empty()) {
tie(x2, y2) = vp[st.back()];
if(dist(x, y, x2, c) < abs(y2 - c)) st.pop_back();
else break;
}
st.push_back(i);
}
}
}
int L = 0, R = n - 1;
int s = sz(st);
rep(i, 0, s) {
int x, y; tie(x, y) = vp[st[i]];
int t = abs(y - c);
// left
d[t]++;
d[t + x - L + 1]--;
// right
d[t + 1]++;
if(i == s - 1) {
d[t + R - x + 1]--;
}else {
int T = (vp[st[i + 1]].first - x + abs(vp[st[i + 1]].second - c) + t + 1) / 2;
T = min(T, R + 1 - x + t);
T = max(T, t + 1);
d[T]--;
L = x + T - t;
}
}
}
rep(i, 1, n + m + 5) d[i] += d[i - 1];
rep(i, 1, n + m + 5) d[i] += d[i - 1];
base ans = 1;
ll mx = min(1LL * n * m, (ll) n + m + 5);
base den = 1;
rep(i, 0, mx) {
ans *= (d[i] - i);
den *= (1LL * n * m - i);
}
ans /= den;
cout << ans << '\n';
}

int main() {
ios::sync_with_stdio(false);
cin.tie(0);
int te;
cin >> te;
while(te--) solve();
}


If you have other approaches or solutions, let’s discuss in comments.