拆贡献。假设分为了 $k$ 组，我们从每组中选择一个元素，设为 $a_1, \dots, a_k$，在这组方案中这 $k$ 个数的贡献就是 $\prod a_i$。

转为先枚举 $k$ 个元素，再将剩下的分组。前面的部分可以 dp，设 $f_{i, j}$ 表示前 $i$ 个数选了 $j$ 个数的方案，可以 $O(n^2m)$ 算出。

对于后面，我们可以转化成以下问题：

有 $n \times m - k$ 个相同的球，要将它们分到 $k$ 个有标号的盒子，使得每个盒子的球数 $\bmod m = -1$。

显然答案为：

$$(nm - k)! \times [x^n] \left(\sum_{i=1}^n \frac{x^i}{(im - 1)!} \right)^k$$

直接模拟即可做到 $O(n^3)$，无法通过。

设 $F(x) = \sum_{i=1}^n \frac{x^i}{(im - 1)!}$。考虑拉格朗日插值，我们维护一系列点值 $v_i = F(i)$，我们只需要用 $v_i^k$ 反演回去就可以了。但 $F^k(x)$ 的次数是 $O(n^2)$ 的，也就是要维护 $O(n^2)$ 个点值，那复杂度还是 $O(n^3)$。

由于我们只需要保留前 $n$ 项的系数，所以可以考虑快速幂。先求出所有 $F^{2^i}(x)$ 的前 $n$ 项系数，求 $F^k(x)$ 的时候只需要将 $k$ 拆成二进制表示，将对应 $\le \log n$ 个 $F^{2^i}(x)$ 乘起来，这样次数就只有 $O(n\log n)$，只需要维护前 $O(n\log n)$ 个点值。

准确来说我们是要从 $k - 1$ 递推到 $k$（即只反转一段低位），这样才不会再多一个 $\log$。

算所有 $F^{2^i}(x)$ 的前 $n$ 项系数可以直接从小到大递推。注意这里只用维护前 $2n$ 项系数就好了，这部分是 $O(n^2 \log n)$ 的。

但还要求点值，$O(\log n)$ 个多项式，每个有 $O(n\log n)$ 个点值，总共是 $O(n^2\log^2 n)$ 的。暂时没想到优化方案。

那么总复杂度就是 $O(n^2(\log^2 n + m))$，实际运行中三个部分 $O(n^2m), O(n^2\log n), O(n^2\log^2 n)$ 用时占比约为 $4 : 1 : 4$。

暂时未能通过。

#include <algorithm>
#include <iostream>
using namespace std;

#define MAXN 2005
#define MAXM 105
#define MAXS 200005
#define MAXL 11
#define MOD 1000000007ll

using ll = long long;

ll fac[MAXS], inf[MAXS], inv[MAXS];

ll qpow(ll b, ll p)
{
    ll r = 1;
    while (p)
    {
        if (p & 1)
        {
            r = r * b % MOD;
        }
        b = b * b % MOD;
        p >>= 1;
    }
    return r;
}

int Tid;

int n, m;

int a[MAXS];

ll f[MAXN];

template <int S>
struct LI
{
    int N;

    ll f[MAXN];

    ll g[S][MAXN];

    void init()
    {
        fill(f + 1, f + n + 1, 0);
        f[0] = 1;
        for (int i = 1; i <= N; i++)
        {
            for (int j = n; j; j--)
            {
                f[j] = (f[j - 1] + MOD - f[j] * i % MOD) % MOD;
            }
            f[0] = f[0] * (MOD - i) % MOD;
        }
        for (int i = 1; i <= N; i++)
        {
            f[0] = f[0] * (MOD - inv[i]) % MOD;
            for (int j = 1; j <= n; j++)
            {
                f[j] = (MOD - inv[i]) * (f[j] + MOD - f[j - 1]) % MOD;
                g[i][j] = f[j] * inf[N - i] % MOD * inf[i - 1] % MOD;
            }
            g[i][0] = f[0] * inf[N - i] % MOD * inf[i - 1] % MOD;
            if ((N - i) & 1)
            {
                for (int j = 0; j <= n; j++)
                {
                    g[i][j] = MOD - g[i][j];
                }
            }
            // clog << "coef " << i << ": ";
            // for (int j = 0; j <= n; j++)
            // {
            //     clog << ' ' << f[j];
            // }
            // clog << endl;
            for (int j = n; j; j--)
            {
                f[j] = (f[j - 1] + MOD - f[j] * i % MOD) % MOD;
            }
            f[0] = f[0] * (MOD - i) % MOD;
        }
    }

