同学打月赛后告诉我有道很离谱的题，遂做，确实挺离谱的，这里讲一下自己做的思路。
题意
- 给定 $n$，平面直角坐标系上有点 $(0,0)\sim(0,n-1),(1,0)\sim(1,n-1)$。对于一个 $n$ 排列 $p$，对每个 $i$ 连线段 $(0,i)-(1,p_i)$。
- 考虑从每个横坐标为 $0$ 的点出发，依附于这些线段各走一条折线（或线段）到达一个横坐标为 $1$ 的点，这些折线（或线段）仅能在原先 $n$ 条线段的交点上重合。
- 求所有 $p$ 的方案数之积${}\bmod{335544323}$。
- $n\le 10^6$。
Step 1
方便起见称原题中完美匹配（简称匹配）的连线为 “线段”，染色的折线（或线段）为 “路径”。
模拟一些例子，发现三点：
- 路径不会 “往回折” 且线段没有任何部分不被经过。
- 选择一条路径在（最后）一个交点的走向，剩余的东西并没有本质不同，感性地说像是有一个乘法原理。
- 三线段交于同一点和三线段交于三个点的方案数是不同的，也就是说这个几何背景是有用的。
可以证明第一点：考虑任一横坐标 $x_0$ 满足原先 $n$ 条线段的交点均不在 $x=x_0$ 上，这时原图有 $n$ 个点在 $x=x_0$ 上，而每一条路径至少穿过一次 $x=x_0$。如果有一条路径往回折，那就穿过某个 $x=x_0$ 大于一次，这就不可能存在 $n$ 条两两交长度 $=0$ 的路径了。
进一步分析，发现一个 “$k$ 线共点” 的情况会对答案贡献 $k!$。因为这个局部相当于一个 $k$ 阶匹配以决定每条路径的走向，且每种方案显然都是合法的。因此最有可能的思路大概就是对于每个可能的交点，求出每种匹配情况下有几条线段经过它，然后阶乘乘起来贡献给答案。当然不能像这样直接算。
显然贡献值只取决于该交点最多被几条线段经过（显然这些线段的端点均不同，所以相当于是一个选择问题），故可以求出 $w_k$ 表示将所有可能的 $n^2$ 条线段画出后 $k$ 线共点的情况数，$p_k$ 表示一个至多作为 $k$ 条线段交点的点在所有匹配情况下对答案的总贡献。就有：
$$ ans=\prod_{k=2}^n p_k^{w_k} $$
Step 2
求 $w_k$。考虑刻画线段交于一点的情况。
（为了美观横过来了）由图中相似关系可见，两边选的相邻点距离要对应成比例。
第一反应是直接硬着对极大共点选择（一些线段交于同一点，且不能再多选一条经过该交点的线段。显然这些线段的端点均不同）计数，官方解法有更方便的容斥做法，这里就算给一个备选方案了。极大共点选择的条件为：
- 两边选的点纵坐标分别成等差数列；
- 不能再往左右扩展；
- 公差互质。
首先发现互质条件无法直接处理，因此莫反，枚举公差的公因数 $g$，贡献乘 $\mu(g)$ 即可。如果考虑左右部分别长为 $x$ 和 $y$ 的极长等差数列，那么画图可得它们能构成的极大共点选择包括：$2$ 个 $2$ 线共点，$2$ 个 $3$ 线共点，……，$2$ 个 $\min(x,y)-1$ 线共点，$\lvert x-y\rvert+1$ 个 $\min(x,y)$ 线共点。
那么只需对于每个 $g$，把每种为 $g$ 倍数的公差所能形成的极长等差数列长度及其数量求出来，统计并扫描处理即可。由于一个公差 $d$ 至多只能产生 $\lfloor (n-1)/d\rfloor+1$ 与 $\lfloor (n-1)/d\rfloor$ 两种长度，而扫描可以优化到线性（枚举长较小者，部分细节略去），故这部分可以做到 $\mathrm{O}(n\log n)$。
Step 3
直接求 $p_k$ 的思路为：对于每个 $i=1,\cdots,k$，求出有多少种匹配恰好出现 $i$ 条线段经过该点，记为 $f_i$，则：
$$ p_k=\prod_{i=1}^k (i!)^{f_i} $$
$f_i$ 无法直接计算，记 $g_i$ 表示钦定 $i$ 条线段经过该点，其他随便匹配的方案数（可能会算重）。则：
$$ \binom ki(n-i)!=g_i=\sum_{j=i}^k\binom jif_i $$
二项式反演得：
$$ f_i=\sum_{j=i}^k(-1)^{j-i}\binom ji\binom kj(n-j)! $$
然后就会发现没法快速求出 $p_i$（如果有大佬会的能教一下吗/kel
Step 4
既然没法快速求 $p_i$，那么考虑利用最终只需求一个积这个性质，把所有 $k$ 放在一起求。记：
$$ c_i=\sum_{k=i}^nw_k\sum_{j=i}^k(-1)^{j-i}\binom ji\binom kj(n-j)! $$
则：
$$ ans=\prod_{i=2}^n(i!)^{c_i} $$
而 $c_i$ 是可以快速求的：
$$ c_i=\sum_{j=i}^n(-1)^{j-i}\binom ji(n-j)!\left(\sum_{k=j}^n w_k\binom kj\right) $$
两次差卷积即可。细节略去。
Step 5
最后一个问题是 $c_i$ 是在指数上的。注意到原模数为质数 $335544323=2^{26}\times 5+3$，因此如果能求出 $c_i$ 模 NTT 质数 $2^{25}\times 5+1$ 与 $c_i$ 的奇偶性即可用 CRT 的思路还原 $c_i\bmod (2^{26}\times 5+2)$。
一个性质是 $w_i\equiv[i=n]\pmod 2$（显然 $w_n=1$）。这是因为除了两部公差均为 $1$ 的情况，其余所有极大共点选择可以对应其两部交换后的另一个极大共点选择。那么：
$$ \begin{aligned}c_i&\equiv\sum_{j=i}^n\binom ji\binom nj(n-j)!\\ &\equiv n\binom{n-1}i+\binom ni\\ &\equiv(n-i+1)[i\subseteq n]\pmod 2\end{aligned} $$
最终时间复杂度为 $\mathrm{O}(n\log n)$。
代码
一个优化：两次卷积分别要求类似于 $\text{DFT}(\mathrm{e}^x)$ 以及 $\text{DFT}(\mathrm{e}^{-x})$ 形式的东西，可以省掉一次。共 $5$ 次 DFT。
void NTT(int n, int *a, int o) {
for (int i=0; i<n; i++)
if (i < f[i]) swap(a[i], a[f[i]]);
for (int i=1, _; i<n; i<<=1)
for (int j=0; j<n; j+=i<<1)
for (int k=j; k<j+i; k++)
_ = 1ll * g[i][k-j] * a[k+i] % P, a[k+i] = R(a[k], P - _), a[k] = R(a[k], _);
if (o) {
reverse (a+1, a+n);
for (int i=0; i<n; i++)
a[i] = 1ll * a[i] * I % P;
}
}
int main() {
cin >> n;
for (int i=mu[1]=1; i<=n; i++) if (mu[i])
for (int j=i<<1; j<=n; j+=i)
mu[j] -= mu[i];
for (int i=1; i<n; i++)
L[i] = (n - 1) / i + 1, C[i] = (n - 1) % i + 1;
for (int g=1; g<n; g++) {
if (! mu[g]) continue;
int m = 0, sl = 0, sc = 0;
for (int i=g; i<n; i+=g) {
if (m && l[m-1] == L[i]) c[m-1] = R(c[m-1], C[i]);
else if (l[m] == L[i]) c[m] = R(c[m], C[i]);
else l[++m] = L[i], c[m] = C[i];
if (C[i] < i && L[i] > 2)
if (l[m] == L[i] - 1) c[m] = R(c[m], i - C[i]);
else l[++m] = L[i] - 1, c[m] = i - C[i];
}
for (int i=1; i<=m; i++) {
w[l[i]] = (w[l[i]] + ((sl + 1ll * (P - l[i] + 1) * sc << 1) + c[i]) % P * (~ mu[g] ? c[i] : P - c[i])) % P;
d[l[i]-1] = (d[l[i]-1] + (2ll * sc + c[i]) * (~ mu[g] ? c[i] : P - c[i])) % P;
sl = (sl + 1ll * l[i] * c[i]) % P, sc = R(sc, c[i]);
}
}
for (int i=n; i>1; i--)
d[i] = R(d[i], d[i+1]), w[i] = R(w[i], R(d[i], d[i]));
for (int i=F[0]=1; i<=n; i++)
F[i] = 1ll * F[i-1] * i % P;
IF[n] = qpow(F[n], P - 2);
for (int i=n-1; ~i; i--)
IF[i] = 1ll * IF[i+1] * (i+1) % P;
while (m < (n << 1) - 3) m <<= 1;
I = qpow(m, P - 2);
for (int i=0; i<m; i++)
f[i] = f[i>>1] >> 1 | (i & 1 ? m >> 1 : 0);
for (int i=1; i<m; i<<=1) {
g[i] = new int [i] {1};
if (i > 1) {
g[i][1] = qpow(3, (P - 1) / (i << 1));
for (int j=2; j<i; j++) g[i][j] = 1ll * g[i][j-1] * g[i][1] % P;
}
}
for (int i=0; i<n-1; i++)
a[i] = 1ll * F[n-i] * w[n-i] % P, b[i] = IF[i];
NTT(m, a, 0), NTT(m, b, 0);
for (int i=0; i<m; i++)
a[i] = 1ll * a[i] * b[i] % P;
NTT(m, a, 1);
for (int i=0; i<m; i++)
a[i] = i < n - 1 ? 1ll * F[i] * a[i] % P : 0;
for (int i=0; i<m>>1; i++)
swap(b[i], b[i|m>>1]);
NTT(m, a, 0);
for (int i=0; i<m; i++)
a[i] = 1ll * a[i] * b[i] % P;
NTT(m, a, 1);
for (int i=0; i<n-1; i++)
c[n-i] = 1ll * IF[n-i] * a[i] % P;
for (int i=2; i<=n; i++)
if ((i & n) == i & n - i - 1 & 1 ^ c[i] & 1)
c[i] += P;
for (int i=n; i>1; i--) {
if ((c[i] += c[i+1]) >= P2 - 1) c[i] -= P2 - 1;
ans = 1ll * ans * qpow_(i, c[i]) % P2;
} cout << ans;
}