任意模数下的积和式并没有多项式解法。考虑积和式的组合意义，$i,A_i$ 都是互不相同的，可以看作是二分图的一个完美匹配。因此任意矩阵的积和式就是 $A_{x,y}$ 看作左侧 $x$ 右侧 $y$ 之间的边权后所有完美匹配边权的积的和。

如果我们给所有 $A_{x,y} \neq 1$ 拆成 $1,A_{x,y} - 1$ 这两个东西，然后所有边上都有了一个长度为 $1$ 的边。于是现在就可以不考虑所有长度为 $1$ 的边了。我们求出除开边权为 $1$ 的边后的任意一个二分图匹配，那么这个匹配的贡献直接乘上未匹配点的个数的阶乘即可。因为剩下的点都可以任意匹配，而且一定可以通过 $1$ 边匹配上。

但是这样点的数量还是 $O(2k)$ 的。我们可以考虑对每个连通块分别计算，然后乘法原理乘起来即可。这样每个连通块内的点数就一定不超过 $k$ 了。一种很容易想到的做法是对每个连通块做状压 $dp$ ，考虑把点数较小的一边压成状态，然后枚举点数较多的一边来做 $dp$ ，这样做一个点数为 $n$ 边数为 $m$ 的连通块的复杂度是 $O((\frac n 2+m)2^{\frac n 2})$ 。这样是过不去的，如果卡满就是 $50 \times 2^{25}$ 。

考虑另外一种做法，求出这个二分图的一棵生成树，然后对于所有环边枚举环边的选取情况，然后在树上做一次背包来解决。这样做的复杂度是 $O(2^{m-n} n^2)$ 。

我们选取其中较小的一种来做，复杂度是 $O(\min\{m 2^{\frac n 2},n^22^{m-n}\})$ 。在 $n < \frac 2 3 m$ 的时候用第一个做法，复杂度就是 $O(n^2 2^{\frac 1 3 m})$ 。

于是最终复杂度就是 $O(k^22^{\frac 1 3 k})$ ，而且非常不满。

cpp

#include "iostream"
#include "algorithm"
#include "cstring"
#include "cstdio"
#include "cmath"
#include "vector"
#include "map"
#include "set"
#include "queue"
using namespace std;
#define MAXN 200006
//#define int long long
#define rep(i, a, b) for (int i = (a), i##end = (b); i <= i##end; ++i)
#define per(i, a, b) for (int i = (a), i##end = (b); i >= i##end; --i)
#define pii pair<int,int>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define eb emplace_back
#define vi vector<int>
#define all(x) (x).begin() , (x).end()
#define mem( a ) memset( a , 0 , sizeof a )
typedef long long ll;
const int P = 1e9 + 7;
int n , k;
int A[MAXN];
int x[MAXN] , y[MAXN] , w[MAXN];
vi X , Y;
int fa[MAXN] , se[MAXN] , sz[MAXN];
int find( int x ) {
	return x == fa[x] ? x : fa[x] = find( fa[x] );
}

inline void Inc( int& x , int y ) {
	x = ( x + y > P ? x + y - P : x + y );
}

vector<pii> G[MAXN];

vi nl , nr;
int bc[MAXN];
int as[106];
namespace lsjtxdy{
	int dp[106][1 << 16 | 66];
	void wkr( ) {
		int sl = nl.size() , sr = nr.size();
		rep( i , 0 , sl + sr ) rep( j , 0 , ( 1 << sl ) - 1 ) dp[i][j] = 0;
		rep( i , 0 , sl + sr ) as[i] = 0;
		dp[0][0] = 1;
		rep( i , 1 , sr ) rep( s , 0 , ( 1 << sl ) - 1 ) {
			Inc( dp[i][s] , dp[i - 1][s] );
			int u = nr[i - 1];
			for( auto t : G[u] ) if( ~s & ( 1 << bc[t.fi] ) ) {
				int v = t.fi , w = t.se;
				int ts = s | ( 1 << bc[v] );
				Inc( dp[i][ts] , dp[i - 1][s] * 1ll * w % P );
			}
		}
		rep( i , 0 , ( 1 << sl ) - 1 ) Inc( as[__builtin_popcount( i )] , dp[sr][i] );
	}
}
namespace zjktsdy {
	vector<pair<int,pii>> fe;
	int dep[MAXN];
	vector<pii> g[MAXN];
	vi nd;
	void dfs( int u , int f ) {
		nd.pb( u );
		dep[u] = dep[f] + 1;
		for( auto t : G[u] ) if( t.fi != f ) {
			if( !dep[t.fi] ) dfs( t.fi , u ) , g[u].pb( t );
			else if( dep[t.fi] > dep[u] ) fe.eb( mp( u , t ) );
		}
	}
	int fk[MAXN];
	int dp[106][2][52] , pd[2][52] , siz[106];
	void afs( int u ) {
		int sz = 1;
		mem( dp[u] );
		dp[u][0][0] = 1;
		for( auto t : g[u] ) {
			int v = t.fi;
			afs( v );
			memcpy( pd , dp[u] , sizeof dp[u] );
			mem( dp[u] );
			rep( i , 0 , sz ) rep( j , 0 , siz[v] ) {
				Inc( dp[u][0][i + j] , pd[0][i] * 1ll * ( dp[v][0][j] + dp[v][1][j] ) % P );
				if( !fk[u] ) {
					Inc( dp[u][1][i + j] , pd[1][i] * 1ll * ( dp[v][0][j] + dp[v][1][j] ) % P );
					if( !fk[v] ) Inc( dp[u][1][i + j + 1] , pd[0][i] * 1ll * dp[v][0][j] % P * t.se % P );
				}
			}
			sz += siz[v];
		}
		siz[u] = sz;
	}
	
