#include <bits/stdc++.h>

using i64 = long long;

constexpr int P = 998244353;

int power(int a, int b) {
    int res = 1;
    for (; b != 0; b /= 2, a = 1LL * a * a % P) {
        if (b % 2 == 1) {
            res = 1LL * res * a % P;
        }
    }
    return res;
}

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);
    
    int c, n, m;
    i64 k;
    std::cin >> c >> n >> m >> k;
    
    std::vector<int> fac(n + 1), invfac(n + 1), inv(n + 1);
    fac[0] = 1;
    for (int i = 1; i <= n; i++) {
        fac[i] = 1LL * fac[i - 1] * i % P;
    }
    invfac[n] = power(fac[n], P - 2);
    for (int i = n; i >= 1; i--) {
        invfac[i - 1] = 1LL * invfac[i] * i % P;
        inv[i] = 1LL * invfac[i] * fac[i - 1] % P;
    }
    
    std::vector p(m, std::vector<int>(n));
    for (int i = 0; i < n * m; i++) {
        int u, v, w;
        std::cin >> u >> v >> w;
        u--, v--, w--;
        p[w][u] = v;
    }
    
    std::vector<std::vector<std::vector<int>>> comp(1), ncomp;
    comp[0].resize(n);
    for (int i = 0; i < n; i++) {
        comp[0][i] = {i};
    }
    
    std::vector<int> vis(n), bel(n);
    for (int w = 0; w < m; w++) {
        ncomp.clear();
        for (auto &c : comp) {
            for (int i = 0; i < c.size(); i++) {
                vis[i] = 0;
                for (auto x : c[i]) {
                    bel[x] = i;
                }
            }
            std::map<std::pair<int, int>, std::vector<std::vector<int>>> mp;
            for (int i = 0; i < c.size(); i++) {
                if (vis[i]) {
                    continue;
                }
                std::vector<int> a;
                int x = c[i][0];
                while (!vis[bel[x]]) {
                    int j = bel[x];
                    std::rotate(c[j].begin(), std::find(c[j].begin(), c[j].end(), x), c[j].end());
                    a.insert(a.end(), c[j].begin(), c[j].end());
                    vis[j] = 1;
                    x = p[w][x];
                }
                int pos = std::find(a.begin(), a.end(), x) - a.begin();
                mp[{a.size(), pos}].push_back(a);
            }
            for (auto &[_, c] : mp) {
                ncomp.emplace_back(std::move(c));
            }
        }
        std::swap(comp, ncomp);
    }
    
    int ans = 1;
    for (auto &c : comp) {
        int n = c.size();
        int siz = c[0].size();
        std::vector<int> f(n + 1);
        for (int s = 1; s <= n; s++) {
            std::vector<int> v;
            int ans = 0;
            auto dfs = [&](auto &&self, int x) -> void {
                v.push_back(x);
                if (x == s * siz) {
                    int m = v.size();
                    int res = 1;
                    for (int i = 1; i < m; i++) {
                        int x = v[i - 1], y = v[i];
                        res = 1LL * res * power(x, s * siz / x - s * siz / y + 1) % P * power(inv[y / x], s * siz / y) % P;
                    }
                    std::vector<int> p(m), q(m + 1);
                    p[m - 1] = 1;
                    q[0] = 1;
                    for (int i = 0; i < m; i++) {
                        for (int j = i + 1; j >= 1; j--) {
                            q[j] = (q[j] + 1LL * q[j - 1] * (P - v[i])) % P;
                        }
                    }
                    i64 t = k;
                    std::vector<int> np, nq;
                    while (t > 0) {
                        np.assign(2 * m, 0);
                        nq.assign(2 * m + 1, 0);
                        for (int i = 0; i < m; i++) {
                            for (int j = 0; j <= m; j++) {
                                np[i + j] = (np[i + j] + 1LL * p[i] * (j % 2 ? P - q[j] : q[j])) % P;
                            }
                        }
                        for (int i = 0; i <= m; i++) {
                            for (int j = 0; j <= m; j++) {
                                nq[i + j] = (nq[i + j] + 1LL * q[i] * (j % 2 ? P - q[j] : q[j])) % P;
                            }
                        }
                        for (int i = 0; i < m; i++) {
                            p[i] = np[2 * i + t % 2];
                        }
                        for (int i = 0; i <= m; i++) {
                            q[i] = nq[2 * i];
                        }
                        t /= 2;
                    }
                    ans = (ans + 1LL * res * p[0]) % P;
                    v.pop_back();
                    return;
                }
                for (int y = x + x; y <= s * siz; y += x) {
                    if (s * siz % y == 0) {
                        self(self, y);
                    }
                }
                v.pop_back();
            };
            dfs(dfs, siz);
            
            f[s] = ans;
        }
        
        std::vector<int> dp(n + 1);
        dp[0] = 1;
        
        for (int s = 1; s <= n; s++) {
            for (int i = n; i >= s; i--) {
                int p = 1;
                for (int c = 1; c * s <= i; c++) {
                    p = 1LL * p * f[s] % P;
                    dp[i] = (dp[i] + 1LL * p * invfac[c] % P * dp[i - c * s]) % P;
                }
            }
        }
        
        ans = 1LL * ans * dp[n] % P * fac[n] % P;
    }
    std::cout << ans << "\n";
    
    return 0;
}