#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <cstring>
#include <functional>
#include <iomanip>
#include <iostream>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <vector>
#include <climits>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
using namespace std;
using namespace __gnu_pbds;

using ll = long long;
using ul = unsigned long long;
using db = long double;
using pi = pair<int, int>;
using vi = vector<int>;
using vl = vector<ll>;
using vpi = vector<pi>;
#define mp make_pair
#define pb push_back
#define eb emplace_back
#define x first
#define y second

template<class T> using V = vector<T>; 
template<class T, size_t SZ> using AR = array<T,SZ>; 

#define FOR(i,a,b) for (int i = (a); i < (b); ++i)
#define F0R(i,a) FOR(i,0,a)
#define ROF(i,a,b) for (int i = (b)-1; i >= (a); --i)
#define R0F(i,a) ROF(i,0,a)
#define each(a,x) for (auto& a: x)

#define sz(x) int((x).size())
#define all(x) (x).begin(), (x).end()
#define rep(i,a,b) for(int i = (a); i < (b); i++)
#define per(i,a,b) for(int i = (b) - 1; i >= (a); i--)

#ifdef LOCAL 
template<class A, class B> auto& operator<<(auto &o, pair<A, B> p) { return o << '(' << p.x << ", " << p.y << ')'; }
auto& operator<<(auto& o, auto a) {
    o << "{";
    for (auto b : a) o << b << ", ";
    return o << "}";
}
void dump(auto... x) { ((cerr << x << ", "), ...) << "\n"; }
#define debug(x...) cerr << "[" #x "]: ", dump(x)
#else
#define debug(...) ;
#endif

template<class T> bool ckmin(T& a, const T& b) { return b < a ? a = b, 1 : 0; }
template<class T> bool ckmax(T& a, const T& b) { return a < b ? a = b, 1 : 0; }
template<class T> int lwb(V<T>& a, const T& b) { return int(lb(all(a),b)-bg(a)); }
template<class T> int upb(V<T>& a, const T& b) { return int(ub(all(a),b)-bg(a)); }
template<class T> void remDup(vector<T>& v) { sort(all(v)); v.erase(unique(all(v)),end(v)); }

template<int MOD, int RT> struct mint {
    static const int mod = MOD;
    static constexpr mint rt() { return RT; }
    int v; explicit operator int() const { return v; }
    mint():v(0) {}
    mint(ll _v) { v = int((-MOD < _v && _v < MOD) ? _v : _v % MOD);
        if (v < 0) v += MOD; }
    bool operator==(const mint& o) const {
        return v == o.v; }
    friend bool operator!=(const mint& a, const mint& b) { 
        return !(a == b); }
    friend bool operator<(const mint& a, const mint& b) { 
        return a.v < b.v; }
    friend ostream& operator<<(ostream& os, mint& a) { return os << a.v; }
    friend istream& operator>>(istream& os, mint& a) { ll x; os >> x; a = mint(x); return os; }

    mint& operator+=(const mint& o) { 
        if ((v += o.v) >= MOD) v -= MOD; 
        return *this; }
    mint& operator-=(const mint& o) { 
        if ((v -= o.v) < 0) v += MOD; 
        return *this; }
    mint& operator*=(const mint& o) { 
        v = int((ll)v*o.v%MOD); return *this; }
    mint& operator/=(const mint& o) { return (*this) *= inv(o); }
    friend mint pow(mint a, ll p) {
        mint ans = 1; assert(p >= 0);
        for (; p; p /= 2, a *= a) if (p&1) ans *= a;
        return ans; }
    friend mint inv(const mint& a) { assert(a.v != 0); 
        return pow(a,MOD-2); }
        
    mint operator-() const { return mint(-v); }
    mint& operator++() { return *this += 1; }
    mint& operator--() { return *this -= 1; }
    friend mint operator+(mint a, const mint& b) { return a += b; }
    friend mint operator-(mint a, const mint& b) { return a -= b; }
    friend mint operator*(mint a, const mint& b) { return a *= b; }
    friend mint operator/(mint a, const mint& b) { return a /= b; }
};

constexpr int MOD = 1e9+7;
using mi = mint<MOD,5>;
using vmi = V<mi>;

// based on https://codeforces.com/gym/421334/submission/210707935
using P = vi;
P operator*(const P &a, const P &b) {
    P c(sz(a));
    F0R(i, sz(a)) c[i] = a[b[i]];
    return c;
}
P inv(const P &a) {
    P res(sz(a));
    F0R(i, sz(a)) res[a[i]] = i;
    return res;
}

struct PermSubgroup {
    vector<P> orbit, S, orbit_inv;
};