	void wkr( int r ) {
		nd.clear();
		rep( i , 0 , k * 2 ) as[i] = dep[i] = 0 , g[i].clear();
		fe.clear();
		dfs( r , r );
		int ta = 0 , ts = 1;
		for( int s = 0 ; s < ( 1 << fe.size() ) ; ++ s ) {
			for( int x : nd ) fk[x] = 0;
			ta = __builtin_popcount( s ) , ts = 1;
			rep( i , 0 , fe.size() - 1 ) if( s & ( 1 << i ) ) {
				int u = fe[i].fi , v = fe[i].se.fi;
				if( fk[u] || fk[v] ) { ts = 0; break; }
				ts = ts * 1ll * fe[i].se.se % P , fk[u] = fk[v] = 1;
			}
			if( !ts ) continue;
			afs( r );
			rep( i , 0 , k ) Inc( as[i + ta] , ( dp[r][0][i] + dp[r][1][i] ) * 1ll * ts % P );
		}
	}
}

int ans[105] , pas[105];
int J[100006];
void solve() {
	cin >> n >> k;
	J[0] = 1;
	rep( i , 1 , n ) J[i] = J[i - 1] * 1ll * i % P;
	rep( i , 1 , k ) {
		scanf("%d%d%d",x + i,y + i,w + i);
		-- w[i];
		X.pb( x[i] ) , Y.pb( y[i] );
	}
	sort( all( X ) ) , sort( all( Y ) );
	X.erase( unique( all( X ) ) , X.end() ) , Y.erase( unique( all( Y ) ) , Y.end() );
	rep( i , 1 , k ) 
		x[i] = lower_bound( all( X ) , x[i] ) - X.begin() + 1 , y[i] = lower_bound( all( Y ) , y[i] ) - Y.begin() + 1;
	int sx = X.size() , sy = Y.size();
	rep( i , 1 , k ) {
		G[x[i]].eb( mp( y[i] + sx , w[i] ) );
		G[y[i] + sx].eb( mp( x[i] , w[i] ) );
//		cout << x[i] << ' ' << y[i] + sx << endl;
	}
	rep( i , 1 , sx + sy ) sz[i] = 1 , fa[i] = i;
	rep( i , 1 , k ) {
		int u = find( x[i] ) , v = find( y[i] + sx );
		++ se[v];
		if( u == v ) continue;
		fa[u] = v , sz[v] += sz[u] , se[v] += se[u];
	}
	ans[0] = 1;
	rep( i , 1 , sx + sy ) if( fa[i] == i && sz[i] > 1 ) {
		if( sz[i] <= 2 * se[i] / 3 ) {
			nl.clear() , nr.clear();
			rep( j , 1 , sx ) if( find( j ) == i ) nl.pb( j ) , bc[j] = nl.size() - 1;
			rep( j , sx + 1 , sx + sy ) if( find( j ) == i ) nr.pb( j ) , bc[j] = nr.size() - 1;
			if( nl.size() > nr.size() ) swap( nl , nr );
			lsjtxdy::wkr( );
		} else {
			zjktsdy::wkr( i );
		}
		memcpy( pas , ans , sizeof ans );
		mem( ans );
		rep( i , 0 , k ) if( as[i] ) 
			rep( j , i , k )
				Inc( ans[j] , as[i] * 1ll * pas[j - i] % P );
	}
	int res = 0;
	rep( i , 0 , k ) Inc( res , ans[i] * 1ll * J[n - i] % P );
	cout << res << endl;
}

signed main() {
//    int T;cin >> T;while( T-- ) solve();
    solve();
}