#include <cstdio>
#include <algorithm>
#include <cstring>
typedef long long ll;
using namespace std;
const int MAXN = 1005, mod = 1004535809, G = 3, invG = (mod + 1) / 3;
int n, m, r, k, q;
ll fac[MAXN], inv[MAXN];
ll Fstpw(ll a, int b) {
    ll res = 1;

    while (b) {
        if (b & 1)
            res = res * a % mod;

        b >>= 1;
        a = a * a % mod;
    }

    return res;
}
void ntt(int *a, int n, int tp) {
    int bit = 0;

    while (1 << bit < n)
        bit++;

    static int rev[MAXN << 2];

    for (int i = 1; i < n; i++) {
        rev[i] = rev[i >> 1] >> 1 | ((i & 1) << bit - 1);

        if (i < rev[i])
            swap(a[i], a[rev[i]]);
    }

    for (int mid = 1; mid < n; mid <<= 1) {
        ll w = 1, w1 = Fstpw(tp == 1 ? G : invG, (mod - 1) / mid / 2);

        for (int j = 0; j < mid; j++, w = w * w1 % mod)
            for (int i = 0; i < n; i += mid * 2) {
                int x = a[i + j], y = w * a[i + j + mid] % mod;
                a[i + j] = (x + y) % mod;
                a[i + j + mid] = (x - y + mod) % mod;
            }
    }

    if (tp == -1) {
        ll t = Fstpw(n, mod - 2);

        for (int i = 0; i < n; i++)
            a[i] = a[i] * t % mod;
    }

    return ;
}
void NTT(int *a, int n, int *b, int m, int *res, int k) {
    int siz = 1;

    while (siz <= n + m)
        siz <<= 1;

    static int f[MAXN << 2], g[MAXN << 2];

    for (int i = 0; i < siz; i++) {
        f[i] = i <= n ? a[i] : 0;
        g[i] = i <= m ? b[i] : 0;
    }

    ntt(f, siz, 1);
    ntt(g, siz, 1);

    for (int i = 0; i < siz; i++)
        f[i] = 1ll * f[i] * g[i] % mod;

    ntt(f, siz, -1);

    for (int i = 0; i <= k; i++)
        res[i] = f[i];

    return ;
}
struct Mat1 {
    int a[5][5];
    int *operator [](int x) {
        return a[x];
    }
    void print() {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++)
                printf("%d ", a[i][j]);

            puts("");
        }

        puts("");
    }
} a;
struct Mat2 {
    int a[5][5][MAXN];
    auto operator [](int x) {
        return a[x];
    }
} b, c;
Mat1 operator *(Mat1 a, Mat1 b) {
    static Mat1 c;

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++) {
            c[i][j] = 0;

            for (int k = 0; k < n; k++)
                c[i][j] = (c[i][j] + 1ll * a[i][k] * b[k][j]) % mod;
        }

    return c;
}
Mat2 operator +(Mat2 a, Mat2 b) {
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k <= r; k++)
                a[i][j][k] = (a[i][j][k] + b[i][j][k]) % mod;

    return a;
}
Mat2 operator *(Mat2 a, Mat1 b) {
    static Mat2 c;

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int l = 0; l <= r; l++) {
                c[i][j][l] = 0;

                for (int k = 0; k < n; k++)
                    c[i][j][l] = (c[i][j][l] + 1ll * a[i][k][l] * b[k][j]) % mod;
            }

    static int f[MAXN];

    for (int i = 0, t = 1; i <= r; i++, t = 1ll * t * m % mod)
        f[i] = t * inv[i] % mod;

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            NTT(c[i][j], r, f, r, c[i][j], r);

    return c;
}
int main() {
    //freopen("in","r",stdin);
    scanf("%d%d%d%d%d", &n, &m, &r, &k, &q);
    fac[0] = 1;

    for (int i = 1; i <= r; i++)
        fac[i] = fac[i - 1] * i % mod;

    inv[r] = Fstpw(fac[r], mod - 2);

    for (int i = r; i; i--)
        inv[i - 1] = inv[i] * i % mod;

    while (m--) {
        int u, v;
        scanf("%d%d", &u, &v);
        u--, v--;
        a[u][v]++;

        if (u != v)
            a[v][u]++;
    }

    for (int i = 0; i < n; i++)
        b[i][i][0] = c[i][i][0] = 1;

    m = 1;

    while (k) {
        if (k & 1)
            c = c * a + b;

        k >>= 1;
        b = b + b * a;
        a = a * a;
        m <<= 1;
        //a.print();
    }

    while (q--) {
        int s, t;
        scanf("%d%d%d", &s, &t, &r);
        int ans = fac[r] * c[s - 1][t - 1][r] % mod;

        if (s == t && r == 0)
            ans--;

        printf("%d\n", ans);
    }

    return 0;
}

