@codechef - [emailprotected] Chef and Bike

@[email?protected]

输入 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

@[email?protected]

首先观察到 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,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 [email?protected]

#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) ? ‘n‘ : ‘ ‘);
}
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);
}

@[email?protected]

这道题的模数，加法刚好溢出 int。
我就说为什么它会出现负数……

2017 冬令营的老题了……