#include <bits/stdc++.h>
using namespace std;
using Tableau = vector<vector<int>>;

// Calculates the Robinson–Schensted–Knuth mapping in O(n^2)
pair<Tableau, Tableau> rskMapping(const vector<int>& pi) {
	Tableau p, q;
	for (int ind = 0; ind < pi.size(); ++ind) {
		int cur = pi[ind];
		bool found = 0;
		for (int i = 0; i < p.size() && !found; ++i) {
			int j = upper_bound(p[i].begin(), p[i].end(), cur) - p[i].begin();
			if (j == p[i].size()) {
				p[i].push_back(cur);
				q[i].push_back(ind);
				found = 1;
			} else {
				swap(p[i][j], cur);
			}
		}
		if (! found) {
			p.emplace_back(vector<int>(1, cur));
			q.emplace_back(vector<int>(1, ind));
		}
	}
	return {p, q};
}

// Calculates the inverse of the Robinson–Schensted–Knuth mapping in O(n^2)
vector<int> rskInverse(int n, Tableau p, Tableau q) {
	vector<int> pi(n);
	for (int ind = n-1; ind >= 0; --ind) {
		int i = 0, j = 0;
		for (; i < q.size(); ++i) {
			for (j = 0; j < q[i].size() && q[i][j] != ind; ++j) {}
			if (j < q[i].size()) break;
		}

		int cur = p[i][j];
		p[i].pop_back();
		q[i].pop_back();

		for (--i; i >= 0; --i) {
			j = upper_bound(p[i].begin(), p[i].end(), cur) - p[i].begin();
			swap(p[i][j - 1], cur);
		}
		pi[ind] = cur;
	}
	return pi;
}

#include <bits/stdc++.h>
using namespace std;
using Tableau = vector<vector<int>>;

// Calculates the first k rows of P in the Robinson-Schensted-Knuth mapping in O(nk)
Tableau partialRSK(const vector<int>& pi, int k) {
	Tableau p(k);
	for (int ind = 0; ind < pi.size(); ++ind) {
		int cur = pi[ind];
		for (int i = 0; i < k; ++i) {
			int j = upper_bound(p[i].begin(), p[i].end(), cur) - p[i].begin();
			if (j < p[i].size()) swap(p[i][j], cur);
			else {
				p[i].push_back(cur);
				break;
			}
		}
	}
	while(p.back().empty()) p.pop_back();
	return p;
}

// Calculates the Robinson-Schensted-Knuth mapping in O(n sqrt(n) log(n)) time
pair<Tableau, Tableau> fastRSKMapping(vector<int> pi) {
	int n = pi.size(), k = 1;
	while(k*k < n) ++k;

	vector<int> inv_pi(n);
	for (int i = 0; i < n; ++i) inv_pi[pi[i]] = i;

	Tableau p = partialRSK(pi, k), q = partialRSK(inv_pi, k);
	reverse(pi.begin(), pi.end()); reverse(inv_pi.begin(), inv_pi.end());
	Tableau p_cols = partialRSK(pi, k), q_cols = partialRSK(inv_pi, k);

	p.resize(p_cols[0].size()); q.resize(q_cols[0].size());
	for (int i = k; i < p.size(); ++i) {
		for (int j = 0; j < p_cols.size() && p_cols[j].size() > i; ++j) {
			p[i].emplace_back(p_cols[j][i]);
			q[i].emplace_back(q_cols[j][i]);
		}
	}
	return {p, q};
}

#include <bits/stdc++.h>
using namespace std;
using Tableau = vector<vector<int>>;
using ll = long long;
const ll MOD = (ll)1e9 + 7;
ll modPow(ll a, ll b) { return (b == 0) ? 1 : ((b & 1) ? a * modPow(a, b ^ 1) % MOD : modPow(a*a % MOD, b >> 1)); }

// Calculates the number of standard Young tableau of shape lambda in O(n) using the hook-length formula
ll hookLength(int n, const vector<int>& lambda) {
	vector<int> lambda_t(lambda[0]);
	for (int i = 0, j = lambda[0]; i < n; ++i) {
		while(j >= 0 && (i + 1 == n || lambda[i+1] <= j)) lambda_t[j] = i + 1;
	}

	ll num = 1, div = 1;
	for (int i = 0; i < n; ++i) num = num * n % MOD;
	for (int i = 0; i < lambda.size(); ++i) {
		for (int j = 0; j < lambda[i]; ++j) {
			ll mult = 1 + (lambda[i] - j - 1) + (lambda_t[j] - i - 1);
			div = div * mult % MOD;
		}
	}
	return num * modPow(div, MOD - 2) % MOD;
}

Denote by $$$\lambda^-$$$ the set of partitions $$$\mu$$$ of $$$n-1$$$ such that $$$Y(\mu) \subseteq Y(\lambda)$$$. These partitions have the same Young diagram as $$$\lambda$$$ with one corner square removed.

Note that $$$c(\lambda)$$$ satisfies the recurrence $$$c(\lambda) = \sum_{\mu \in \lambda^-} c(\mu)$$$ as the maximum value has to occur in some corner cell.

