Problem Link:
Difficulty:
Hard
Pre-requisites:
Math, Combinatorics
Problem:
You are given N. Count how many NxN matrices are there, satisfying the following conditions:
(a) M[x][x] = 0
(b) M[x][y] = M[y][x] > 0
(c) M[x][y] <= max(M[x][z], M[z][y])
(d) M[x][y] (x!= y) is taken from the set “[N-2]”:= {1, 2, …, N-2}, and further every element from the set is chosen for some M[x][y] (x != y).
The constraint is: 3 <= N <= 10^7.
Quick explanation:
It seems that many users find the solution writing brute-force for small values of N, searching for the corresponding sequence at OEIS, finding it here and using hint about fractional iteration of e^x - 1, which leads to using higher-order Bell numbers as explained here. This approach leads to precalculation where 4 * N + O(1) “heavy” operations A * B % MOD is used and it can get AC only after proper optimizations.
The problem-setter was aware of presence on this sequence in OEIS while he creates the problem, but we decided that information provided there is not enough to solve the problem. Well we were wrong 
One of the purposes of this problem is to emphasize importance of minimizing the number of these heavy operations. The problem-setter solution described below use completely different approach and leads to a solution with only 2 * N + O(sqrt(N)) heavy operations in precalc step, while after this, the answer for each N can be found in 3 additional heavy operations. This leads to 2 * N + O(sqrt(N)) + 3 * T heavy operations in all and allows to pass the TL using only 0.300 seconds from 0.666. Which is quite far from strict. Also we have a solution that precalc all answers up to N using 3 * N + O(1) heavy operations and one more that use only 2 * N + O(1) heavy operations for such precalc, but due to cache-misses even it can’t pass the TL. The final formula used in all 3 solutions is:
N! * (N-1)! / 2^(N-1) * (N/2 - 2/3 - H(N-1)/3),
where H(K) = 1 + 1/2 + 1/3 + … + 1/K - the K-th Harmonic number.
Explanation:
The matrix you are supposed to find is very similar to a kind of distance function (condition (c) clearly implies triangle inequality). In fact, such a metric space (one that satisfies conditions (a) - (c)) is called an “ultrametric space”.
The spectrum of an ultrametric space X, denoted by Sp(X), is defined as the set of non-zero values that the distance can take. Note that, given a particular spectrum set S, we can transform it into an equivalent spectrum set {1, 2, …, |S|}, just by keeping relative order. The question thus asks to find, given N, how many ultrametric spaces are there, having spectrum of size N-2.
We call the “diameter” of the space, denoted as diam(X), the largest distance between any two points in the space. Now, let us derive some useful properties of the ultrametric space. Also, we further denote M[x][y] as “d(x, y)” or “dist(x, y)”.
Intuitively, we would like to fill in values for dist(x, y) in a particular order. We could either start from 1 to diam(X), or we can go from diam(X) to 1. We end up choosing to fill in values from diam(X) to 1.
Consider the relation R(d), (d > 0) on the set of points, where R(d)(x, y) iff dist(x, y) < d. This relation R(d) is an equivalence relation as can be verified:
Reflexivity: dist(x, x) = 0 < d => R(d)(x, x).
Symmetry: R(d)(x, y) => dist(x, y) < d => dist(y, x) < d => R(d)(y, x).
Transitivity: Given R(d)(x, y) and R(d)(y, z), we have dist(x, y) < d and dist(y, z) < d.
Now, dist(x, z) <= max(dist(x, y), dist(y, z)) < d => R(d)(x, z).
Hence, R(d) is an equivalence relation, for any d > 0.
This further gives us, that if we were to fill in values for diam(X) in dist(x, y), we would have to partition the set of points X (this partition would define our equivalence classes for R(diam(X))), and then we’d get that all d(x, y) = diam(X) ACROSS such partitions!
We also get from the above, that |Sp(X)| <= |X| - 1. (You can prove this “directly” by considering relations R(d) : dist(x, y) < d, and letting d go from 1 to |Sp(X)|+1. Each time, R(i) merges atleast two equivalence classes of R(d-1).)
We also provide an inductive proof whose logic will be used in further parts of the Editorial.
Clearly, for |X| = 1, this is true, since Sp(X) = phi.
Now, consider R(diam(X)), which partitions X into X1, X2, …, Xk, each of which is an ultrametric space. We thus have,
Sp(X) = Sp(X1) U Sp(X2) U … U Sp(Xk) U {diam(X)}.
Now, using the inequality, we get
|Sp(X)| <= sum(|Sp(Xi)|) + 1 (this is because |A U B| <= |A| + |B|) (ineq 1)
<= sum(|Xi|-1) + 1 (this is by induction) (ineq 2)
= |X| - k + 1 (Since Xi are disjoint, and Union Xi = X)
<= |X| - 1 (Since k >= 2). (ineq 3)
The important thing to keep in mind in the above is:
|Sp(X)| = |X| - 1 iff at each stage, we break up our set into exactly two partitions, and further, the spectra of each partition are all disjoint.
Let f(N, K) = number of ultrametric spaces X, where |X| = N, and |Sp(X)| = k.
Our ultimate goal is to find f(N, N-2).
Now, suppose we have our space X, and we need |Sp(X)| = |X|-2. This can happen either due to
(a) There are three equivalence classes for R(N-2), and all three subspaces have disjoint spectra (this is due to a loss of “1” from ineq. 3)
(b) There are two equivalence classes for R(N-2), and exactly one value is common to both subspace (this is due to a loss of 1 from ineq. 1)
(c) There are two equivalence classes, both of them have disjoint spectra, but one of them has a spectra of size |Xi|-2 (this is due to a loss of 1 from ineq. 2).
We denote the answer to case (a) as A1, to case (b) as A2, to case (c) as A3.
We can now deduce recurrences for A1, A2, A3, and one important term that would arise in all of them is f(N, N-1) : [in A1, our spectra all have size |Xi|-1. in A2, both spectra have such size. in A3, one spectrum has |Xi|-1].
f(N, N-1) can be thought of as: first partitioning the set X into two sets of size k and N-k, then dividing the set of Spectra from the union {1, 2, …, N-2} into parts of sizes k-1, and N-k-1 and then, by re-numbering the spectra so that they are consecutive integers (1 to k-1, or 1 to N-k-1 respectively), we get the term corresponding to f(k, k-1) and f(N-k, N-k-1).
That is, f(N, N-1) = sum[Comb(N-1, k-1) * Comb(N-2, k-1) * f(k, k-1) * f(N-k, N-k-1), 1 <= k <= N-1]. The first term corresponds to “the partition of size k containing a fixed element, say 1”. The second term corresponds to dividing the set of distances {1, 2, …, N-2} into the part consisting of size k-1, the rest going to the other part, and f(k, k-1) and f(N-k, N-k-1) correspond to the number of ways of getting ultrametric spaces satisfying the required conditions.
Using induction, one can get that f(N, N-1) = N! * (N-1)! / 2^(N-1).
Apart from induction, we use the following counting argument to get at the above.
This can be modified slightly to prove/derive values for A1, A2 later.
First fix two permutations, a permutation of the points from 1 to N (p-permutation), and a permutation of “distances between consecutive points” from 1 to N-1 (d-permutation).
That is, we visualize our pair (p, d) as follows:
p(1) ------ p(2) ------ p(3) ------ ... -------- p(N)
d(1) d(2) d(3) ... d(N-1)
Lets call this pair (p, d). We now map this permutation-pair to an ultrametric space as follows:
Look at where N-1 is in the d-permutation. Give the distance N-1 to all pairs of points going across this ‘edge’ in the d-permutation. This divides our set of points into those to the left, and those to the right, and further divides our set of distances as those to the left and those to the right of distance N-1. Now recurse.
In other words, given this (p, d) pair, one can get the ultrametric space as: dist(x, y) is got by finding the largest value of the d-permutation that lies between the positions of x and y in the p-permutation.
Now further, we get that we will be double counting because, at each point, we are dividing our set of points and distances into two parts, but the manner in which we do this does not depend on whether we have (left, right) or (right, left). Basically, for each value of d(i), we can swap the left and right parts of the permutations p and d, and still get the same ultrametric space.
To illustrate this, consider N = 4, and the (p, d) pair as
1 ----- 2 ----- 3 ----- 4
2 1 3
The space that its mapped to will be the same as (by swapping about “d=3”)
4 ----- 3 ----- 2 ----- 1
3 1 2
will be same as (by swapping about “d=1”)
4 ----- 2 ----- 3 ----- 1
3 1 2
etc.
Indeed, we will be double counting for each choice of d, whether to have it as left-right or right-left permutations. This gives us our 2^(N-1) in the denominator. Thus, f(N, N-1) = N! * (N-1)! / 2^(N-1).
Now, let us get back to calculating A1.
We first divide our set of points into three parts X1, X2, X3, having size n1, n2, n3, and then divide our distances {1, 2, …, N-3} into three parts of size n1-1, n2-1, n3-1, and then assign ultrametric spaces on those subspaces.
This gives us, (similar to f(N, N-1)):
A1 = sum[Comb(N-1, n1-1) * Comb(N-n1-1, n2-1) * Comb(N-3, n1-1) * Comb(N-n1-2, n2-1) * f(n1, n1-1) * f(n2, n2-1) * f(n3, n3-1), n1+n2+n3 = N, n1, n2, n3>=1].
For this, one can compute this directly to arrive at the answer using rather robust calculation. The better way is to modify the above mapping with (p,d)-pair.
We now take our same two permutations p and d. Replace N-1 with N-2 in d. We thus get three “parts” for the value N-2, and this corresponds to a case of A1. We now have to re-analyze our double-counting.
For distances < N-2, they will partition it only into two parts, thus giving us a double-count of 2 each. If we look at parts divided by N-2, say
part1 ----- part2 ----- part3
N-2 N-2
we can permute part1, part2, part3 arbitrarily among themselves, giving us a double-count of 3!. Finally, since the two d-permutations
... N-1 ... N-2 ...
... N-2 ... N-1 ...
yield the same distances, after replacing N-1 with N-2, we get a double-count of 2 from here.
Putting it all together gives,
A1 = N! * (N-1)! / 2^(N-3) / 3! / 2 = N! * (N-1)! / 2^(N-1) / 3.
For A2, we divide our set of points into sets of sizes k and N-k, choose k-1 distances from {1, 2, …, N-3} for the first part, choose on these k-1 distance as common distance of both parts and assign it and all other distances to the second part, and multiply the answer by the number of ultrametric subspaces for each part.
This gives,
A2 = sum[Comb(N-1, k-1) * Comb(N-3, k-1) * (k-1) * f(k, k-1) * f(N-k, N-k-1), 2 <= k <= N-2].
This simplifies to N! * (N-1)! / 2^(N-1) * (N-3)/6, which can be also derived from the (p,d)-pair approach as follows.
Let (p,d)-pair be fixed. Choose some k from 1 to N-3 and replace the value N-1 with k, in the d-permutation. Then get your ultrametric space using the given mapping.
For double-counting, notice that in case A2, we need that the value N-2 should be between k and N-1 in the original d-permutation. Further, we have that
... k ... N-2 ... N-1 ...
gives the same space as
... N-1 ... N-2 ... k ...
Thus, we take just one out of the 3! possible ways of ordering k, N-2, N-1 in the d-permutation. This gives us a double-count of 3!. Further, each time we divide our set into two parts, and this can be double-counted twice in the left-right manner.
Hence, we get A2 = (N-3) * [N! * (N-1)!/3!/2^(N-1)].
Finally, we consider A3. Here, we pick one set of size k, and give it a spectrum of size k-2, and give the remaining N-k elements a disjoint spectrum of size N-k-1.
Thus,
A3 = sum[Comb(N, k) * Comb(N-3, k-2) * f(k, k-2) * f(N-k, N-k-1), 3 <= k <= N-1].
Unfortunately, while we can still use our (p, d) permutation-pair and replace N-1 with k to give us a mapped ultrametric space, finding out how much we double-count this is hard.
Combining all three counted numbers together we get the final recurrence
f(N,N-2) = A1 + A2 + A3 = N! * (N-1)! / 2^N * (N-1)/3 + N! * (N-3)! / 2^(N-1) * Sum[2^k / k! / (k-2)! * f(k,k-2) : 3 <= k <= N-1].
Putting g(N) = f(N,N-2) / N! / (N-1)! * 2^(N-1) we get after straightforward calculations
g(N) = (N-1)/6 + 2/(N-1)/(N-2) * Sum[(k-1) * g(k) : 3 <= k <= N-1].
Using standard approach we can get rid of sum eliminating it using expression for g(N+1) and g(N):
g(N+1) = N/6 + 2/N/(N-1) * Sum[(k-1) * g(k) : 3 <= k <= N].
g(N) = (N-1)/6 + 2/(N-1)/(N-2) * Sum[(k-1) * g(k) : 3 <= k <= N-1].
Multiply the first equation by N, and the second by (N-2) and take difference, which will lead to:
N * g(N+1) - (N-2) * g(N) = N * N/6 - (N-1) * (N-2)/6 + 2/(N-1) * ((N-1) * g(N))
This leads to
g(N+1) = 1/2 - 1/3/N + g(N),
with initial data g(3) = 1/3.
This clearly leads to
g(N) = N/2 - 2/3 - H(N-1)/3,
where H(K) = 1 + 1/2 + 1/3 + … + 1/K - the K-th Harmonic number.
Finally we get
F(N) := f(N, N-2) = N! * (N-1)! / 2^(N-1) * (N/2 - 2/3 - H(N-1)/3).
Computing the value of F(N)
Given this, let W(N) = N! * (2 + H(N)). Then,
W(0) = 2,
W(N) = N * W(N-1) + (N-1)!
Also, F(N) = (1/2)^N * N! * (N! - 2 * (W(N-1)/3)).
Thus, given N, we will need to have pw2[N] := 1/2^N, fact[N] := N! and W(N) with us.
The easiest way is to precalculate pw2[], fact[] and W[] arrays completely, using 3 * maxN + O(1) heavy operations and then use 2 heavy operations to calculate each F(N). And this way should easily pass the TL.
However, we can exploit the fact that T is much smaller than maxN and do better. So as above we precalculate values of fact[] and W[] upto maxN, but values of pw2[n] only for n upto sqN := ceiling(sqrt(maxN)). Also we need values pw22[n] = 1/2^(sqN * n) for n upto sqN.
Then for each N upto maxN we can represent it as N = sqN * u + v, with u = N / sqN and v = N % sqN and calculate 1/2^N as 1/2^(u * sqN + v) = (1/2^sqN)^u * 1/2^v = pw22[u] * pw2[v]. This leads only to 2 * N + O(sqrt(N)) heavy operations in precalc and 3 heavy operations for each F(N) as already mentioned in Quick explanation.
You can find the second solution mentioned in Quick explanation here. Though it uses only 2 * N + O(1) heavy operations in all it gets TLE. It is a good exercise to understand why 
The third mentions solution is here. It uses 3 * N + O(1) heavy operations in all but almost two times faster than the previous TL solution. It is based on the above formula for F(N) but use completely different idea for computing it. If you ask, problem-setter could explain it in detail.
Setter’s Solution:
Can be found here
Tester’s Solution:
Will be updated soon.
