You love array, don't you?

You are given two integer N and M, along with an array B of length N containing pairwise distinct integers. You have to find the number of array A of length N satisfying the follow conditions:

1. 1 ≤ A_{i} ≤ M.
2. A_{i} \neq A_{j} for all 1 ≤ i < j ≤ N.
3. A_{i} \neq B_{i} for all 1 ≤ i ≤ N.

Since the answer can be very large, output it under modulo 10^{9} + 7.

------ Input Format ------ 

- The first line of the input contains two space-separated integers N and M.
- The second line of the input contains N separated integers B_{1}, B_{2}, \dots, B_{N} - the array B.

------ Output Format ------ 

Output on a single line the number of satisfactory array modulo 10^{9} + 7.

------ Constraints ------ 

$1 ≤ N ≤ 2 \cdot 10^{5}$
$1 ≤ M ≤ 3 \cdot 10^{5}$
$1 ≤ B_{i} ≤ 10^{9}$
$B_{i} \neq B_{j}$ for all $1 ≤i < j ≤N$.

------ subtasks ------ 

Subtask $1$ (30 points): $1 ≤ R ≤ 10\,000$
Subtask $2$ (70 points): $1 ≤ R ≤ 10^{9}$

----- Sample Input 1 ------ 
3 4
5 1 2
----- Sample Output 1 ------ 
14
----- explanation 1 ------ 
All satisfactory arrays are:
- $[1, 2, 3]$
- $[2, 3, 1]$
- $[3, 2, 1]$
- $[1, 2, 4]$
- $[2, 4, 1]$
- $[4, 2, 1]$
- $[1, 3, 4]$
- $[1, 4, 3]$
- $[3, 4, 1]$
- $[4, 3, 1]$
- $[2, 3, 4]$
- $[2, 4, 3]$
- $[3, 2, 4]$
- $[4, 2, 3]$

----- Sample Input 2 ------ 
2 2
1 3
----- Sample Output 2 ------ 
1
----- explanation 2 ------ 
The only satisfactory array is:
- $[2, 1]$

Difficulty: VERY_HARD
Tags: Combinatorics

Mathematics

Starter code:

Solutions:
import sys
input = sys.stdin.readline
mod = 10 ** 9 + 7

def ModInverse(n):
	return pow(n, mod - 2, mod)

def NChooseK(n, k):
	return Factorial[n] * ModInverse(Factorial[n - k] * Factorial[k]) % mod
(size, cap) = [int(x) for x in input().strip().split()]
arr = [int(x) for x in input().strip().split()]
if cap < size:
	print(0)
	exit()
Factorial = [1]
for n in range(1, cap + 1):
	Factorial += [Factorial[-1] * n % mod]
num_elements_in_bounds = sum((1 for x in arr if x <= cap))
ans = NChooseK(cap, size) * Factorial[size] % mod
sign = -1
for min_num_overlaps in range(1, num_elements_in_bounds + 1):
	delta = NChooseK(num_elements_in_bounds, min_num_overlaps) * NChooseK(cap - min_num_overlaps, size - min_num_overlaps) * Factorial[size - min_num_overlaps]
	delta %= mod
	ans = (ans + sign * delta) % mod
	sign *= -1
print(ans)


import sys
from bisect import bisect_left as lb
from bisect import bisect_right as rb
from collections import deque
sys.setrecursionlimit(10 ** 8)
from queue import PriorityQueue as pq
from math import gcd
input_ = lambda : sys.stdin.readline().strip('\r\n')
ii = lambda : int(input_())
il = lambda : list(map(int, input_().split()))
ilf = lambda : list(map(float, input_().split()))
lii = lambda : list(map(int, list(ip())))
ip = lambda : input_()
fi = lambda : float(input_())
ap = lambda ab, bc, cd: ab[bc].append(cd)
li = lambda : list(input_())
pr = lambda x: print(x)
prinT = lambda x: print(x)
f = lambda : sys.stdout.flush()
inv = lambda x: pow(x, mod - 2, mod)
dx = [0, 0, 1, -1]
dy = [1, -1, 0, 0]
mod = 10 ** 9 + 7
mod1 = 998244353
fact = [1]
for i in range(1, 4 * 10 ** 5 + 10):
	fact.append(fact[-1] * i % mod)
inv = [pow(fact[-1], mod - 2, mod)]
for i in range(4 * 10 ** 5 + 9, 0, -1):
	inv.append(inv[-1] * i % mod)