We define a random walk on the Young diagram $$$Y(\lambda)$$$ of the partition. The random walk starts at a uniformly random square of $$$Y(\lambda)$$$, and in every step moves from $$$(i, j)$$$ to a uniformly random square in $$$H_\lambda(i, j) \setminus \{(i, j)\}$$$ where $$$H_\lambda(i, j)$$$ is the set of squares in $$$Y(\lambda)$$$ directly below or to the right of $$$(i, j)$$$. Thus the probability of ending up in square $$$(i', j') \in H_\lambda(i, j) \setminus \{(i, j)\}$$$ is $$$\frac{1}{h_\lambda(i, j) - 1}$$$. When $$$h_\lambda(i, j) = 1$$$, the random walk stops.

After the random walk stops, remove the square it ended up at from $$$Y$$$. We'll show that the probability that the resulting partition corresponding to the diagram after removing the square is $$$\mu \in \lambda^-$$$ with probability $$$t(\mu) / t(\lambda)$$$. Since the process always results in some $$$\mu \in \lambda^-$$$, the probabilities must sum to one, so we must then have $$$t(\lambda) = \sum_{\mu \in \lambda^-} t(\mu)$$$. Since $$$t(\mu) = c(\mu) = 1$$$ for the only partition $$$\mu$$$ of $$$1$$$, the result follows by induction.

First, note that if $$$Y(\lambda) = Y(\mu) + (i, j)$$$, we have \begin{equation} \frac{t(\mu)}{t(\lambda)} = \frac{1}{n} \prod_{k < i} \frac{h_\lambda(k, j)}{h_\lambda(k, j) — 1} \prod_{k < j} \frac{h_\lambda(i, k)}{h_\lambda(i, k) — 1} = \frac{1}{n} \prod_{k < i} \left(1 + \frac{1}{h_\lambda(k, j) — 1}\right) \prod_{k < j} \left(1 + \frac{1}{h_\lambda(i, k) — 1}\right) \end{equation}

For a random walk $$$((a_1, b_1), \dots, (a_k, b_k))$$$, let $$$A = \{a_1, \dots, a_k\}$$$ be the set of visited $$$y$$$-coordinates, and $$$B = \{b_1, \dots, b_k\}$$$ be the set of visited $$$x$$$-coordinates. Since the square the walk ends up in is $$$(a_k, b_k) = (\max A, \max B)$$$, the pair $$$(A, B)$$$ uniquely determines where the cell ends up in (though it does not uniquely determine the walk).

By $$$p(A, B | a, b)$$$, denote the probability that a random walk starting at $$$(a, b)$$$ has the given sets $$$A$$$ and $$$B$$$. If $$$\alpha = \max A$$$ and $$$\beta = \max B$$$, we claim that \begin{equation} p(A, B\ | a, b) = \prod_{i \in A \setminus \{\alpha\}} \frac{1}{h_\lambda(i, \beta) — 1} \prod_{j \in B \setminus \{\beta\}} \frac{1}{h_\lambda(\alpha, j) — 1} \end{equation} this is a easy proof by induction on $$$|A| + |B|$$$: if $$$A = \{a_1, \dots, a_k\}$$$ and $$$b = \{b_1, \dots, b_k\}$$$, $$$a_1 = a, b_1 = b$$$ and $$$a_k = \alpha, b_k = \beta$$$, \begin{equation} p(A, B\ | a, b) = \frac{1}{h_\lambda(a, b) — 1} \left(p(A \setminus \{a\}, B | a_2, b) + p(A, B \setminus \{b\} | a, b_2)\right) \end{equation} which by induction equals \begin{equation} \frac{1}{h_\lambda(a, b) — 1} \left((h_\lambda(a, \beta) — 1) + (h_\lambda(\alpha, b) — 1)\right)\prod_{i \in A \setminus \{\alpha\}} \frac{1}{h_\lambda(i, \beta) — 1} \prod_{j \in B \setminus \{\beta\}} \frac{1}{h_\lambda(\alpha, j) — 1} \end{equation} We have $$$h_\lambda(a, b) - 1 = (h_\lambda(a, \beta) - 1) + (h_\lambda(\alpha, b) - 1)$$$, thus the claimed equation for $$$p(A, B\ | a, b)$$$ holds.

Now, we are done, as \begin{equation} \sum_{a \leq \alpha, b \leq \beta} \frac{1}{n} \sum_{\stackrel{A \subseteq [n]}{\stackrel{min A = a}{max A = \alpha}}} \sum_{\stackrel{B \subseteq [n]}{\stackrel{min B = b}{max B = \beta}}} p(A, B | a, b) = \frac{1}{n} \prod_{k < \alpha} \left(1 + \frac{1}{h_\lambda(k, \beta) — 1}\right) \prod_{k < \beta} \left(1 + \frac{1}{h_\lambda(\alpha, k) — 1}\right) \end{equation} you can see this by looping $$$y$$$ from $$$\alpha$$$ to $$$1$$$ and deciding if $$$y$$$ should be in $$$A$$$, then doing the same for $$$x, \beta$$$ and $$$B$$$. The $$$\frac{1}{n}$$$ factor comes from the random selection of the starting square.