struct PermGroup {
    int N;
    vector<PermSubgroup> G;
    PermGroup(int _N) : N(_N) {
        G.resize(N);
        P id(N);
        iota(all(id), 0);
        F0R(i, N) {
            G[i].orbit.resize(i + 1);
            G[i].orbit.back() = id;
            G[i].orbit_inv.resize(i + 1);
            G[i].orbit_inv.back() = id;
        }
    }
    bool contains(int i, const P &p) {
        if (i == 0) return true;
        assert(0 <= p[i] && p[i] <= i);
        if (!sz(G[i].orbit[p[i]])) return false;
        return contains(i - 1, G[i].orbit_inv[p[i]] * p);
    }
    void insert(int i, const P &p) {
        if (contains(i, p)) return;
        assert(i > 0);
        G[i].S.pb(p);
        assert(sz(G[i].orbit) == i + 1);
        F0R(j, i + 1)
        if (sz(G[i].orbit[j])) add_to_orbit(i, p * G[i].orbit[j]);
    }
    void add_to_orbit(int i, const P &p) {
        assert(0 <= p[i] && p[i] <= i);
        if (sz(G[i].orbit[p[i]]))
            insert(i - 1, G[i].orbit_inv[p[i]] * p);
        else {
            G[i].orbit[p[i]] = p;
            G[i].orbit_inv[p[i]] = inv(p);
            for (const P &m : G[i].S) add_to_orbit(i, m * p);
        }
    }
    mi size() {
        mi b{1};
        F0R(i, N) {
            int mul = 0;
            F0R(j, i + 1) mul += !!sz(G[i].orbit[j]);
            b *= mul;
        }
        return b;
    }
};

const int N = 3000;
mi dp[2][N][N];
vector<vi> gen;
int n, k;

void solve_schreier_sims() {
    PermGroup group(n);

    rep(i,0,k) {
        group.insert(n - 1, gen[i]);
    }

    mi ord = group.size();

    for (const auto &p : group.G[n-1].orbit_inv) {
        if (p.empty()) continue;
        rep(a,0,n) rep(b,a+1,n) {
            dp[(n-1)%2][p[a]][p[b]] += 1;
        }
    }

    per(i,0,n-1) {
        int cur = i%2;
        int old = (i+1)%2;
        rep(j,0,n) memcpy(dp[cur][j], dp[old][j], n * sizeof(mi));
        for (const auto &p : group.G[i].orbit) {
            if (p.empty()) continue;
            bool any = false;
            rep(l,0,n) {
                if (p[l] != l) {
                    any = true;
                    break;
                }
            }
            if (!any) continue;
            rep(a,0,n) rep(b,0,n) {
                dp[cur][a][b] += dp[old][p[a]][p[b]];
            }
        }
    }

    mi res = 0;
    rep(i,0,n) rep(j,i+1,n) {
        res += dp[0][j][i];
    }

    res /= ord;
    cout << res << '\n';
}

P p;

template<class T>
using ordered_set = tree<T, null_type, greater<T>, rb_tree_tag, tree_order_statistics_node_update>;

vector<vi> get_cycles() {
    vector<vi> res;
    vector<bool> vis(n, false);
    rep(i,0,n) {
        if (vis[i]) continue;
        res.eb();
        int j = i;
        do {
            res.back().pb(j);
            vis[j] = true;
            j = p[j];
        } while (j != i);
    }
    return res;
}

mi cycle_inversions_1(const vi& cycle) {
    int m = sz(cycle);
    mi cur = 0;
    rep(i,0,m) {
        vpi fragment;
        rep(j,0,m) {
            int pos = cycle[j % m];
            int val = p[cycle[(i + j) % m]];
            fragment.eb(pos, val);
        }
        sort(all(fragment));
        ordered_set<int> oset;
        for (auto [_, v] : fragment) {
            cur += oset.order_of_key(v);
            oset.insert(v);
        }
    }
    cur /= mi(m);
    return cur;
}

mi cycle_inversions_2(const vi& c1, const vi& c2) {
    mi res = 0;
    auto d = __gcd(sz(c1), sz(c2));
    mi len = sz(c1) * sz(c2) / d;

    rep(off,0,d) {
        mi sorted_pos = 0;
        mi sorted_val = 0;

        rep(idx,0,len) {
            int i = idx % sz(c1);
            int j = (off + idx) % sz(c2);
            if (c1[i] < c2[j]) sorted_pos += 1;
            if (p[c1[i]] < p[c2[j]]) sorted_val += 1;
        }

        res += sorted_pos * (len - sorted_val);
        res += sorted_val * (len - sorted_pos);
    }

    res /= len;
    return res;
}

void solve1() {
    p = gen[0];
    mi res = 0;
    
    auto cycles = get_cycles();
    for (const auto &cycle : cycles) {
        auto delta = cycle_inversions_1(cycle);
        res += delta;
    }

    rep(i,0,sz(cycles)) {
        rep(j,i+1,sz(cycles)) {
            const auto &c1 = cycles[i];
            const auto &c2 = cycles[j];
            auto delta = cycle_inversions_2(c1, c2);
            res += delta;
        }
    }

    cout << res << '\n';
}

signed main() {
    ios_base::sync_with_stdio(false); cin.tie(nullptr);

    cin >> n >> k;
    
    gen.resize(k, vi(n));
    rep(i,0,k) rep(j,0,n) {
        cin >> gen[i][j];
        gen[i][j]--;
    }

    if (k == 1) {
        solve1();
    }
    else {
        solve_schreier_sims();
    }

    return 0;
}