inv.reverse()

def ncr(n, r):
	return fact[n] * inv[r] * inv[n - r] % mod
for _ in range(1):
	(n, m) = il()
	a = il()
	if m < n:
		print(0)
		continue
	t = n
	x = 0
	for i in range(n):
		if a[i] > m:
			t -= 1
			x += 1
	dp = [0 for i in range(t + 1)]
	y = m - t
	dp[t] = ncr(y, x) * fact[x] % mod
	y += 1
	x += 1
	for i in range(t - 1, -1, -1):
		dp[i] = (mod + ncr(y, x) * ncr(t, i) * fact[x] % mod - dp[i + 1]) % mod
		y += 1
		x += 1
	print(dp[0])


import math, sys
input = sys.stdin.readline
S = lambda : input().rstrip()
I = lambda : int(S())
M = lambda : map(int, S().split())
L = lambda : list(M())
mod1 = 1000000007
mod2 = 998244353
fact = [1] * 300005
for i in range(1, 300005):
	fact[i] = i * fact[i - 1] % mod1

def C(n, r):
	return fact[n] * pow(fact[r], mod1 - 2, mod1) * pow(fact[n - r], mod1 - 2, mod1) % mod1 if n >= r else 0
(n, m) = M()
a = L()
p = ans = 0
for i in range(n):
	p += a[i] <= m
for i in range(p + 1):
	ans = (ans + pow(-1, i) * C(p, i) * C(m - i, n - i) * fact[n - i]) % mod1
print(ans)


class Solution:

	def __init__(self):
		self.perm = [1] * (3 * 100005)
		self.inverse = [1] * (3 * 100005)
		self.mod = 1000000007

	def pre(self):
		for i in range(2, 3 * 100005):
			self.perm[i] = self.perm[i - 1] * i % self.mod

	def inversepre(self):
		for i in range(2, 3 * 100005):
			self.inverse[i] = pow(self.perm[i], self.mod - 2, self.mod)

	def comb(self, n, c):
		ans = self.perm[n] * self.inverse[n - c]
		ans = ans * self.inverse[c]
		return ans % self.mod
s = Solution()
s.pre()
s.inversepre()
(n, m) = [int(x) for x in input().split()]
l = [int(x) for x in input().split()]
if n > m:
	print(0)
	exit()
ans = 0
p = 0
for i in l:
	if i <= m:
		p += 1
ans += s.comb(m, n) * s.perm[n]
ans = ans % s.mod
for i in range(1, p + 1):
	if i % 2 == 1:
		ans = ans - s.comb(p, i) * s.comb(m - i, n - i) * s.perm[n - i]
		ans = ans % s.mod
	else:
		ans = ans + s.comb(p, i) * s.comb(m - i, n - i) * s.perm[n - i]
		ans = ans % s.mod
print(ans)


import math, bisect, heapq, random, sys, itertools
ints = lambda : list(map(int, input().split()))
alpha = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
p = 10 ** 9 + 7
fact = [1]
for i in range(3 * 10 ** 5 + 100):
	fact.append(fact[-1] * (i + 1) % p)
inv = []
for i in fact:
	inv.append(pow(i, p - 2, p))

def permutate(a, b):
	return fact[a] * inv[a - b] % p

def comb(a, b):
	return fact[a] * inv[a - b] * inv[b] % p
(n, m) = ints()
a = ints()
while len(a) != n:
	a += ints()
if n > m:
	print(0)
	quit()
x = 0
for i in a:
	if i <= m:
		x += 1
ans = permutate(m, n)
for i in range(1, x + 1):
	a = m - i
	b = n - i
	perm = fact[a] * inv[a - b] % p
	a = x
	b = i
	comb = fact[a] * inv[a - b] * inv[b] % p
	k = perm * comb
	if i % 2:
		ans = (ans - k) % p
	else:
		ans = (ans + k) % p
print(ans % p)


from sys import stdin
input = stdin.readline
mod = 10 ** 9 + 7

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

def sub(a, b):
	return (a - b + mod) % mod

def mul(a, b):
	return a % mod * (b % mod) % mod

def ncr(n, r):
	return mul(fact[n], mul(facti[n - r], facti[r]))

def div(a, b):
	return mul(a, pow(b, mod - 2, mod))
(fact, facti) = ([1], [1])
size = 5 * 10 ** 5 + 1
for i in range(1, size):
	fact.append(mul(fact[-1], i))
	facti.append(pow(fact[-1], mod - 2, mod))

