Monocarp is playing a game "Assimilation IV". In this game he manages a great empire: builds cities and conquers new lands.

Monocarp's empire has $n$ cities. In order to conquer new lands he plans to build one Monument in each city. The game is turn-based and, since Monocarp is still amateur, he builds exactly one Monument per turn.

Monocarp has $m$ points on the map he'd like to control using the constructed Monuments. For each point he knows the distance between it and each city. Monuments work in the following way: when built in some city, a Monument controls all points at distance at most $1$ to this city. Next turn, the Monument controls all points at distance at most $2$, the turn after — at distance at most $3$, and so on. Monocarp will build $n$ Monuments in $n$ turns and his empire will conquer all points that are controlled by at least one Monument.

Monocarp can't figure out any strategy, so during each turn he will choose a city for a Monument randomly among all remaining cities (cities without Monuments). Monocarp wants to know how many points (among $m$ of them) he will conquer at the end of turn number $n$. Help him to calculate the expected number of conquered points!


-----Input-----

The first line contains two integers $n$ and $m$ ($1 \le n \le 20$; $1 \le m \le 5 \cdot 10^4$) — the number of cities and the number of points.

Next $n$ lines contains $m$ integers each: the $j$-th integer of the $i$-th line $d_{i, j}$ ($1 \le d_{i, j} \le n + 1$) is the distance between the $i$-th city and the $j$-th point.


-----Output-----

It can be shown that the expected number of points Monocarp conquers at the end of the $n$-th turn can be represented as an irreducible fraction $\frac{x}{y}$. Print this fraction modulo $998\,244\,353$, i. e. value $x \cdot y^{-1} mod 998244353$ where $y^{-1}$ is such number that $y \cdot y^{-1} mod 998244353 = 1$.


-----Examples-----

Input
3 5
1 4 4 3 4
1 4 1 4 2
1 4 4 4 3
Output
166374062


-----Note-----

Let's look at all possible orders of cities Monuments will be build in:

$[1, 2, 3]$:

the first city controls all points at distance at most $3$, in other words, points $1$ and $4$;

the second city controls all points at distance at most $2$, or points $1$, $3$ and $5$;

the third city controls all points at distance at most $1$, or point $1$.

In total, $4$ points are controlled.

$[1, 3, 2]$: the first city controls points $1$ and $4$; the second city — points $1$ and $3$; the third city — point $1$. In total, $3$ points.

$[2, 1, 3]$: the first city controls point $1$; the second city — points $1$, $3$ and $5$; the third city — point $1$. In total, $3$ points.

$[2, 3, 1]$: the first city controls point $1$; the second city — points $1$, $3$ and $5$; the third city — point $1$. In total, $3$ points.

$[3, 1, 2]$: the first city controls point $1$; the second city — points $1$ and $3$; the third city — points $1$ and $5$. In total, $3$ points.

$[3, 2, 1]$: the first city controls point $1$; the second city — points $1$, $3$ and $5$; the third city — points $1$ and $5$. In total, $3$ points.

The expected number of controlled points is $\frac{4 + 3 + 3 + 3 + 3 + 3}{6}$ $=$ $\frac{19}{6}$ or $19 \cdot 6^{-1}$ $\equiv$ $19 \cdot 166374059$ $\equiv$ $166374062$ $\pmod{998244353}$

Source: codeforces | Difficulty: HARD
Tags: Probability

Combinatorics

Amortized analysis

Dynamic programming

Mathematics


Solutions:
ans = 0
(n, m) = map(int, input().split())
mod = 998244353
a = [list(map(int, input().split())) for i in range(n)]
fac = 1
for i in range(1, n + 1):
	fac *= i
inv = pow(fac, mod - 2, mod)
for j in range(m):
	na = sorted([a[i][j] for i in range(n)])
	now = 1
	able = 0
	for i in range(n):
		while len(na) > 0 and na[-1] > n - i:
			del na[-1]
			able += 1
		now *= able
		able -= 1
	ans += now * inv
	ans %= mod
print((m - ans) % mod)


import sys
import math
from math import factorial, inf, gcd, sqrt
from heapq import *
from functools import *
from itertools import *
from collections import *
from typing import *
from bisect import *
import random
sys.setrecursionlimit(10 ** 5)

def rarray():
	return [int(i) for i in input().split()]
(n, m) = rarray()
mod = 998244353
a = [rarray() for i in range(n)]
a = list(zip(*a))
fac = [[0] * (n + 1) for i in range(n + 1)]
fac[0][0] = 1
for i in range(1, n + 1):
	fac[i][0] = 1
	for j in range(1, i + 1):
		fac[i][j] = fac[i][j - 1] * (i - j + 1) % mod
ans = 0
for i in range(m):
	a[i] = sorted(a[i])
	j = 0
	res = 1
	while j < n:
		k = j
		while j < n and a[i][j] == a[i][k]:
			j += 1
		t = j - k
		res = res * fac[a[i][k] - 1 - k][t] % mod
	ans = (ans + fac[n][n] - res + mod) % mod
print(ans * pow(fac[n][n], mod - 2, mod) % mod)


from bisect import bisect, bisect_left
from collections import *
from heapq import *
from math import gcd, ceil, sqrt, floor, inf
from itertools import *
from operator import add, mul, sub, xor, truediv, floordiv
from functools import *
import os
import sys
from io import BytesIO, IOBase
BUFSIZE = 8192

