## Sunday, May 28, 2017

### Codeforces Round #415 (Div. 2) E. Find a car [or Div. 1 C]

Problem Statement:
Codeforces Round #415 (Div. 2) E. Find a car

Paraphrase:
An integer matrix M is defined as:
M[i][j] = the minimum integer x >= 1 such that:
- M[i][k] != x for all k < j
- M[k][j] != x for all k < i

Example:
1 2 3 4
2 1 4 3
3 4 1 2
4 3 2 1

Goal: Given a rectangular block in M {(x1, y1), (x2, y2)} and an integer K, compute the sum of elements in that rectangle, but exclude elements whose value is larger than K.

Solution:

Let's construct the matrices with sizes 1, 2, 4, 8, ...
Size 1:
1

Size 2:
1 2
2 1

Size 4:
1 2 3 4
2 1 4 3
3 4 1 2
4 3 2 1

Size 8:
1 2 3 4 5 6 7 8
2 1 4 3 6 5 8 7
3 4 1 2 7 8 5 6
4 3 2 1 8 7 6 5
5 6 7 8 1 2 3 4
6 5 8 7 2 1 4 3
7 8 5 6 3 4 1 2
8 7 6 5 4 3 2 1

Notice the pattern yet?

Yes. So given a matrix M[k] with size 2^k, we can construct M[k+1] by appending M[k] and N[k] as follows:
M[k] N[k]
N[k] M[k]
where N[k] is M[k] with all elements incremented by 2^k. (Which can be proven inductively, an exercise for the reader :) )

Now how is this useful, you ask. This is useful because now you can perform a technique similar to quad tree traversal, which if done correctly, will allow us to answer the query in O(log N) where N is the size of the matrix.

Notice that each row of M[k+1] consists of values [1, 2, 3, ..., 2 ^ {k+1}] (in some sudoku-like permutation).
Also, M[k] only consists of values in [1, ..., 2^k] while N[k] only of [2^k + 1, ...,  2^{k+1}].
Hence when you split M[k+1] into 4 equal sized quadrants, you get 4 similar construct, they just differ in ranges of values they contain.
You can also check that when you split N[k] into 4 similar sectors, you would also see this pattern holds true.

From this point, you may already know what to do next:
1. Start from the whole matrix (you can define it as [1, ..., 2^31] for this problem)
2. If our current space is still not fully inside the target rectangle, split up into 4 quadrants and then recursively proceed (don't forget to carry the information on the range of values our space currently contains). Each quadrant will be our next "space" to check.
3. If our current space does not even intersect with the target rectangle, throw away this quadrant.
4. If our current space is fully contained in our target rectangle, you can quite easily compute the sum of elements in our space that does not exceed K.

The approach above is correct, but there is one problem. Consider a query such as {(1, 1), (1, 1e9)} for any K. You will end up summing up each of the cells that falls in that rectangle 1e9 times!
One easy fix to this, is to relax our requirement of having our space to be "fully inside" the target rectangle. Instead, as long as one of the side of the space is fully contained in the target rectangle, we can already compute the required sum (Exercise to the reader! Or you can figure it out in the code below). With that, we are sure to have an O(log N) query time.

Code:

#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;
constexpr int64_t MOD = (int64_t) 1e9 + 7LL;
int64_t x1, y1;
int64_t x2, y2;
return (x1 <= other.x1 && other.x2 <= x2) ||
(y1 <= other.y1 && other.y2 <= y2);
}
return !(x2 < other.x1 || x1 > other.x2
|| y2 < other.y1 || y1 > other.y2);
}
};

int64_t contribute(int64_t L, int64_t R, int64_t k, int64_t n) {
if (k < L) return 0;
int64_t A = L + min(k, R);
int64_t B = min(k,R) - L + 1;
if (A % 2LL) B /= 2LL;
else A /= 2LL;
A %= MOD;
B %= MOD;
return (((A * B) % MOD) * n) % MOD;
}

int64_t dfs(quad space, int64_t L, int64_t R, quad target, int64_t k) {
if (!target.intersects(space)) {
return 0;
}
if (target.contains(space)) {
int64_t n = min(
min(target.x2, space.x2) - max(target.x1, space.x1) + 1,
min(target.y2, space.y2) - max(target.y1, space.y1) + 1
);
return contribute(L, R, k, n % MOD);
}
int64_t ans = 0;
int64_t M = (L + R) / 2LL;
int64_t m = (R - L) / 2LL;
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
next_space.x1 = space.x1 + (i ? m + 1: 0);
next_space.x2 = (i ? space.x2 : space.x1 + m);
next_space.y1 = space.y1 + (j ? m + 1 : 0);
next_space.y2 = (j ? space.y2 : space.y1 + m);
if (i ^ j) {
ans += dfs(next_space, M+1, R, target, k);
ans %= MOD;
} else {
ans += dfs(next_space, L, M, target, k);
ans %= MOD;
}
}
}
return ans;
}

void solve(int x1, int x2, int y1, int y2, int64_t k) {
int64_t L = 1;
int64_t R = 1LL << 31;
space.x1 = space.y1 = L;
space.x2 = space.y2 = R;
target.x1 = x1;
target.x2 = x2;
target.y1 = y1;
target.y2 = y2;
cout << dfs(space, L, R, target, k) << endl;
}

int main(){
int q;
scanf("%d",&q);
for (int i = 0; i < q; ++i) {
int64_t x1,x2,y1,y2,k;
cin >> x1 >> y1 >> x2 >> y2 >> k;
solve(x1,x2,y1,y2,k);
}
return 0;
}