[From this](https://www.luogu.com.cn/blog/_post/254278)
考虑经典容斥，枚举哪些人的路径相交了。
相交的路径可以转化为，在相交点将两个人之后的路径取反，最终是第一个人走到第二个人的终点去了，第二个人走到第一个人的终点去了。
直接这样容斥复杂度是 $O(n!)$ 的。
进一步观察可以发现，容斥的系数就是逆序对的奇偶性，这和行列式的定义完全一样。
我们可以求出每个起点到每个终点的方案数，组成一个矩阵，求出来行列式就是答案。
妙啊！
顺便打磨了一下之前的模板。
尝试修改码风为大括号换行...
```cpp
namespace FAC
{
const int MAXN=2000005;
struct _Factor {
int fac[MAXN], ifac[MAXN];
_Factor() {
fac[0]=1;
for(int i=1; i<MAXN; ++i) fac[i]=mul(fac[i-1], i);
ifac[MAXN-1]=mdinv(fac[MAXN-1]);
for(int i=MAXN-1; i>=1; --i) ifac[i-1]=mul(ifac[i], i);
}
il int C(const int n, const int m) {return (n<m||m<0)?0:mul(fac[n], ifac[m], ifac[n-m]);}
il int inv(const int n) {return mul(ifac[n], fac[n-1]);}
} Factor;
il int C(const int n, const int m) {return Factor.C(n, m);}
il int inv(const int n) {return Factor.inv(n);}
} using namespace FAC;
namespace MAT {
const int SZ=101;
struct Mat {
int sz;
int m[SZ][SZ];
il void clear() {mem(m, 0);}
il Mat () {sz = 0; clear();}
il int *operator[](const int x) {return this->m[x];}
il Mat operator*(const Mat &x)const {
Mat res; res.sz = sz;
for (int i = 1; i <= sz; i++)
for (int k = 1; k <= sz; k++)
for (int j = 1; j <= sz; j++)
inc(res.m[i][j], mul(m[i][k], x.m[k][j]));
return res;
}
il void operator*=(const Mat &x) {*this = (*this) * x;}
il void operator+=(const Mat &x) {
for (int i = 1; i <= sz; ++i)
for (int j = 1; j <= sz; ++j)
inc(m[i][j], x.m[i][j]);
}
il Mat operator+(const Mat &x)const {Mat res = *this; res += x; return res;}
il void operator-=(const Mat &x) {
for (int i = 1; i <= sz; ++i)
for (int j = 1; j <= sz; ++j)
dec(m[i][j], x.m[i][j]);
}
il Mat operator-(const Mat &x)const {Mat res = *this; res -= x; return res;}
il void print() {for(int i=1; i<=sz; ++i) prt(m[i]+1, sz);}
il Mat e() {clear(); for (int i = 1; i <= sz; i++) m[i][i] = 1; return *this;}
il Mat qpow(int x) {
if (x == 0) return Mat().e();
Mat res, mul = *this;
res = *this; --x;
for (; x; x >>= 1, mul *= mul) if (x & 1) res *= mul;
return res;
}
il int det() {
int ans = sz == 1 ? m[1][1] : 1;
for (int i = 1; i <= sz; ++i) {
for (int j = i + 1, t; j <= sz; ++j) {
while (m[j][i]) {
t = m[i][i] / m[j][i];
for (int k = i; k <= sz; ++k) {
dec(m[i][k], mul(m[j][k], t));
swap(m[i][k], m[j][k]);
}
ans = md - ans;
}
}
if (m[i][i] == 0) return 0;
else ans=mul(ans, m[i][i]);
}
return ans;
}
};
} using namespace MAT;
const int N=105;
Mat mat;
int n, m;
int calc(int a, int b) {
int x=b-a, y=n-1;
return C(x+y, x);
}
void clear() {
mat.clear();
}
int a[N], b[N];
il void solve() {
mat.sz=m;
for(int i=1; i<=m; ++i)
for(int j=1; j<=m; ++j)
mat[i][j]=calc(a[i], b[j]);
out(mat.det());
clear();
}
```