@description@
输入 n(n ≤ 22) 个点，m(m ≤ 8000) 个边。每个边连接着点 (si, ei)，有两个长度 fi, ri。
问对于每个点 k，有多少条路径（不一定是简单路径）由 t (t ≤ 10^9) 条边组成，从 k 开始，并且以 k 结束；并且路径上所有边 f 的和 mod n 为 x；并且路径上所有边 r 的和 mod (n − 1) 为 y。
对于每一个 (x, y) 都要计算。
方案数 mod 1163962801 输出。
input
输入的第一行包含三个整数 n、m 和 t。n, m, t 的含义如上，注意 n 是点数，也是模数。如果重复经过同一条边，需要计算两次。
之后的 m 行里，每一行表示一条边。其中第 i 行包含四个整数，分别代表 si, ei, fi, ri。含义如上。
可能有重边自环。
output
输出被分为 n 个部分。其中第 k 个部分 (1 ≤ k ≤ n) 包含 n 行，而每一行包含 n − 1 个整数。
第 k 个输出部分的第 i 行的第 j 个数应当表示从点 k 出发，走过 t 条边之后回到点 k，路径上所有边 f 的和 mod n 为 i；并且路径上所有边 r 的和 mod (n − 1) 为 j。
输出的所有数字都应当对 1163962801 取模。
sample input
3 6 6
1 2 0 0
2 1 1 1
1 3 1 2
3 1 2 1
2 3 5 5
3 2 10 10
sample output
1 5
0 10
1 5
1 5
0 10
1 5
1 8
0 10
1 2
@solution@
首先观察到 n 的范围很小，t 的范围很大，并且可以经过重复的边。不难想到使用矩阵加速。
一个很 naive 的想法是，我们对于矩阵的每一个位置 (i, j) 记录 n*(n-1) 个二元组 (x, y)，表示从 i 到 j 路径上所有边 f 的和 mod n 为 x，r 的和 mod (n − 1) 为 y 的方案数。
这样做矩阵乘法 (C = A*B) 的时候，(C(i, j) = sum A(i, k)*B(k, j)) 实际上就是 (C(i, j, (x+p) mod n, (y+q) mod (n-1)) = sum(A(i, k, x, y) * B(k, j, p, q)))。
要做 (O(n^3)) 次乘法，每次乘法 (O(n^4))。时间复杂度 (O(n^7*log k))。虽然 n 很小但还是太慢了。
观察我们乘法的转移形式，不难发现它是一个循环卷积的形式，可以使用 fft 优化。
优化过后每次乘法只需要 (O(n^2*log n)) 的时间，时间复杂度降为 (O(n^5*log n*log k))。
但是注意到模数不是 998244353 而是 1163962801，所以这里的 fft 必须要使用一些常数较大的技巧来避免溢出。
所以实际运行时效果远不如 (O(n^5*log n*log k)) 那么优秀。
考虑我们平常使用多项式作快速幂，是将其转为点值形式，使用点值做快速幂。这个地方是否也可以运用这一方法呢。
对于长度为 n 的循环卷积，无法将其长度拓展为 2 的幂，所以无法使用一般的 fft。
但是，代入 n 次单位根很多性质依然可以满足。或者说，除了不能使用分治进行 fft，其他性质都可以满足。
注意到这个地方的 n 很小，我们可以直接暴力 (O(n^2)) 求值，暴力 (O(n^2)) 还原。
观察模数 mod = 1163962801。想一想为什么可以使用 p = 998244353 进行数论变换，因为 p-1 恰好是 2^23 的倍数，所以它有 2^x(x ≤ 23) 次单位根。
对 mod - 1 进行质因数分解，它等于 (2^4*3^2*5^2*7*11*13*17*19)。我们发现它是 1~22 所有数字的公倍数，因此它有 1~22 次单位根。
最后一个问题，这个地方是一个二维的卷积形式。
因此，我们求值时是代入一个二元组 ((w_n^i, w_{n-1}^j))，一共有 n*(n-1) 对二元组要求值。
最后逆变换还原时，插值 ((w_n^{-i}, w_{n-1}^{-j})) 即可。直观理解起来不难。
时间复杂度为 (O(n^5*log k))。
@accepted code@
#include<cstdio>
const int G = 46;
const int MAXN = 22 + 5;
const int MAXM = 10000 + 5;
const int MOD = 1163962801;
inline int add(int x, int y) {return (1LL*x + 1LL*y)%MOD;}
inline int mul(int x, int y) {return 1LL*x*y%MOD;}
int pow_mod(int b, int p) {
int ret = 1;
while( p ) {
if( p & 1 ) ret = mul(ret, b);
b = mul(b, b);
p >>= 1;
}
return ret;
}
struct matrix{
int m[MAXN][MAXN], r, c;
friend matrix operator * (const matrix &A, const matrix &B) {
matrix C; C.r = A.r, C.c = B.c;
for(int i=0;i<C.r;i++)
for(int j=0;j<C.c;j++) {
C.m[i][j] = 0;
for(int k=0;k<A.c;k++) {
C.m[i][j] = add(C.m[i][j], mul(A.m[i][k], B.m[k][j]));
}
}
return C;
}
};
matrix mpow(matrix A, int p) {
matrix ret; ret.r = ret.c = A.r;
for(int i=0;i<A.r;i++)
for(int j=0;j<A.c;j++)
ret.m[i][j] = (i == j);
while( p ) {
if( p & 1 ) ret = ret * A;
A = A * A;
p >>= 1;
}
return ret;
}
void func(int a[], int n) {
a[0] = 1, a[1] = pow_mod(G, (MOD - 1)/n);
for(int i=2;i<n;i++) a[i] = mul(a[i - 1], a[1]);
}
int wx[MAXN], wy[MAXN];
int u[MAXM], v[MAXM], x[MAXM], y[MAXM];
int t, n, m;
matrix get_matrix(int i, int j) {
matrix A; A.r = A.c = n;
for(int k=0;k<A.r;k++)
for(int l=0;l<A.c;l++)
A.m[k][l] = 0;
for(int k=1;k<=m;k++)
A.m[u[k]][v[k]] = add(A.m[u[k]][v[k]], mul(pow_mod(wx[i], x[k]), pow_mod(wy[j], y[k])));
return A;
}
int f[MAXN][MAXN][MAXN];
int ans[MAXN][MAXN];
void func2(int i) {
for(int j=0;j<n;j++)
for(int k=0;k<n-1;k++)
ans[j][k] = 0;
for(int j=0;j<n;j++)
for(int k=0;k<n-1;k++)
for(int p=0;p<n;p++)
for(int q=0;q<n-1;q++) {
ans[j][k] = add(ans[j][k], mul(f[i][p][q], mul(pow_mod(wx[j], n-p), pow_mod(wy[k], (n-1)-q))));
}
int inv = pow_mod(n*(n - 1), MOD - 2);
for(int j=0;j<n;j++)
for(int k=0;k<n-1;k++)
printf("%d%c", mul(ans[j][k], inv), (k + 1 == n - 1) ? '
' : ' ');
}
int main() {
scanf("%d%d%d", &n, &m, &t);
for(int i=1;i<=m;i++) {
scanf("%d%d%d%d", &u[i], &v[i], &x[i], &y[i]);
u[i]--, v[i]--, x[i] %= n, y[i] %= (n-1);
}
func(wx, n), func(wy, n - 1);
for(int i=0;i<n;i++)
for(int j=0;j<n-1;j++) {
matrix M = mpow(get_matrix(i, j), t);
for(int k=0;k<n;k++)
f[k][i][j] = M.m[k][k];
}
for(int i=0;i<n;i++)
func2(i);
}
@details@
这道题的模数，加法刚好溢出 int。
我就说为什么它会出现负数……
2017 冬令营的老题了……