def answer():
	if m < n:
		return 0
	total = mul(ncr(m, n), fact[n])
	count = 0
	for i in range(n):
		if b[i] <= m:
			count += 1
	minus = 0
	for i in range(1, count + 1):
		x = mul(ncr(m - i, n - i), fact[n - i])
		x = mul(x, ncr(count, i))
		if i & 1:
			minus = add(minus, x)
		else:
			minus = sub(minus, x)
	ans = sub(total, minus)
	return ans
for T in range(1):
	(n, m) = map(int, input().split())
	b = list(map(int, input().split()))
	print(answer())


import sys
import bisect
from bisect import bisect_left as lb
from bisect import bisect_right as rb
input_ = lambda : sys.stdin.readline().strip('\r\n')
from math import log
from math import gcd
from math import atan2, acos
from random import randint
sa = lambda : input_()
sb = lambda : int(input_())
sc = lambda : input_().split()
sd = lambda : list(map(int, input_().split()))
sflo = lambda : list(map(float, input_().split()))
se = lambda : float(input_())
sf = lambda : list(input_())
flsh = lambda : sys.stdout.flush()
sys.setrecursionlimit(3 * 10 ** 5)
mod = 10 ** 9 + 7
mod1 = 998244353
gp = []
cost = []
dp = []
mx = []
ans1 = []
ans2 = []
special = []
specnode = []
a = 0
kthpar = []

def dfs2(root, par):
	if par != -1:
		dp[root] = dp[par] + 1
	for i in range(1, 20):
		if kthpar[root][i - 1] != -1:
			kthpar[root][i] = kthpar[kthpar[root][i - 1]][i - 1]
	for child in gp[root]:
		if child == par:
			continue
		kthpar[child][0] = root
		dfs(child, root)
ans = 0
a = []
(n, k) = (0, 0)
b = []
vis = []
tot = 0
time = []
time1 = []
adj = []
mx = -1
eps = 1e-07
fact = [0] * (3 * 10 ** 5 + 5)
inv = [0] * (3 * 10 ** 5 + 5)
fact[0] = 1
fact[1] = 1
inv[0] = 1
inv[1] = 1
for i in range(2, 3 * 10 ** 5 + 4):
	fact[i] = fact[i - 1] * i % mod
	inv[i] = pow(fact[i], mod - 2, mod)

def ncr(n, r):
	if r > n or n < 0:
		return 0
	if n == 0 or r == 0:
		return 1
	ans = fact[n] * inv[r] * inv[n - r]
	ans %= mod
	return ans

def hnbhai(tc):
	(n, m) = sd()
	b = sd()
	tot = n
	for i in b:
		if i > m:
			tot -= 1
	dp = [0] * (tot + 1)
	left = n - tot
	dp[tot] = ncr(m - tot, left) * fact[left] % mod
	for i in range(tot - 1, -1, -1):
		dp[i] = ncr(m - i, n - i) * ncr(tot, i) * fact[n - i] % mod - dp[i + 1]
		dp[i] %= mod
	print(dp[0])
for _ in range(1):
	hnbhai(_ + 1)


import sys

def input():
	return sys.stdin.readline().replace('\n', '').strip()
sys.setrecursionlimit(10 ** 9)
mod = 10 ** 9 + 7

def inv(z):
	return infc[z]
sz = 4 * 10 ** 5 + 3
fc = [0] * sz
fc[0] = 1
fc[1] = 1
infc = [0] * sz
infc[0] = 1
infc[1] = 1
for i in range(2, sz):
	fc[i] = fc[i - 1] * i % mod
for i in range(2, sz):
	infc[i] = pow(fc[i], mod - 2, mod)

def fact(x):
	return fc[x]

def ncr(a, b):
	if b > a:
		return 0
	else:
		return fact(a) * inv(a - b) % mod * inv(b) % mod
for _ in range(1):
	(n, m) = list(map(int, input().split()))
	b = list(map(int, input().split()))
	if m < n:
		print(0)
		continue
	no_r = 0
	for i in b:
		if 1 <= i <= m:
			no_r += 1
	ans = fact(m) * inv(m - n) % mod
	for j in range(1, no_r + 1):
		if j % 2 == 0:
			ans += ncr(no_r, j) * fact(m - j) * inv(m - n)
		else:
			ans -= ncr(no_r, j) * fact(m - j) * inv(m - n)
	print(ans % mod)