地牢路径之谜：路径长度幂次和求解

题目分析

本题要求计算从起点$S_i$到终点$T_i$，所有长度$\leq k$的不同路径的路径长度的$r_i$次方和。路径由长度为1的边构成，图中可能存在重边和自环（但不存在长度为0的自环路径）。由于$k$最大可达$10^9$，直接枚举路径长度计算是不可能的，需要结合矩阵快速幂、生成函数和数论变换等高级算法。

核心思路

1. 路径计数与邻接矩阵

• 图中长度为$d$的路径数可通过邻接矩阵的幂表示：设邻接矩阵为$A$（$A[u][v]$表示$u$到$v$的边数），则$A^d[u][v]$表示从$u$到$v$长度为$d$的路径数。
• 目标是计算$\sum_{d=1}^{\min(k, d_{\text{max}})} A^d[S][T] \cdot d^r$，即路径数乘以路径长度的$r$次方之和。

2. $d^r$的生成函数表示

由于$r$最大为1000，直接维护$d^r$的累加困难，需利用生成函数将$d^r$转化为可高效计算的形式：

• 利用组合数学性质：$d^r = \sum_{m=0}^r S(r, m) \cdot m! \cdot \binom{d}{m}$，其中$S(r, m)$为斯特林数。
• 通过生成函数系数维护$\binom{d}{m}$的累加，再通过卷积（NTT）转换为$d^r$的系数。

3. 矩阵快速幂求和

由于$k$极大，采用类似等比数列求和的快速幂思想，维护一个包含生成函数系数的矩阵结构，高效计算$\sum_{d=0}^k A^d$相关的系数：

• 定义特殊矩阵结构Mat2（$5 \times 5 \times (r+1)$），其中Mat2[i][j][l]表示从$i$到$j$的路径中，与$\binom{d}{l}$相关的累加系数。
• 通过矩阵乘法和加法，结合快速幂迭代，将求和范围从$k$分解为二进制幂次，降低时间复杂度。

算法步骤

1. 预处理

• 计算阶乘fac和逆阶乘inv，用于生成函数系数转换。
• 初始化邻接矩阵A，记录图中边的数量。

2. 矩阵定义与运算

• Mat1：普通$5 \times 5$矩阵，用于表示邻接矩阵及其幂次（路径数）。
• Mat2：扩展矩阵，每个元素是长度为$r+1$的数组，存储生成函数系数。
• 矩阵乘法：Mat2 * Mat1通过卷积（NTT）维护生成函数系数，实现路径数与系数的结合。
• 矩阵加法：Mat2 + Mat2实现系数的累加，对应路径长度范围的合并。

3. 快速幂迭代

• 模仿二进制快速幂，将$k$分解为$2^0, 2^1, ..., 2^t$的组合。
• 迭代维护b（当前幂次的系数矩阵）和c（累加的系数矩阵），通过c = c * a + b和b = b + b * a更新，最终c包含$\sum_{d=0}^k A^d$的生成函数系数。

4. 询问处理

• 对于每个询问$(S, T, r_i)$，从c中提取对应系数，乘以$r_i!$（还原$d^r$的系数），得到结果。
• 特殊处理：当$S=T$且$r_i=0$时，需减去长度为0的路径（题目规定不存在）。

代码解析

关键函数与结构

• ntt/NTT：快速数论变换，用于高效计算生成函数的卷积。
• Mat1乘法：普通矩阵乘法，计算邻接矩阵的幂次（路径数）。
• Mat2运算：结合卷积的扩展矩阵运算，维护生成函数系数。
• 快速幂循环：分解$k$并更新矩阵，实现大范围求和。

复杂度分析

• 预处理：$O(r)$（阶乘、逆元）。
• 矩阵运算：由于$n=5$，单次Mat1乘法为$O(n^3)$，Mat2运算结合NTT为$O(n^3 \cdot r \log r)$。
• 快速幂迭代：$O(\log k)$次矩阵运算，总复杂度为$O(\log k \cdot n^3 \cdot r \log r)$，可应对$k=10^9$和$q=25000$的规模。

算法标签

• 矩阵快速幂
• 生成函数
• 快速数论变换（NTT）
• 组合数学
• 动态规划

通过上述方法，成功将大$k$的路径求和问题转化为可高效计算的矩阵运算，结合生成函数和数论变换处理幂次项，最终在规定时间内求解所有询问。