    ll h[MAXN];

    void recov(ll *f)
    {
        fill(h, h + n + 1, 0);
        for (int i = 1; i <= N; i++)
        {
            for (int j = 0; j <= n; j++)
            {
                h[j] = (h[j] + f[i] * g[i][j]) % MOD;
            }
        }
    }

    ll cn(ll *f)
    {
        ll res = 0;
        for (int i = 1; i <= N; i++)
        {
            res = (res + f[i] * g[i][n]) % MOD;
        }
        return res;
    }
};

LI<MAXN << 1> lis;
LI<MAXN * 10> lib;

ll g[MAXN << 1];
ll pw[MAXL][MAXN * 10];
ll h[MAXL][MAXN * 10];
ll ih[MAXL][MAXN * 10];

ll cur[MAXN * 10];

void solve()
{
    cin >> n >> m;
    for (int i = 1; i <= n * m; i++)
    {
        cin >> a[i];
    }
    fill(f + 1, f + n + 1, 0);
    f[0] = 1;
    for (int i = 1; i <= n * m; i++)
    {
        for (int j = n; j; j--)
        {
            f[j] = (f[j] + f[j - 1] * a[i]) % MOD;
        }
    }

    // for (int i = 1; i <= n; i++)
    // {
    //     clog << f[i] << ' ';
    // }
    // clog << endl;

    lis.N = 2 * n + 1;
    lib.N = 10 * n + 1;
    lis.init();
    lib.init();
    for (int i = 1; i <= lib.N; i++)
    {
        pw[0][i] = i;
        h[0][i] = 0;
        for (int j = n; j; j--)
        {
            h[0][i] = (h[0][i] + inf[j * m - 1]) * i % MOD;
        }
        ih[0][i] = qpow(h[0][i], MOD - 2);
        if (i <= lis.N)
        {
            g[i] = h[0][i];
            // clog << g[i] << ' ';
        }
    }
    // clog << endl;
    for (int i = 1; i <= 10; i++)
    {
        for (int j = 1; j <= lis.N; j++)
        {
            g[j] = g[j] * g[j] % MOD;
        }
        lis.recov(g);

        // clog << i << ':';
        // for (int j = 1; j <= lis.N; j++)
        // {
        //     clog << ' ' << g[j];
        // }
        // clog << endl;
        // clog << i << ':';
        // for (int j = 0; j <= n; j++)
        // {
        //     clog << ' ' << lis.h[j];
        // }
        // clog << endl;

        for (int j = 1; j <= lib.N; j++)
        {
            pw[i][j] = pw[i - 1][j] * pw[i - 1][j] % MOD;
            h[i][j] = 0;
            for (int k = n; k >= (1 << i); k--)
            {
                h[i][j] = (h[i][j] * j + lis.h[k]) % MOD;
            }
            h[i][j] = h[i][j] * pw[i][j] % MOD;
            ih[i][j] = qpow(h[i][j], MOD - 2);
            if (j <= lis.N)
            {
                g[j] = h[i][j];
            }
        }
        // clog << i << ':';
        // for (int j = 1; j <= lis.N; j++)
        // {
        //     clog << ' ' << g[j];
        // }
        // clog << endl;
    }
    for (int i = 1; i <= lib.N; i++)
    {
        cur[i] = 1;
    }
    ll res = 0;
    for (int i = 1; i <= n; i++)
    {
        int v = i ^ (i - 1);
        for (int k = 0; (1 << k) <= v; k++)
        {
            for (int j = 1; j <= lib.N; j++)
            {
                cur[j] = (cur[j] * ((1 << k) == v ? h : ih)[k][j]) % MOD;
            }
            v ^= (1 << k);
        }
        res = (res + f[i] * lib.cn(cur) % MOD * fac[n * m - i]) % MOD;
        // clog << lib.cn(cur) % MOD * fac[n * m - i] << ' ';
    }
    // clog << endl;
    cout << res << '\n';
}

int main()
{
    freopen("divide.in", "r", stdin);
    freopen("divide.out", "w", stdout);
    fac[0] = 1;
    for (int i = 1; i < MAXS; i++)
    {
        fac[i] = fac[i - 1] * i % MOD;
    }
    inf[MAXS - 1] = qpow(fac[MAXS - 1], MOD - 2);
    for (int i = MAXS - 1; i; i--)
    {
        inf[i - 1] = inf[i] * i % MOD;
        inv[i] = inf[i] * fac[i - 1] % MOD;
    }
    int t;
    cin >> Tid >> t;
    while (t--)
    {
        solve();
    }
    return 0;
}