#include <bits/stdc++.h>
using i64 = long long;
using namespace std;
const int N = 5e5 + 5;
// const int B = 400;
// const int M = 4e6 + 5;
// const int base = 15171;
// const int mod = 998244353;
// const i64 mod = (i64)1e18 + 7;
// const double pi = acos(-1);

const int P = 998244353;
using i64 = long long;
using Poly = vector<int>;
/*---------------------------------------------------------------------------*/
#define MUL(a, b) (i64(a) * (b) % P)
#define ADD(a, b) (((a) += (b)) >= P ? (a) -= P : 0) // (a += b) %= P
#define SUB(a, b) (((a) -= (b)) < 0 ? (a) += P : 0)  // ((a -= b) += P) %= P
Poly getInv(int L)
{
    Poly inv(L);
    inv[1] = 1;
    for (int i = 2; i < L; i++)
        inv[i] = MUL((P - P / i), inv[P % i]);
    return inv;
}
int POW(i64 a, int b = P - 2, i64 x = 1)
{
    for (; b; b >>= 1, a = a * a % P)
        if (b & 1)
            x = x * a % P;
    return x;
}
auto inv = getInv(N);
namespace NTT
{
    const int g = 3;
    Poly Omega(int L)
    {
        int wn = POW(g, P / L);
        Poly w(L);
        w[L >> 1] = 1;
        for (int i = L / 2 + 1; i < L; i++)
            w[i] = MUL(w[i - 1], wn);
        for (int i = L / 2 - 1; i >= 1; i--)
            w[i] = w[i << 1];
        return w;
    }
    auto W = Omega(1 << 20); // Length
    void DIF(int *a, int n)
    {
        for (int k = n >> 1; k; k >>= 1)
            for (int i = 0, y; i < n; i += k << 1)
                for (int j = 0; j < k; ++j)
                    y = a[i + j + k], a[i + j + k] = MUL(a[i + j] - y + P, W[k + j]), ADD(a[i + j], y);
    }
    void IDIT(int *a, int n)
    {
        for (int k = 1; k < n; k <<= 1)
            for (int i = 0, x, y; i < n; i += k << 1)
                for (int j = 0; j < k; ++j)
                    x = a[i + j], y = MUL(a[i + j + k], W[k + j]),
                    a[i + j + k] = x - y < 0 ? x - y + P : x - y, ADD(a[i + j], y);
        int Inv = P - (P - 1) / n;
        for (int i = 0; i < n; i++)
            a[i] = MUL(a[i], Inv);
        reverse(a + 1, a + n);
    }
}
/*---------------------------------------------------------------------------*/
namespace Polynomial
{
    // basic operator
    int norm(int n) { return 1 << (__lg(n - 1) + 1); }
    void norm(Poly &a)
    {
        if (!a.empty())
            a.resize(norm(a.size()), 0);
        else
            a = {0};
    }
    void DFT(Poly &a) { NTT::DIF(a.data(), a.size()); }
    void IDFT(Poly &a) { NTT::IDIT(a.data(), a.size()); }
    Poly &dot(Poly &a, Poly &b)
    {
        for (int i = 0; i < a.size(); i++)
            a[i] = MUL(a[i], b[i]);
        return a;
    }

    // mul / div int
    Poly &operator*=(Poly &a, int b)
    {
        for (auto &x : a)
            x = MUL(x, b);
        return a;
    }
    Poly operator*(Poly a, int b) { return a *= b; }
    Poly operator*(int a, Poly b) { return b * a; }
    Poly &operator/=(Poly &a, int b) { return a *= POW(b); }
    Poly operator/(Poly a, int b) { return a /= b; }

    // Poly add / sub
    Poly &operator+=(Poly &a, Poly b)
    {
        a.resize(max(a.size(), b.size()));
        for (int i = 0; i < b.size(); i++)
            ADD(a[i], b[i]);
        return a;
    }
    Poly operator+(Poly a, Poly b) { return a += b; }
    Poly &operator-=(Poly &a, Poly b)
    {
        a.resize(max(a.size(), b.size()));
        for (int i = 0; i < b.size(); i++)
            SUB(a[i], b[i]);
        return a;
    }
    Poly operator-(Poly a, Poly b) { return a -= b; }

    // Poly mul
    Poly operator*(Poly a, Poly b)
    {
        int n = a.size() + b.size() - 1, L = norm(n);
        if (a.size() <= 8 || b.size() <= 8)
        {
            Poly c(n);
            for (int i = 0; i < a.size(); i++)
                for (int j = 0; j < b.size(); j++)
                    c[i + j] = (c[i + j] + (i64)a[i] * b[j]) % P;
            return c;
        }
        a.resize(L), b.resize(L);
        DFT(a), DFT(b), dot(a, b), IDFT(a);
        return a.resize(n), a;
    }
}
using namespace Polynomial;
int n, k, fa[N], siz[N], fac[N], ifac[N];
int binom(int x, int y)
{
    if (x < 0 || y < 0 || x < y)
        return 0;
    return (i64)fac[x] * ifac[y] % P * ifac[x - y] % P;
}
int find(int x)
{
    if (x == fa[x])
        return x;
    else
        return fa[x] = find(fa[x]);
}
void merge(int x, int y)
{
    int fx = find(x), fy = find(y);
    if (fx ^ fy)
    {
        fa[fy] = fx;
        siz[fx] += siz[fy];
    }
}
void solve()
{
    cin >> n >> k;
    iota(fa + 1, fa + n + 1, 1);
    fill(siz + 1, siz + n + 1, 1);
    for (int i = 1, x; i <= n; i++)
        cin >> x, merge(i, x);
    vector<Poly> tmp;
    for (int i = 1; i <= n; i++)
    {
        if (find(i) != i)
            continue;
        int cur = siz[i];
        Poly f(min(n / 2 + 1, cur + 1));
        for (int j = 0; j < f.size(); j++)
            f[j] = (binom(cur - j, j) + binom(cur - j - 1, j - 1)) % P;
        tmp.push_back(f);
    }
    auto solve = [&](auto &solve, int l, int r) -> Poly
    {
        if (l == r)
            return tmp[l];
        int mid = l + r >> 1;
        return solve(solve, l, mid) * solve(solve, mid + 1, r);
    };
    auto res = solve(solve, 0, tmp.size() - 1);
    res.resize(n + 1);
    cout << res[k] << '\n';
}
signed main()
{
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);
    fac[0] = 1;
    for (int i = 1; i < N; i++)
        fac[i] = (i64)fac[i - 1] * i % P;
    ifac[N - 1] = POW(fac[N - 1]);
    for (int i = N - 2; i >= 0; i--)
        ifac[i] = (i64)ifac[i + 1] * (i + 1) % P;
    int t = 1;
    cin >> t;
    // cout << fixed << setprecision(10);
    while (t--)
        solve();
}