class FastIO(IOBase):
	newlines = 0

	def __init__(self, file):
		self._fd = file.fileno()
		self.buffer = BytesIO()
		self.writable = 'x' in file.mode or 'r' not in file.mode
		self.write = self.buffer.write if self.writable else None

	def read(self):
		while True:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			if not b:
				break
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines = 0
		return self.buffer.read()

	def readline(self):
		while self.newlines == 0:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			self.newlines = b.count(b'\n') + (not b)
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines -= 1
		return self.buffer.readline()

	def flush(self):
		if self.writable:
			os.write(self._fd, self.buffer.getvalue())
			(self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(IOBase):

	def __init__(self, file):
		self.buffer = FastIO(file)
		self.flush = self.buffer.flush
		self.writable = self.buffer.writable
		self.write = lambda s: self.buffer.write(s.encode('ascii'))
		self.read = lambda : self.buffer.read().decode('ascii')
		self.readline = lambda : self.buffer.readline().decode('ascii')
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')

def RL():
	return map(int, sys.stdin.readline().split())

def RLL():
	return list(map(int, sys.stdin.readline().split()))

def N():
	return int(input())

def A(n):
	return [0] * n

def AI(n, x):
	return [x] * n

def A2(n, m):
	return [[0] * m for i in range(n)]

def G(n):
	return [[] for i in range(n)]

def GP(it):
	return [[ch, len(list(g))] for (ch, g) in groupby(it)]
from types import GeneratorType

def bootstrap(f, stack=[]):

	def wrappedfunc(*args, **kwargs):
		if stack:
			return f(*args, **kwargs)
		else:
			to = f(*args, **kwargs)
			while True:
				if type(to) is GeneratorType:
					stack.append(to)
					to = next(to)
				else:
					stack.pop()
					if not stack:
						break
					to = stack[-1].send(to)
			return to
	return wrappedfunc
mod = 10 ** 9 + 7
farr = [1]
ifa = []

def fact(x, mod=0):
	if mod:
		while x >= len(farr):
			farr.append(farr[-1] * len(farr) % mod)
	else:
		while x >= len(farr):
			farr.append(farr[-1] * len(farr))
	return farr[x]

def ifact(x, mod):
	global ifa
	fact(x, mod)
	ifa.append(pow(farr[-1], mod - 2, mod))
	for i in range(x, 0, -1):
		ifa.append(ifa[-1] * i % mod)
	ifa.reverse()

def per(i, j, mod=0):
	if i < j:
		return 0
	if not mod:
		return fact(i) // fact(i - j)
	return farr[i] * ifa[i - j] % mod

def com(i, j, mod=0):
	if i < j:
		return 0
	if not mod:
		return per(i, j) // fact(j)
	return per(i, j, mod) * ifa[j] % mod

def catalan(n):
	return com(2 * n, n) // (n + 1)

def isprime(n):
	for i in range(2, int(n ** 0.5) + 1):
		if n % i == 0:
			return False
	return True

def floorsum(a, b, c, n):
	if a == 0:
		return b // c * (n + 1)
	if a >= c or b >= c:
		return floorsum(a % c, b % c, c, n) + b // c * (n + 1) + a // c * n * (n + 1) // 2
	m = (a * n + b) // c
	return n * m - floorsum(c, c - b - 1, a, m - 1)

def inverse(a, m):
	a %= m
	if a <= 1:
		return a
	return (1 - inverse(m, a) * m) // a % m

def lowbit(n):
	return n & -n

class BIT:

	def __init__(self, arr):
		self.arr = arr
		self.n = len(arr) - 1

	def update(self, x, v):
		while x <= self.n:
			self.arr[x] += v
			x += x & -x

	def query(self, x):
		ans = 0
		while x:
			ans += self.arr[x]
			x &= x - 1
		return ans

class ST:

	def __init__(self, arr):
		n = len(arr)
		mx = n.bit_length()
		self.st = [[0] * mx for i in range(n)]
		for i in range(n):
			self.st[i][0] = arr[i]
		for j in range(1, mx):
			for i in range(n - (1 << j) + 1):
				self.st[i][j] = max(self.st[i][j - 1], self.st[i + (1 << j - 1)][j - 1])

	def query(self, l, r):
		if l > r:
			return -inf
		s = (r + 1 - l).bit_length() - 1
		return max(self.st[l][s], self.st[r - (1 << s) + 1][s])

class DSU:

	def __init__(self, n):
		self.c = [-1] * n

	def same(self, x, y):
		return self.find(x) == self.find(y)

	def find(self, x):
		if self.c[x] < 0:
			return x
		self.c[x] = self.find(self.c[x])
		return self.c[x]

	def union(self, u, v):
		(u, v) = (self.find(u), self.find(v))
		if u == v:
			return False
		if self.c[u] > self.c[v]:
			(u, v) = (v, u)
		self.c[u] += self.c[v]
		self.c[v] = u
		return True

	def size(self, x):
		return -self.c[self.find(x)]

class UFS:

	def __init__(self, n):
		self.parent = [i for i in range(n)]
		self.ranks = [0] * n

	def find(self, x):
		if x != self.parent[x]:
			self.parent[x] = self.find(self.parent[x])
		return self.parent[x]

	def union(self, u, v):
		(pu, pv) = (self.find(u), self.find(v))
		if pu == pv:
			return False
		if self.ranks[pu] >= self.ranks[pv]:
			self.parent[pv] = pu
			if self.ranks[pv] == self.ranks[pu]:
				self.ranks[pu] += 1
		else:
			self.parent[pu] = pv

class UF:

	def __init__(self, n):
		self.parent = [i for i in range(n)]
		self.ranks = [0] * n
		self.size = AI(n, 1)
		self.edge = A(n)

	def find(self, x):
		if x != self.parent[x]:
			self.parent[x] = self.find(self.parent[x])
		return self.parent[x]

	def union(self, u, v):
		(pu, pv) = (self.find(u), self.find(v))
		if pu == pv:
			self.edge[pu] += 1
			return False
		if self.ranks[pu] >= self.ranks[pv]:
			self.parent[pv] = pu
			self.edge[pu] += self.edge[pv] + 1
			self.size[pu] += self.size[pv]
			if self.ranks[pv] == self.ranks[pu]:
				self.ranks[pu] += 1
		else:
			self.parent[pu] = pv
			self.edge[pv] += self.edge[pu] + 1
			self.size[pv] += self.size[pu]

def Prime(n):
	c = 0
	prime = []
	flag = [0] * (n + 1)
	for i in range(2, n + 1):
		if not flag[i]:
			prime.append(i)
			c += 1
		for j in range(c):
			if i * prime[j] > n:
				break
			flag[i * prime[j]] = prime[j]
			if i % prime[j] == 0:
				break
	return prime

def dij(s, graph):
	d = AI(n, inf)
	d[s] = 0
	heap = [(0, s)]
	vis = A(n)
	while heap:
		(dis, u) = heappop(heap)
		if vis[u]:
			continue
		vis[u] = 1
		for (v, w) in graph[u]:
			if d[v] > d[u] + w:
				d[v] = d[u] + w
				heappush(heap, (d[v], v))
	return d

def bell(s, g):
	dis = AI(n, inf)
	dis[s] = 0
	for i in range(n - 1):
		for (u, v, w) in edge:
			if dis[v] > dis[u] + w:
				dis[v] = dis[u] + w
	change = A(n)
	for i in range(n):
		for (u, v, w) in edge:
			if dis[v] > dis[u] + w:
				dis[v] = dis[u] + w
				change[v] = 1
	return dis

def lcm(a, b):
	return a * b // gcd(a, b)

def lis(nums):
	res = []
	for k in nums:
		i = bisect.bisect_left(res, k)
		if i == len(res):
			res.append(k)
		else:
			res[i] = k
	return len(res)

def RP(nums):
	n = len(nums)
	s = set(nums)
	d = {}
	for (i, k) in enumerate(sorted(s), 1):
		d[k] = i
	bi = BIT([0] * (len(s) + 1))
	ans = 0
	for i in range(n - 1, -1, -1):
		ans += bi.query(d[nums[i]] - 1)
		bi.update(d[nums[i]], 1)
	return ans

class DLN:

	def __init__(self, val):
		self.val = val
		self.pre = None
		self.next = None

def nb(i, j, n, m):
	for (ni, nj) in [[i + 1, j], [i - 1, j], [i, j - 1], [i, j + 1]]:
		if 0 <= ni < n and 0 <= nj < m:
			yield (ni, nj)

def topo(n):
	q = deque()
	res = []
	for i in range(1, n + 1):
		if ind[i] == 0:
			q.append(i)
			res.append(i)
	while q:
		u = q.popleft()
		for v in g[u]:
			ind[v] -= 1
			if ind[v] == 0:
				q.append(v)
				res.append(v)
	return res

@bootstrap
def gdfs(r, p):
	for ch in g[r]:
		if ch != p:
			yield gdfs(ch, r)
	yield None

@bootstrap
def dfs(r, p):
	for (ch, w) in g[r]:
		if ch != p:
			res[ch] = res[r] ^ w
			yield dfs(ch, r)
	yield None
t = 1
mod = 998244353
for i in range(t):
	(n, m) = RL()
	d = []
	for i in range(n):
		d.append(RLL())
	ifact(n, mod)
	ans = 0
	for j in range(m):
		res = []
		for i in range(n):
			ma = d[i][j] - 1
			res.append(ma)
		res.sort()
		tmp = 1
		for (i, v) in enumerate(res):
			if v <= i:
				tmp = 0
				break
			else:
				tmp = tmp * (v - i) % mod
		ans += 1 - tmp * ifa[n] % mod
		ans %= mod
	print(ans)


import sys
from sys import stdin
ans = 0
(n, m) = map(int, stdin.readline().split())
mod = 998244353
a = [list(map(int, stdin.readline().split())) for i in range(n)]
fac = 1
for i in range(1, n + 1):
	fac *= i
inv = pow(fac, mod - 2, mod)
for j in range(m):
	na = [a[i][j] for i in range(n)]
	na.sort()
	now = 1
	able = 0
	for i in range(n):
		while len(na) > 0 and na[-1] > n - i:
			del na[-1]
			able += 1
		now *= able
		able -= 1
	ans += now * inv
	ans %= mod
print((m - ans) % mod)


MOD = 998244353

def inv_f(x):
	return pow(x, MOD - 2, MOD)
inv = [0] * 50
for i in range(1, 50):
	inv[i] = inv_f(i)
(n, m) = map(int, input().split())
board = []
for i in range(n):
	board.append(list(map(int, input().split())))
out = 0
for j in range(m):
	count = [0] * (n + 1)
	for i in range(n):
		count[board[i][j] - 1] += 1
	p = 0
	tot = n
	for i in range(n):
		tot -= count[~i]
		pp = tot * inv[n - i] % MOD
		p = (p + (1 - p) * pp) % MOD
	out += p
	out %= MOD
print(out)


from sys import stdin, stdout
import heapq
from collections import defaultdict
import math
import bisect
import io, os

def main():
	MOD = 998244353
	(n, m) = list(map(int, stdin.readline().split()))
	dp_div = [0] * (n + 2)
	dp_div[n + 1] = 1
	for x in range(n, 1, -1):
		dp_div[x] = dp_div[x + 1] * x % MOD
	arr = []
	for _ in range(n):
		arr.append(list(map(int, stdin.readline().split())))
	ans = 0
	for i in range(m):
		v = [0] * (n + 2)
		for j in range(n):
			x = arr[j][i]
			v[x] += 1
		pos = 0
		cnt = 0
		for j in range(1, n + 1):
			if pos == -1:
				break
			while v[j] > 0:
				v[j] -= 1
				if j == 1:
					pos = -1
					break
				if cnt == 0:
					pos = n - j + 1
					cnt += 1
					continue
				xx = n - j + 1
				yy = n - cnt
				if xx >= yy:
					pos = -1
					break
				pos = (pos * yy + xx * dp_div[yy + 1] - pos * xx) % MOD
				cnt += 1
		if pos == -1:
			ans = (ans + dp_div[2]) % MOD
		else:
			while cnt != n:
				pos = pos * (n - cnt) % MOD
				cnt += 1
			ans = (ans + pos) % MOD
		y_inv = 1
		for z in range(2, n + 1):
			zz = pow(z, MOD - 2, MOD)
			y_inv = y_inv * zz % MOD
	print(ans * y_inv % MOD)
main()


mod = 998244353

def gcd(a, b):
	if a == 0 or b == 0:
		return 1
	while b:
		(a, b) = (b, a % b)
	return a

def gcdExtended(a, b):
	if a == 0:
		return (b, 0, 1)
	(gcd, x1, y1) = gcdExtended(b % a, a)
	x = y1 - b // a * x1
	y = x1
	return (gcd, x, y)

def simpl(a):
	gecd = gcd(a[0], a[1])
	return (a[0] // gecd, a[1] // gecd)

def mulmod(a, b):
	a = summod(a, 0)
	b = summod(b, 0)
	return a % mod * (b % mod) % mod

def summod(a, b):
	return (a % mod + b % mod + mod) % mod

def addfrac(a, b):
	a = (summod(a[0], 0), summod(a[1], 0))
	b = (summod(b[0], 0), summod(b[1], 0))
	gecd = gcd(a[1], b[1])
	lcm = a[1] // gecd * b[1]
	x = mulmod(a[0], b[1] // gecd)
	y = mulmod(b[0], a[1] // gecd)
	return simpl((summod(x, y), lcm))

def mulfrac(a, b):
	return simpl((mulmod(a[0], b[0]), mulmod(a[1], b[1])))
(n, m) = map(int, input().split())
dists = [[0 for j in range(n)] for i in range(m)]
out = (0, 1)
for i in range(n):
	dist = list(map(int, input().split()))
	for j in range(m):
		dists[j][i] = dist[j]
for i in range(m):
	dists[i].sort()
out = 0
nfac = 1
for i in range(2, n + 1):
	nfac = mulmod(nfac, i)
for i in range(m):
	cur = 1
	for j in range(n):
		cur = mulmod(cur, max(0, dists[i][j] - (j + 1)))
	out = summod(out, cur)
out = mulfrac((-1, 1), (out, nfac))
out = addfrac((m, 1), out)
out = mulmod(out[0], (gcdExtended(out[1], mod)[1] + mod) % mod)
print(out)


mod = 998244353
eps = 10 ** (-9)

def main():
	import sys
	from bisect import bisect_left
	input = sys.stdin.readline
	(N, M) = map(int, input().split())
	A = [[0] * N for _ in range(M)]
	for i in range(N):
		line = list(map(int, input().split()))
		for j in range(M):
			A[j][i] = line[j]
	F = 1
	for i in range(1, N + 1):
		F = F * i % mod
	invF = pow(F, mod - 2, mod)
	ans = 0
	for B in A:
		B.sort()
		tmp = 1
		for i in range(N, 0, -1):
			j = bisect_left(B, i + 1)
			k = max(0, N - j - (N - i))
			tmp = tmp * k % mod
		ans = (ans + (F - tmp)) % mod
	print(ans * invF % mod)
main()


import sys
input = sys.stdin.readline
(fact, inv_fact) = ([0] * 21, [0] * 21)
fact[0] = 1
mod = 998244353

def make_nCr_mod(max_n=20, mod=998244353):
	global fact
	global inv_fact
	for i in range(max_n):
		fact[i + 1] = fact[i] * (i + 1) % mod
	inv_fact[-1] = pow(fact[-1], mod - 2, mod)
	for i in reversed(range(max_n)):
		inv_fact[i] = inv_fact[i + 1] * (i + 1) % mod
make_nCr_mod()

def nCr_mod(n, r):
	global fact
	global inv_fact
	res = 1
	while n or r:
		(a, b) = (n % mod, r % mod)
		if a < b:
			return 0
		res = res * fact[a] % mod * inv_fact[b] % mod * inv_fact[a - b] % mod
		n //= mod
		r //= mod
	return res
(n, m) = map(int, input().split())
a = [list(map(int, input().split())) for _ in range(n)]
freq = [[0 for j in range(n + 2)] for i in range(m)]
for j in range(m):
	for i in range(n):
		freq[j][a[i][j]] += 1
ans = 0
for i in range(m):
	counted = 0
	bad = 1
	for j in range(1, n + 2):
		if freq[i][j] == 0:
			continue
		if freq[i][j] >= j - counted:
			bad = 0
			break
		bad *= nCr_mod(j - 1 - counted, freq[i][j]) * fact[freq[i][j]]
		counted += freq[i][j]
	ans += (fact[n] - bad) % mod
print(ans * pow(fact[n], mod - 2, mod) % mod)


import sys
import math
from math import factorial, inf, gcd, sqrt
from heapq import *
from functools import *
from itertools import *
from collections import *
from typing import *
from bisect import *
import random
sys.setrecursionlimit(10 ** 5)

def rarray():
	return [int(i) for i in input().split()]
(n, m) = rarray()
mod = 998244353
a = [rarray() for i in range(n)]
a = list(zip(*a))
ans = 0
fac = 1
for i in range(1, n + 1):
	fac = fac * i % mod
for i in range(m):
	c = [0] * (n + 1)
	for j in a[i]:
		c[j - 1] += 1
	s = 0
	res = 1
	for j in reversed(range(1, n + 1)):
		s += c[j]
		res = res * s % mod
		s -= 1
	ans = (ans + fac - res + mod) % mod
print(ans * pow(fac, mod - 2, mod) % mod)


import sys
input = sys.stdin.readline
MOD = 998244353

def main():
	(n, m) = map(int, input().split())
	dlst = [list(map(int, input().split())) for _ in range(n)]
	ans = 0
	fact = 1
	for i in range(2, n + 1):
		fact *= i
		fact %= MOD
	for i in range(m):
		lst = [dlst[j][i] for j in range(n)]
		lst.sort(reverse=True)
		lst += [0]
		ret = 1
		s = 0
		pos = 0
		for j in range(n, 0, -1):
			while lst[pos] > j:
				s += 1
				pos += 1
			ret *= s
			s -= 1
		ans += fact - ret
		ans %= MOD
	print(ans * pow(fact, MOD - 2, MOD) % MOD)
for _ in range(1):
	main()


p = 998244353

def f(x, y):
	r = 1
	x = x % p
	if x == 0:
		return 0
	while y > 0:
		if y % 2 == 1:
			r = r * x % p
		y = y >> 1
		x = x * x % p
	return r
(n, m) = map(int, input().split())
l = []
for _ in range(n):
	l.append(list(map(int, input().split())))
a = 1
for i in range(1, n + 1):
	a = a * i % p
q = 0
for i in range(m):
	v = [0] * n
	for j in range(n):
		v[j] = l[j][i]
	v.sort()
	x = 1
	for j in range(n):
		x = x * (min(v[j] - 1, n) - j) % p
	q = (q + a - x) % p
print(q * f(a, p - 2) % p)


from sys import stdin
mod = 998244353
(n, m) = map(int, input().split())
arr = []
for i in range(m):
	arr.append([])
for i in range(n):
	temp = list(map(int, stdin.readline().split()))
	for j in range(m):
		arr[j].append(temp[j])
for i in range(m):
	arr[i].sort()
ans = 1
for i in range(1, n + 1):
	ans *= i
den = ans
ans *= m
ans %= mod
for i in range(m):
	temp = 1
	for j in range(n):
		temp *= arr[i][j] - (j + 1)
	ans -= temp
	ans %= mod
ans *= pow(den, mod - 2, mod)
ans %= mod
print(ans)


def divisors(M):
	d = []
	i = 1
	while M >= i ** 2:
		if M % i == 0:
			d.append(i)
			if i ** 2 != M:
				d.append(M // i)
		i = i + 1
	return d

def popcount(x):
	x = x - (x >> 1 & 1431655765)
	x = (x & 858993459) + (x >> 2 & 858993459)
	x = x + (x >> 4) & 252645135
	x = x + (x >> 8)
	x = x + (x >> 16)
	return x & 127

def eratosthenes(n):
	res = [0 for i in range(n + 1)]
	prime = set([])
	for i in range(2, n + 1):
		if not res[i]:
			prime.add(i)
			for j in range(1, n // i + 1):
				res[i * j] = 1
	return prime

def factorization(n):
	res = []
	for p in prime:
		if n % p == 0:
			while n % p == 0:
				n //= p
			res.append(p)
	if n != 1:
		res.append(n)
	return res

def euler_phi(n):
	res = n
	for x in range(2, n + 1):
		if x ** 2 > n:
			break
		if n % x == 0:
			res = res // x * (x - 1)
			while n % x == 0:
				n //= x
	if n != 1:
		res = res // n * (n - 1)
	return res

def ind(b, n):
	res = 0
	while n % b == 0:
		res += 1
		n //= b
	return res

def isPrimeMR(n):
	d = n - 1
	d = d // (d & -d)
	L = [2, 3, 5, 7, 11, 13, 17]
	for a in L:
		t = d
		y = pow(a, t, n)
		if y == 1:
			continue
		while y != n - 1:
			y = y * y % n
			if y == 1 or t == n - 1:
				return 0
			t <<= 1
	return 1

def findFactorRho(n):
	from math import gcd
	m = 1 << n.bit_length() // 8
	for c in range(1, 99):
		f = lambda x: (x * x + c) % n
		(y, r, q, g) = (2, 1, 1, 1)
		while g == 1:
			x = y
			for i in range(r):
				y = f(y)
			k = 0
			while k < r and g == 1:
				ys = y
				for i in range(min(m, r - k)):
					y = f(y)
					q = q * abs(x - y) % n
				g = gcd(q, n)
				k += m
			r <<= 1
		if g == n:
			g = 1
			while g == 1:
				ys = f(ys)
				g = gcd(abs(x - ys), n)
		if g < n:
			if isPrimeMR(g):
				return g
			elif isPrimeMR(n // g):
				return n // g
			return findFactorRho(g)

def primeFactor(n):
	i = 2
	ret = {}
	rhoFlg = 0
	while i * i <= n:
		k = 0
		while n % i == 0:
			n //= i
			k += 1
		if k:
			ret[i] = k
		i += 1 + i % 2
		if i == 101 and n >= 2 ** 20:
			while n > 1:
				if isPrimeMR(n):
					(ret[n], n) = (1, 1)
				else:
					rhoFlg = 1
					j = findFactorRho(n)
					k = 0
					while n % j == 0:
						n //= j
						k += 1
					ret[j] = k
	if n > 1:
		ret[n] = 1
	if rhoFlg:
		ret = {x: ret[x] for x in sorted(ret)}
	return ret

def divisors(n):
	res = [1]
	prime = primeFactor(n)
	for p in prime:
		newres = []
		for d in res:
			for j in range(prime[p] + 1):
				newres.append(d * p ** j)
		res = newres
	res.sort()
	return res

def xorfactorial(num):
	if num == 0:
		return 0
	elif num == 1:
		return 1
	elif num == 2:
		return 3
	elif num == 3:
		return 0
	else:
		x = baseorder(num)
		return 2 ** x * ((num - 2 ** x + 1) % 2) + function(num - 2 ** x)

def xorconv(n, X, Y):
	if n == 0:
		res = [X[0] * Y[0] % mod]
		return res
	x = [X[i] + X[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	y = [Y[i] + Y[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	z = [X[i] - X[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	w = [Y[i] - Y[i + 2 ** (n - 1)] for i in range(2 ** (n - 1))]
	res1 = xorconv(n - 1, x, y)
	res2 = xorconv(n - 1, z, w)
	former = [(res1[i] + res2[i]) * inv for i in range(2 ** (n - 1))]
	latter = [(res1[i] - res2[i]) * inv for i in range(2 ** (n - 1))]
	former = list(map(lambda x: x % mod, former))
	latter = list(map(lambda x: x % mod, latter))
	return former + latter

def merge_sort(A, B):
	(pos_A, pos_B) = (0, 0)
	(n, m) = (len(A), len(B))
	res = []
	while pos_A < n and pos_B < m:
		(a, b) = (A[pos_A], B[pos_B])
		if a < b:
			res.append(a)
			pos_A += 1
		else:
			res.append(b)
			pos_B += 1
	res += A[pos_A:]
	res += B[pos_B:]
	return res

class UnionFindVerSize:

	def __init__(self, N):
		self._parent = [n for n in range(0, N)]
		self._size = [1] * N
		self.group = N

	def find_root(self, x):
		if self._parent[x] == x:
			return x
		self._parent[x] = self.find_root(self._parent[x])
		stack = [x]
		while self._parent[stack[-1]] != stack[-1]:
			stack.append(self._parent[stack[-1]])
		for v in stack:
			self._parent[v] = stack[-1]
		return self._parent[x]

	def unite(self, x, y):
		gx = self.find_root(x)
		gy = self.find_root(y)
		if gx == gy:
			return
		self.group -= 1
		if self._size[gx] < self._size[gy]:
			self._parent[gx] = gy
			self._size[gy] += self._size[gx]
		else:
			self._parent[gy] = gx
			self._size[gx] += self._size[gy]

	def get_size(self, x):
		return self._size[self.find_root(x)]

	def is_same_group(self, x, y):
		return self.find_root(x) == self.find_root(y)

class WeightedUnionFind:

	def __init__(self, N):
		self.parent = [i for i in range(N)]
		self.size = [1 for i in range(N)]
		self.val = [0 for i in range(N)]
		self.flag = True
		self.edge = [[] for i in range(N)]

	def dfs(self, v, pv):
		stack = [(v, pv)]
		new_parent = self.parent[pv]
		while stack:
			(v, pv) = stack.pop()
			self.parent[v] = new_parent
			for (nv, w) in self.edge[v]:
				if nv != pv:
					self.val[nv] = self.val[v] + w
					stack.append((nv, v))

	def unite(self, x, y, w):
		if not self.flag:
			return
		if self.parent[x] == self.parent[y]:
			self.flag = self.val[x] - self.val[y] == w
			return
		if self.size[self.parent[x]] > self.size[self.parent[y]]:
			self.edge[x].append((y, -w))
			self.edge[y].append((x, w))
			self.size[x] += self.size[y]
			self.val[y] = self.val[x] - w
			self.dfs(y, x)
		else:
			self.edge[x].append((y, -w))
			self.edge[y].append((x, w))
			self.size[y] += self.size[x]
			self.val[x] = self.val[y] + w
			self.dfs(x, y)

class Dijkstra:

	class Edge:

		def __init__(self, _to, _cost):
			self.to = _to
			self.cost = _cost

	def __init__(self, V):
		self.G = [[] for i in range(V)]
		self._E = 0
		self._V = V

	@property
	def E(self):
		return self._E

	@property
	def V(self):
		return self._V

	def add_edge(self, _from, _to, _cost):
		self.G[_from].append(self.Edge(_to, _cost))
		self._E += 1

	def shortest_path(self, s):
		import heapq
		que = []
		d = [10 ** 15] * self.V
		d[s] = 0
		heapq.heappush(que, (0, s))
		while len(que) != 0:
			(cost, v) = heapq.heappop(que)
			if d[v] < cost:
				continue
			for i in range(len(self.G[v])):
				e = self.G[v][i]
				if d[e.to] > d[v] + e.cost:
					d[e.to] = d[v] + e.cost
					heapq.heappush(que, (d[e.to], e.to))
		return d

def Z_algorithm(s):
	N = len(s)
	Z_alg = [0] * N
	Z_alg[0] = N
	i = 1
	j = 0
	while i < N:
		while i + j < N and s[j] == s[i + j]:
			j += 1
		Z_alg[i] = j
		if j == 0:
			i += 1
			continue
		k = 1
		while i + k < N and k + Z_alg[k] < j:
			Z_alg[i + k] = Z_alg[k]
			k += 1
		i += k
		j -= k
	return Z_alg

class BIT:

	def __init__(self, n, mod=0):
		self.BIT = [0] * (n + 1)
		self.num = n
		self.mod = mod

	def query(self, idx):
		res_sum = 0
		mod = self.mod
		while idx > 0:
			res_sum += self.BIT[idx]
			if mod:
				res_sum %= mod
			idx -= idx & -idx
		return res_sum

	def update(self, idx, x):
		mod = self.mod
		while idx <= self.num:
			self.BIT[idx] += x
			if mod:
				self.BIT[idx] %= mod
			idx += idx & -idx
		return

class dancinglink:

	def __init__(self, n, debug=False):
		self.n = n
		self.debug = debug
		self._left = [i - 1 for i in range(n)]
		self._right = [i + 1 for i in range(n)]
		self.exist = [True for i in range(n)]

	def pop(self, k):
		if self.debug:
			assert self.exist[k]
		L = self._left[k]
		R = self._right[k]
		if L != -1:
			if R != self.n:
				(self._right[L], self._left[R]) = (R, L)
			else:
				self._right[L] = self.n
		elif R != self.n:
			self._left[R] = -1
		self.exist[k] = False

	def left(self, idx, k=1):
		if self.debug:
			assert self.exist[idx]
		res = idx
		while k:
			res = self._left[res]
			if res == -1:
				break
			k -= 1
		return res

	def right(self, idx, k=1):
		if self.debug:
			assert self.exist[idx]
		res = idx
		while k:
			res = self._right[res]
			if res == self.n:
				break
			k -= 1
		return res

class SparseTable:

	def __init__(self, A, merge_func, ide_ele):
		N = len(A)
		n = N.bit_length()
		self.table = [[ide_ele for i in range(n)] for i in range(N)]
		self.merge_func = merge_func
		for i in range(N):
			self.table[i][0] = A[i]
		for j in range(1, n):
			for i in range(0, N - 2 ** j + 1):
				f = self.table[i][j - 1]
				s = self.table[i + 2 ** (j - 1)][j - 1]
				self.table[i][j] = self.merge_func(f, s)

	def query(self, s, t):
		b = t - s + 1
		m = b.bit_length() - 1
		return self.merge_func(self.table[s][m], self.table[t - 2 ** m + 1][m])

class BinaryTrie:

	class node:

		def __init__(self, val):
			self.left = None
			self.right = None
			self.max = val

	def __init__(self):
		self.root = self.node(-10 ** 15)

	def append(self, key, val):
		pos = self.root
		for i in range(29, -1, -1):
			pos.max = max(pos.max, val)
			if key >> i & 1:
				if pos.right is None:
					pos.right = self.node(val)
					pos = pos.right
				else:
					pos = pos.right
			elif pos.left is None:
				pos.left = self.node(val)
				pos = pos.left
			else:
				pos = pos.left
		pos.max = max(pos.max, val)

	def search(self, M, xor):
		res = -10 ** 15
		pos = self.root
		for i in range(29, -1, -1):
			if pos is None:
				break
			if M >> i & 1:
				if xor >> i & 1:
					if pos.right:
						res = max(res, pos.right.max)
					pos = pos.left
				else:
					if pos.left:
						res = max(res, pos.left.max)
					pos = pos.right
			elif xor >> i & 1:
				pos = pos.right
			else:
				pos = pos.left
		if pos:
			res = max(res, pos.max)
		return res

def solveequation(edge, ans, n, m):
	x = [0] * m
	used = [False] * n
	for v in range(n):
		if used[v]:
			continue
		y = dfs(v)
		if y != 0:
			return False
	return x

	def dfs(v):
		used[v] = True
		r = ans[v]
		for (to, dire, id) in edge[v]:
			if used[to]:
				continue
			y = dfs(to)
			if dire == -1:
				x[id] = y
			else:
				x[id] = -y
			r += y
		return r

class SegmentTree:

	def __init__(self, init_val, segfunc, ide_ele):
		n = len(init_val)
		self.segfunc = segfunc
		self.ide_ele = ide_ele
		self.num = 1 << (n - 1).bit_length()
		self.tree = [ide_ele] * 2 * self.num
		self.size = n
		for i in range(n):
			self.tree[self.num + i] = init_val[i]
		for i in range(self.num - 1, 0, -1):
			self.tree[i] = self.segfunc(self.tree[2 * i], self.tree[2 * i + 1])

	def update(self, k, x):
		k += self.num
		self.tree[k] = x
		while k > 1:
			self.tree[k >> 1] = self.segfunc(self.tree[k], self.tree[k ^ 1])
			k >>= 1

	def query(self, l, r):
		if r == self.size:
			r = self.num
		res = self.ide_ele
		l += self.num
		r += self.num
		while l < r:
			if l & 1:
				res = self.segfunc(res, self.tree[l])
				l += 1
			if r & 1:
				res = self.segfunc(res, self.tree[r - 1])
			l >>= 1
			r >>= 1
		return res

	def bisect_l(self, l, r, x):
		l += self.num
		r += self.num
		Lmin = -1
		Rmin = -1
		while l < r:
			if l & 1:
				if self.tree[l] <= x and Lmin == -1:
					Lmin = l
				l += 1
			if r & 1:
				if self.tree[r - 1] <= x:
					Rmin = r - 1
			l >>= 1
			r >>= 1
		if Lmin != -1:
			pos = Lmin
			while pos < self.num:
				if self.tree[2 * pos] <= x:
					pos = 2 * pos
				else:
					pos = 2 * pos + 1
			return pos - self.num
		elif Rmin != -1:
			pos = Rmin
			while pos < self.num:
				if self.tree[2 * pos] <= x:
					pos = 2 * pos
				else:
					pos = 2 * pos + 1
			return pos - self.num
		else:
			return -1

def cmb(n, r, mod):
	if r < 0 or r > n:
		return 0
	r = min(r, n - r)
	return g1[n] * g2[r] * g2[n - r] % mod
mod = 998244353
N = 100
g1 = [1] * (N + 1)
g2 = [1] * (N + 1)
inverse = [1] * (N + 1)
for i in range(2, N + 1):
	g1[i] = g1[i - 1] * i % mod
	inverse[i] = -inverse[mod % i] * (mod // i) % mod
	g2[i] = g2[i - 1] * inverse[i] % mod
inverse[0] = 0
import sys, random, bisect
from collections import deque, defaultdict
from heapq import heapify, heappop, heappush
from itertools import permutations
from math import gcd, log
input = lambda : sys.stdin.buffer.readline()
mi = lambda : map(int, input().split())
li = lambda : list(mi())
(n, m) = mi()
d = [li() for i in range(n)]
res = m
for i in range(m):
	cnt = [0 for j in range(n)]
	flag = True
	for j in range(n):
		if d[j][i] == 1:
			flag = False
			break
		else:
			cnt[d[j][i] - 2] += 1
	if not flag:
		continue
	tmp = 1
	rest = 0
	for j in range(n):
		rest += 1
		for k in range(cnt[j]):
			tmp *= rest
			tmp %= mod
			rest -= 1
	res -= tmp * g2[n]
	res %= mod
print(res)


import sys, io, os
BUFSIZE = 8192

class FastIO(io.IOBase):
	newlines = 0

	def __init__(self, file):
		self._fd = file.fileno()
		self.buffer = io.BytesIO()
		self.writable = 'x' in file.mode or 'r' not in file.mode
		self.write = self.buffer.write if self.writable else None

	def read(self):
		while True:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			if not b:
				break
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines = 0
		return self.buffer.read()

	def readline(self):
		while self.newlines == 0:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			self.newlines = b.count(b'\n') + (not b)
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines -= 1
		return self.buffer.readline()

	def flush(self):
		if self.writable:
			os.write(self._fd, self.buffer.getvalue())
			(self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(io.IOBase):

	def __init__(self, file):
		self.buffer = FastIO(file)
		self.flush = self.buffer.flush
		self.writable = self.buffer.writable
		self.write = lambda s: self.buffer.write(s.encode('ascii'))
		self.read = lambda : self.buffer.read().decode('ascii')
		self.readline = lambda : self.buffer.readline().decode('ascii')

def print(*args, **kwargs):
	(sep, file) = (kwargs.pop('sep', ' '), kwargs.pop('file', sys.stdout))
	at_start = True
	for x in args:
		if not at_start:
			file.write(sep)
		file.write(str(x))
		at_start = False
	file.write(kwargs.pop('end', '\n'))
	if kwargs.pop('flush', False):
		file.flush()
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')
mod = 998244353
(N, M) = map(int, input().split())
D = [list(map(int, input().split())) for _ in range(N)]
ans = 0
for t in zip(*D):
	t = sorted(t)
	p = 1
	for (i, d) in enumerate(t, 1):
		d -= i
		if d <= 0:
			break
		p = p * d % mod
	else:
		ans += p
denom = 1
for i in range(1, N + 1):
	denom = denom * i % mod
ans = (M - ans * pow(denom, mod - 2, mod)) % mod
print(ans)


from math import factorial
from fractions import Fraction
(n, m) = map(int, input().split())
dji = [[] for _ in range(m)]
for i in range(n):
	for (j, x) in enumerate(map(int, input().split())):
		dji[j].append(x)
out = 0
for dj in dji:
	dj.sort()
	p = 1
	for (i, d) in enumerate(dj):
		if i + 1 >= d:
			p = 0
			break
		p *= d - i - 1
	out += p
f = m - Fraction(out, factorial(n))
print(f.numerator * pow(f.denominator, -1, 998244353) % 998244353)


import sys
import getpass
import math, random
import functools, itertools, collections, heapq, bisect
from collections import Counter, defaultdict, deque
input = sys.stdin.readline
M9 = 998244353
(yes, no) = ('YES', 'NO')
MAXINT = sys.maxsize
OFFLINE_TEST = getpass.getuser() == 'hkmac'

def log(*args):
	if OFFLINE_TEST:
		print('\x1b[36m', *args, '\x1b[0m', file=sys.stderr)

def solve(*args):
	if OFFLINE_TEST:
		log('----- solving ------')
		log(*args)
		log('----- ------- ------')
	return solve_(*args)

def read_matrix(rows):
	return [list(map(int, input().split())) for _ in range(rows)]

def read_strings(rows):
	return [input().strip() for _ in range(rows)]

def minus_one(arr):
	return [x - 1 for x in arr]

def minus_one_matrix(mrr):
	return [[x - 1 for x in row] for row in mrr]

def factorial_mod_p(n, p):
	val = 1
	for i in range(1, n + 1):
		val = val * i % p
	return val

def solve_(mrr, n, m):
	val = 0
	for col in zip(*mrr):
		col = sorted(col)
		numer = 1
		for (i, x) in enumerate(col):
			dem = n - i
			num = dem - max(0, n - x + 1)
			numer *= num
		val += math.factorial(n) - numer
	return val * pow(factorial_mod_p(n, M9), M9 - 2, M9) % M9
for case_num in [0]:
	(n, m) = list(map(int, input().split()))
	mrr = read_matrix(n)
	res = solve(mrr, n, m)
	print(res)


import sys
input = sys.stdin.readline
get_int = lambda : int(input().rstrip())
get_arr = lambda : [int(w) for w in input().split()]
get_str = lambda : input().rstrip()
(n, m) = get_arr()
d = [get_arr() for _ in range(n)]
mod = 998244353
inv = 1
for t in range(1, n + 1):
	inv *= pow(t, -1, mod)
	inv %= mod
ans = 0
for j in range(m):
	arr = sorted((d[i][j] - 1 for i in range(n)))
	p = 1
	for k in range(n):
		val = arr[k] - k
		if val <= 0:
			p = 0
			break
		p *= val
	ans += (1 - p * inv) % mod
	ans %= mod
print(ans)


import sys
input = sys.stdin.buffer.readline
mod = 998244353
ans = 0
(n, m) = map(int, input().split())
fact = 1
for i in range(2, n + 1):
	fact = fact * i % mod
inv_fact = pow(fact, mod - 2, mod)
dist = [list(map(int, input().split())) for i in range(n)]
for point in range(m):
	time = [n + 1 - dist[city][point] + 1 for city in range(n)]
	add = [0] * (n + 2)
	for i in time:
		add[i] += 1
	dp = [[0] * (n + 1) for i in range(n + 2)]
	dp[0][0] = 1
	for t in range(n + 1):
		for to_place in range(n + 1):
			if dp[t][to_place]:
				dp[t + 1][to_place + add[t + 1]] = (dp[t + 1][to_place + add[t + 1]] + dp[t][to_place]) % mod
				if to_place:
					dp[t + 1][to_place - 1 + add[t + 1]] = (dp[t + 1][to_place - 1 + add[t + 1]] + to_place * dp[t][to_place]) % mod
	ans = (ans + (fact - dp[n + 1][0]) * inv_fact) % mod
print(ans)


import sys, io, os
try:
	Z = io.BytesIO(os.read(0, os.fstat(0).st_size)).readline
except:
	Z = lambda : sys.stdin.readline().encode()
Y = lambda : map(int, Z().split())
M = 998244353
try:
	import __pypy__
	int_add = __pypy__.intop.int_add
	int_sub = __pypy__.intop.int_sub
	int_mul = __pypy__.intop.int_mul

	def make_mod_mul(mod=M):
		fmod_inv = 1.0 / mod

		def mod_mul(a, b, c=0):
			res = int_sub(int_add(int_mul(a, b), c), int_mul(mod, int(fmod_inv * a * b + fmod_inv * c)))
			if res >= mod:
				return res - mod
			elif res < 0:
				return res + mod
			else:
				return res
		return mod_mul
	mod_mul = make_mod_mul()
except:

	def mod_mul(a, b):
		return a * b % M

def mod_add(a, b):
	v = a + b
	if v >= M:
		v -= M
	if v < 0:
		v += M
	return v

def mod_sum(a):
	v = 0
	for i in a:
		v = mod_add(v, i)
	return v
(n, m) = Y()
d = [[*Y()] for i in range(n)]
s = 0
inv = 1
for i in range(n):
	inv = mod_mul(inv, i + 1)
inv = pow(inv, M - 2, M)

def f(l):
	v = 1
	n = len(l)
	p = n - 1
	c = 0
	for i in range(n):
		while p >= 0 and l[p] > n - i:
			c += 1
			p -= 1
		if not c:
			return 0
		v = mod_mul(v, c)
		c -= 1
	return v
for i in range(m):
	l = sorted((d[j][i] for j in range(n)))
	p = mod_mul(f(l), inv)
	s = mod_add(s, 1 - p)
print(s)


def power(a, n, mod):
	answer = 1
	while n > 0:
		if n % 2 == 0:
			n //= 2
			a = a * a % mod
		else:
			answer *= a
			answer %= mod
			n -= 1
	return answer

def putin():
	return map(int, input().split())

def main():
	(n, m) = putin()
	mod = 998244353
	fac = 1
	for i in range(1, n + 1):
		fac *= i
		fac %= mod
	inv = power(fac, mod - 2, mod)
	assert inv * fac % mod == 1
	A = []
	for i in range(n):
		A.append(list(putin()))
	bad_sum = 0
	for i in range(m):
		B = []
		for j in range(n):
			B.append(A[j][i])
		B.sort(reverse=True)
		cnt = 0
		bad = 1
		for j in range(n, 0, -1):
			while cnt < len(B) and B[cnt] > j:
				cnt += 1
			bad *= cnt + j - n
		bad_sum += bad
	print((fac * m - bad_sum) * inv % mod)
main()


import os, sys
from io import BytesIO, IOBase

def main():
	mod = 998244353
	(n, m) = map(int, input().split())
	arr = [list(map(int, input().split())) for _ in range(n)]
	(ans, fac) = (0, 1)
	for i in range(1, n + 1):
		fac = fac * i % mod
	for tt in range(m):
		x = sorted([arr[i][tt] for i in range(n)])
		val = 1
		for i in range(n):
			val = val * (x[i] - i - 1) % mod
		val = (fac - val) % mod
		ans = (ans + val * pow(fac, mod - 2, mod)) % mod
	print(ans)
BUFSIZE = 8192

class FastIO(IOBase):
	newlines = 0

	def __init__(self, file):
		self._fd = file.fileno()
		self.buffer = BytesIO()
		self.writable = 'x' in file.mode or 'r' not in file.mode
		self.write = self.buffer.write if self.writable else None

	def read(self):
		while True:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			if not b:
				break
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines = 0
		return self.buffer.read()

	def readline(self):
		while self.newlines == 0:
			b = os.read(self._fd, max(os.fstat(self._fd).st_size, BUFSIZE))
			self.newlines = b.count(b'\n') + (not b)
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines -= 1
		return self.buffer.readline()

	def flush(self):
		if self.writable:
			os.write(self._fd, self.buffer.getvalue())
			(self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(IOBase):

	def __init__(self, file):
		self.buffer = FastIO(file)
		self.flush = self.buffer.flush
		self.writable = self.buffer.writable
		self.write = lambda s: self.buffer.write(s.encode('ascii'))
		self.read = lambda : self.buffer.read().decode('ascii')
		self.readline = lambda : self.buffer.readline().decode('ascii')
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')
main()


def get_one(ns, m):
	n = len(ns)
	ans = 1
	for i in range(n):
		t = ns[i] - i
		if t <= 0:
			return 0
		ans = ans * t % m
	return ans

def get_one_pre(ns: list):
	ans = [0] * len(ns)
	ns.sort()
	l = 0
	for i in range(len(ns)):
		for j in range(l, len(ns)):
			if i + 1 < ns[j]:
				l = j
				ans[i] = len(ns) - l
				break
		else:
			l = j
			ans[i] = 0
	return list(reversed(ans))

def pow(a, b, m):
	tb = [a]
	for i in range(35):
		tb.append(tb[-1] * tb[-1] % m)
	ans = 1
	for i in range(35):
		if b >> i & 1 == 1:
			ans = ans * tb[i] % m
	return ans

def inv(a, m):
	return pow(a, m - 2, m)

def functoria(x, m):
	ans = 1
	for i in range(1, x + 1):
		ans = ans * i % m
	return ans

def gns():
	return list(map(int, input().split()))

def gn():
	return int(input())
(n, m) = gns()
md = 998244353
ms = [[] for i in range(m)]
for i in range(n):
	ns = gns()
	for i in range(m):
		ms[i].append(ns[i])
sm = 0
for i in range(m):
	x = get_one_pre(ms[i])
	y = get_one(x, md)
	sm = (sm + y) % md
f = functoria(n, md)
fv = inv(f, md)
ans = (f * m - sm + md) * fv % md
print(ans)


import time
import os, sys
from io import BytesIO, IOBase
from types import GeneratorType
from bisect import bisect_left, bisect_right
from collections import defaultdict as dd, deque as dq, Counter as dc
import math, string, heapq as h
BUFSIZE = 8192

class FastIO(IOBase):
	newlines = 0

	def __init__(self, file):
		import os
		self.os = os
		self._fd = file.fileno()
		self.buffer = BytesIO()
		self.writable = 'x' in file.mode or 'r' not in file.mode
		self.write = self.buffer.write if self.writable else None

	def read(self):
		while True:
			b = self.os.read(self._fd, max(self.os.fstat(self._fd).st_size, BUFSIZE))
			if not b:
				break
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines = 0
		return self.buffer.read()

	def readline(self):
		while self.newlines == 0:
			b = self.os.read(self._fd, max(self.os.fstat(self._fd).st_size, BUFSIZE))
			self.newlines = b.count(b'\n') + (not b)
			ptr = self.buffer.tell()
			(self.buffer.seek(0, 2), self.buffer.write(b), self.buffer.seek(ptr))
		self.newlines -= 1
		return self.buffer.readline()

	def flush(self):
		if self.writable:
			self.os.write(self._fd, self.buffer.getvalue())
			(self.buffer.truncate(0), self.buffer.seek(0))

class IOWrapper(IOBase):

	def __init__(self, file):
		self.buffer = FastIO(file)
		self.flush = self.buffer.flush
		self.writable = self.buffer.writable
		self.write = lambda s: self.buffer.write(s.encode('ascii'))
		self.read = lambda : self.buffer.read().decode('ascii')
		self.readline = lambda : self.buffer.readline().decode('ascii')
(sys.stdin, sys.stdout) = (IOWrapper(sys.stdin), IOWrapper(sys.stdout))
input = lambda : sys.stdin.readline().rstrip('\r\n')

def getInt():
	return int(input())

def getStrs():
	return input().split()

def getInts():
	return list(map(int, input().split()))

def getStr():
	return input()

def listStr():
	return list(input())

def getMat(n):
	return [getInts() for _ in range(n)]

def getBin():
	return list(map(int, list(input())))

def isInt(s):
	return '0' <= s[0] <= '9'

def ceil_(a, b):
	return a // b + (a % b > 0)
MOD = 998244353
primes = [2, 3, 5, 7, 11, 13, 17, 19]

def solve():
	(N, M) = getInts()
	dist = [getInts() for _ in range(N)]
	num = 0
	for j in range(M):
		tmp = [dist[i][j] for i in range(N)]
		tmp.sort()
		mult = 1
		for i in range(N):
			mult *= max(tmp[i] - 1 - i, 0)
		num += mult
	denom = 1
	for i in range(2, N + 1):
		denom *= i
	num = denom * M - num
	for i in primes:
		if i > N:
			break
		while num % i == 0 and denom % i == 0:
			num //= i
			denom //= i
	num %= MOD
	denom %= MOD
	return num * pow(denom, MOD - 2, MOD) % MOD
print(solve())


import sys
input = sys.stdin.readline
(n, m) = map(int, input().split())
A = [list(map(int, input().split())) for _ in range(n)]
mod = 998244353
invF = 1
for i in range(1, n + 1):
	invF = invF * i % mod
invF = pow(invF, mod - 2, mod)
ans = 0
for j in range(m):
	cnt = [0] * (n + 1)
	for i in range(n):
		cnt[n + 1 - A[i][j]] += 1
	tot = 0
	cur = 1
	for i in range(n):
		tot += cnt[i]
		cur = cur * tot % mod
		tot = max(0, tot - 1)
	ans = (ans + 1 - cur * invF) % mod
print(ans)


(n, m) = map(int, input().split())
mod = 998244353
dis = [[0] * (m + 1) for i in range(n + 1)]
jc = 1
ans = 0

def ksm(a, b, p):
	a = a % p
	b = b % (p - 1)
	cheng = a
	ret = 1
	while b > 0:
		if b % 2:
			ret = ret * cheng % p
		cheng = cheng * cheng % p
		b //= 2
	return ret

def inv(a, p=mod):
	return ksm(a, p - 2, p)
for i in range(n):
	jc = jc * (i + 1) % mod
for i in range(n):
	t = list(map(int, input().split()))
	for j in range(m):
		dis[i][j] = t[j]
for i in range(m):
	c = [0] * n
	for j in range(n):
		c[j] = dis[j][i] - 1
	c.sort()
	ret = 1
	for (index, j) in enumerate(c):
		ret = ret * max(0, j - index)
	ret = ret * inv(jc) % mod
	ret = ((1 - ret) % mod + mod) % mod
	ans = (ans + ret) % mod
print(ans)