第一眼想到了minmax容斥可还行。。。
注意到 max 套个 min 很不好，我们考虑把 max 容斥掉，考虑：
\[\max(S)=\sum_{T\subseteq S} (-1)^{|T|+1} \min(T)
\]
注意到一个集合最小值的期望只和它的大小有关，对于一个大小为 \(k\) 的集合，最小值期望显然是 \({1\over x^k}\sum_{i=1}^x i^k\)，我们只需要对于每个长度，统计出并出它的方案的系数和。考虑 dp。
首先，如果一个区间包含令一个区间，那它显然没用，我们直接删去它。这样我们所有区间都没有包含关系。然后按右端点排序，设 \(f_{i,j,0/1}\) 表示考虑前 \(i\) 个区间，强制选第 \(i\) 个，区间并长度为 \(j\)，选了偶数/奇数个区间的方案数，那么枚举上一个选的区间是什么，\(O(q^3)\) 的背包是显然的。
考虑优化，注意到本质不同的转移只有两种，设当前区间为 \([l,r]\)，上一个区间为 \([l',r']\)，如果 \(r'<l\)，显然会从 \(k-(r-l+1)\) 转移过来，这是个定值，由于 \(r\) 递增，满足条件的上一个区间是一段前缀，可以记录前缀和转移。如果 \(r'>l\)，则区间有交。由于区间互不包含，新增的段一定是最右边的一段，会从 \(k-r+r'\) 转移过来。我们记 \(h_{i,j,k-r}=f_{i,j,k}\)，显然会从 \(h\) 的 \(k-r\) 转移过来，这也可以做前缀和转移。
#include <bits/stdc++.h>
using namespace std;
const int N = 2005, B = 2100, mod = 666623333;
inline int read() {
register int s = 0, f = 1; register char ch = getchar();
while (!isdigit(ch)) f = (ch == '-' ? -1 : 1), ch = getchar();
while (isdigit(ch)) s = (s * 10) + (ch & 15), ch = getchar();
return s * f;
}
inline int power(int a, int b) {
int t = 1, y = a, k = b;
while (k) {
if (k & 1) t = 1ll * t * y % mod;
y = 1ll * y * y % mod; k >>= 1;
} return t;
}
struct node {
int l, r;
} p[N];
inline bool cmp(node a, node b) {
return a.r < b.r;
}
bool vis[N];
int f[2005][2][2005], g[2005][2][2005], n, m, q, pw[N], t[N], res[N];
unordered_map<int, int> h[2005][2];
int main() {
n = read(); m = read(); q = read();
for (int i = 1; i <= q; ++i) {
p[i].l = read(); p[i].r = read();
} int top = 0;
for (int i = 1; i <= q; ++i) {
for (int j = 1; j <= q; ++j) {
if (i == j) continue;
if (p[i].l <= p[j].l && p[j].r <= p[i].r) {
if (vis[j] && p[i].l == p[j].l && p[i].r == p[j].r) continue;
vis[i] = 1; break;
}
}
}
for (int i = 1; i <= q; ++i)
if (!vis[i]) p[++top] = p[i];
int cnt = 0;
sort(p + 1, p + top + 1, cmp);
f[0][0][0] = 1; g[0][0][0] = 1; h[0][0][0] = 1;
for (int i = 1; i <= top; ++i) {
for (int j = 0; j <= 1; ++j) {
int L = 0;
for (int k = 1; k < i; ++k)
if (p[k].r < p[i].l) L = k;
else break;
memcpy(g[i][j], g[i - 1][j], sizeof g[i][j]);
h[i][j] = h[i - 1][j];
for (int k = p[i].r - p[i].l + 1; k <= p[i].r; ++k) {
f[i][j][k] += g[L][j ^ 1][k - (p[i].r - p[i].l + 1)];
f[i][j][k] += h[i - 1][j ^ 1][k - p[i].r] - h[L][j ^ 1][k - p[i].r];
if (f[i][j][k] >= mod) f[i][j][k] -= mod;
if (f[i][j][k] < 0) f[i][j][k] += mod;
g[i][j][k] += f[i][j][k];
if (g[i][j][k] >= mod) g[i][j][k] -= mod;
h[i][j][k - p[i].r] += f[i][j][k];
if (h[i][j][k - p[i].r] >= mod) h[i][j][k - p[i].r] -= mod;
if (j) res[k] += f[i][j][k];
else res[k] -= f[i][j][k];
if (res[k] >= mod) res[k] -= mod;
if (res[k] < 0) res[k] += mod;
}
}
} int ans = 0;
for (int i = 1; i <= m; ++i) pw[i] = 1;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= m; ++j) {
pw[j] = 1ll * j * pw[j] % mod;
t[i] += pw[j]; if (t[i] >= mod) t[i] -= mod;
} t[i] = 1ll * t[i] * power(power(m, mod - 2), i) % mod;
ans += 1ll * t[i] * res[i] % mod;
if (ans >= mod) ans -= mod;
} printf("%d\n", ans); return 0;
}