Problem:
 Score : 500 points    Problem Statement Count the pairs of length- N sequences consisting of integers between 1 and M (inclusive), A_1, A_2, \cdots, A_{N} and B_1, B_2, \cdots, B_{N} , that satisfy all of the following conditions:   A_i \neq B_i , for every i such that 1\leq i\leq N .  A_i \neq A_j and B_i \neq B_j , for every (i, j) such that 1\leq i < j\leq N .   Since the count can be enormous, print it modulo (10^9+7) .      Constraints  1\leq N \leq M \leq 5\times10^5  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M       Output Print the count modulo (10^9+7) .        Sample Input 1 2 2      Sample Output 1 2  A_1=1,A_2=2,B_1=2,B_2=1 and A_1=2,A_2=1,B_1=1,B_2=2 satisfy the conditions.       Sample Input 2 2 3      Sample Output 2 18       Sample Input 3 141421 356237      Sample Output 3 881613484

Solution (Python):
def examA():
    N = I()
    ans = N + N**2 + N**3
    print(ans)
    return

def examB():
    S = SI()
    T = SI()
    ans = 0
    for s,t in zip(S,T):
        if s!=t:
            ans += 1
    print(ans)
    return

def examC():
    N, M, K = LI()
    A = LI()
    B = LI()
    SA = 0
    b = 0
    SB = 0
    for i in range(M):
        if SB+B[i]>K:
            break
        SB += B[i]
        b += 1
    
    ans = 0
    ans = max(ans,b)
    for i in range(N):
        SA += A[i]
        if SA+SB>K:
            while(b>=0):
                if SA+SB<=K:
                    break
                SB -= B[b-1]
                b -= 1
        if b==-1:
            break
        

        ans = max(ans,b+i+1)
    print(ans)
    return

def examD():
    def primes(n):
        is_prime = [True] * (n + 1)
        is_prime[0] = False
        is_prime[1] = False
        for i in range(2, int(n ** 0.5) + 1):
            if not is_prime[i]:
                continue
            for j in range(i * 2, n + 1, i):
                is_prime[j] = False
        return [i for i in range(n + 1) if is_prime[i]]
    N = I()
    P = primes(N)
    
    D = [1]*(N+1)
    for p in P:
        for i in range(1,1+N//p):
            cnt = 1
            cur = i
            while(cur%p==0):
                cnt += 1
                cur //= p
            D[p*i] *= (1+cnt)

    ans = 0
    for i in range(1,N+1):
        ans += i*D[i]
    print(ans)
    return

def examE():
    class combination():
        
        def __init__(self, n, mod):
            self.n = n
            self.mod = mod
            self.fac = [1] * (n + 1)
            self.inv = [1] * (n + 1)
            for j in range(1, n + 1):
                self.fac[j] = self.fac[j - 1] * j % mod

            self.inv[n] = pow(self.fac[n], mod - 2, mod)
            for j in range(n - 1, -1, -1):
                self.inv[j] = self.inv[j + 1] * (j + 1) % mod

        def comb(self, n, r):
            if r > n or n < 0 or r < 0:
                return 0
            return self.fac[n] * self.inv[n - r] * self.inv[r] % self.mod
    N, M = LI()
    
    dp = [0] * (N + 1)
    dp[2] = 1
    for n in range(3, N + 1):
        dp[n] += (dp[n - 1] + dp[n - 2]) * (n - 1)
        dp[n] %= mod
    
    C = combination(M,mod)
    Acur = C.fac[M]
    Acur %= mod
    
    for i in range(1,N+1):
        cur = C.comb(N,i)*C.fac[M-i]*C.inv[M-N]
        cur %= mod
        if i%2==1:
            Acur = Acur + mod - cur
        else:
            Acur = Acur + mod + cur
        Acur %= mod
        
    if N%2==1:
        Acur += 1
    ans = dp[N] * ((C.fac[M] * (C.inv[M-N])))
    ans %= mod
    print(ans)
    return

def examE2():
    class combination():
        
        def __init__(self, n, mod):
            self.n = n
            self.mod = mod
            self.fac = [1] * (n + 1)
            self.inv = [1] * (n + 1)
            for j in range(1, n + 1):
                self.fac[j] = self.fac[j - 1] * j % mod

            self.inv[n] = pow(self.fac[n], mod - 2, mod)
            for j in range(n - 1, -1, -1):
                self.inv[j] = self.inv[j + 1] * (j + 1) % mod

        def comb(self, n, r):
            if r > n or n < 0 or r < 0:
                return 0
            return self.fac[n] * self.inv[n - r] * self.inv[r] % self.mod

    N, M = LI()
    C = combination(M,mod)
    cnt = 0
    for i in range(N+1):
        cur = C.comb(N,i)*C.fac[M-i]*C.inv[M-N] % mod
        if i%2==1:
            cnt -= cur
            cnt += mod
            cnt %= mod
        else:
            cnt += cur
            cnt %= mod
    ans = cnt * (C.fac[M]*C.inv[M-N])
    ans %= mod
    print(ans)
    return

def examF():
    ans = 0
    print(ans)
    return

from decimal import getcontext,Decimal as dec
import sys,bisect,itertools,heapq,math,random
from copy import deepcopy
from heapq import heappop,heappush,heapify
from collections import Counter,defaultdict,deque
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
def I(): return int(input())
def LI(): return list(map(int,sys.stdin.readline().split()))
def DI(): return dec(input())
def LDI(): return list(map(dec,sys.stdin.readline().split()))
def LSI(): return list(map(str,sys.stdin.readline().split()))
def LS(): return sys.stdin.readline().split()
def SI(): return sys.stdin.readline().strip()
global mod,mod2,inf,alphabet,_ep
mod = 10**9 + 7
mod2 = 998244353
inf = 10**18
_ep = dec("0.000000000001")
alphabet = [chr(ord('a') + i) for i in range(26)]
alphabet_convert = {chr(ord('a') + i): i for i in range(26)}

getcontext().prec = 28

sys.setrecursionlimit(10**7)

if __name__ == '__main__':
    examE2()

"""
142
12 9 1445 0 1
asd dfg hj o o
aidn
"""

Solution (Python):
mod = 10 ** 9 + 7
N, M = map(int, input().split())

class Combination:
    def __init__(self, size):
        self.size = size + 2
        self.fact = [1, 1] + [0] * size
        self.factInv = [1, 1] + [0] * size
        self.inv = [0, 1] + [0] * size

        for i in range(2, self.size):
            self.fact[i] = self.fact[i - 1] * i % mod
            self.inv[i] = -self.inv[mod % i] * (mod // i) % mod
            self.factInv[i] = self.factInv[i - 1] * self.inv[i] % mod

    def npr(self, n, r):
        if n < r or n < 0 or r < 0:
            return 0
        return self.fact[n] * self.factInv[n - r] % mod

    def ncr(self, n, r):
        if n < r or n < 0 or r < 0:
            return 0
        return self.fact[n] * (self.factInv[r] * self.factInv[n - r] % mod) % mod

    def nhr(self, n, r):  
        return self.ncr(n + r - 1, n - 1)

comb = Combination(M + 100)

A = comb.npr(M, N) % mod

B = 0
for i in range(N + 1):
    B += (-1) ** i * comb.ncr(N, i) * comb.npr(M - i, N - i) % mod

print(A * B % mod)

Solution (Python):
import sys
import math
import heapq
sys.setrecursionlimit(10**7)
INTMAX = 9223372036854775807
INTMIN = -9223372036854775808
DVSR = 1000000007
def POW(x, y): return pow(x, y, DVSR)
def INV(x, m=DVSR): return pow(x, m - 2, m)
def DIV(x, y, m=DVSR): return (x * INV(y, m)) % m
def LI(): return map(int, sys.stdin.readline().split())
def LF(): return [float(x) for x in sys.stdin.readline().split()]
def LS(): return sys.stdin.readline().split()
def II(): return int(sys.stdin.readline())


def FLIST(n):
    res = [1]
    for i in range(1, n+1):
        res.append(res[i-1]*i % DVSR)
    return res


N, M = LI()
FACT = FLIST(M)


def NCR(n, r):
    res = 1
    res *= FACT[n-r]
    res *= FACT[r]
    res %= DVSR
    res = INV(res)
    res *= FACT[n]
    res %= DVSR
    return res


def NPR(n, r):
    res = 1
    res *= FACT[n]
    res *= INV(FACT[n-r])
    res %= DVSR
    return res


ALL = NPR(M, N)
ALL %= DVSR
ALL *= ALL
ALL %= DVSR


PAT = 0
for i in reversed(range(1, N+1)):
    v = NCR(N, i)
    v *= NPR(M, i)
    v %= DVSR
    v *= NPR(M-i, N-i)
    v %= DVSR
    v *= NPR(M-i, N-i)
    v %= DVSR
    PAT = v - PAT
    PAT %= DVSR

print((ALL - PAT) % DVSR)

Solution (Python):
N, M = map(int, input().split())
MOD = 10 ** 9 + 7



factorial = [1, 1]
inverse = [1, 1]
invere_base = [0, 1]
for i in range(2, N + M + 1):
    factorial.append((factorial[-1] * i) % MOD)
    invere_base.append((-invere_base[MOD % i] * (MOD // i)) % MOD)
    inverse.append((inverse[-1] * invere_base[-1]) % MOD)


def nCr(n, r):
    if not 0 <= r <= n:
        return 0
    return factorial[n] * inverse[r] * inverse[n - r] % MOD


def nPr(n, r):
    if not 0 <= r <= n:
        return 0
    return nCr(n, r) * factorial[r] % MOD


ans = pow(nPr(M, N), 2, MOD)
for k in range(1, N + 1):
    x = nCr(M, k) * nPr(N, k) * pow(nPr(M - k, N - k), 2, MOD) % MOD
    if k % 2 == 1:
        x = -x
    ans += x
    ans %= MOD

print(ans)

Solution (Python):
from functools import lru_cache

class Combination:
    def __init__(self, size):
        self.size = size + 2
        self.fact = [1, 1] + [0] * size
        self.factInv = [1, 1] + [0] * size
        self.inv = [0, 1] + [0] * size

        for i in range(2, self.size):
            self.fact[i] = self.fact[i - 1] * i % MOD
            self.inv[i] = -self.inv[MOD % i] * (MOD // i) % MOD
            self.factInv[i] = self.factInv[i - 1] * self.inv[i] % MOD

    def npr(self, n, r):
        if n < r or n < 0 or r < 0:
            return 0
        return self.fact[n] * self.factInv[n - r] % MOD

    def ncr(self, n, r):
        if n < r or n < 0 or r < 0:
            return 0
        return self.fact[n] * (self.factInv[r] * self.factInv[n - r] % MOD) % MOD

    def nhr(self, n, r):  
        return self.ncr(n + r - 1, n - 1)

N, M = map(int, input().split())
MOD = 10**9 + 7
comb = Combination(M + 100)

ans = 0
for k in range(N + 1):
    ans += comb.ncr(N, k) * ((-1)**k) * comb.npr(M, k) * (comb.npr(M - k, N - k)**2)
    ans %= MOD
print(ans)

Solution (Python):
import math
N, M = map(int, input().split())

MOD = 10 ** 9 + 7

factorial = [1, 1]
inverse = [1, 1]
invere_base = [0, 1]
for i in range(2, N + M + 1):
    factorial.append((factorial[-1] * i) % MOD)
    invere_base.append((-invere_base[MOD % i] * (MOD // i)) % MOD)
    inverse.append((inverse[-1] * invere_base[-1]) % MOD)

def nCr(n, r):
    if not (0 <= r <= n):
        return 0
    return factorial[n] * inverse[r] * inverse[n - r] % MOD

def nPr(n, r):
    if not (0 <= r <= n):
        return 0
    return nCr(n, r) * factorial[r] % MOD

mod = MOD
ans = 0

n = N
m = M

if N == M:
    for i in range(N+1):
        
        ans += nCr(n,i)*factorial[n]*factorial[n-i]*(-1)**(i)
        
        ans %= mod
        
    print(ans)
else:
    for i in range(N+1):
        ans += nCr(n,i)*nPr(m,n)*nPr(m-i,n-i)*(-1)**(i)
        ans %= mod
        
    print(ans)

Solution (Python):
class Comb():
    def __init__(self, N, mod=1):
        self.mod = mod
        fct = [1]
        invfct = [1]
        for i in range(1, N+1):
            fct.append(fct[i-1] * i % mod)
            invfct.append(pow(fct[i], mod-2, mod))
        self.fct = fct
        self.invfct = invfct
    
    def cmb(self, n, k):
        """ Calculate nCk = n! / (k! * (n-k)!) mod self.mod """
        return self.fct[n] * self.invfct[k] * self.invfct[n-k]

    def pam(self, n, k):
        """ Calculate nPk = n! / (n-k)! mod self.mod """
        return self.fct[n] * self.invfct[n-k]
N, M = map(int, input().split())
mod = 10**9 + 7
comb = Comb(N+M, mod)
ans = 0
for k in range(N+1):
    ans += pow(-1, k) * comb.cmb(N, k) * comb.pam(M, k) * pow(comb.pam(M-k, N- k), 2)
    ans %= mod
print(ans)

Solution (Python):
def nPk(n, k):
  return fact[n] * invfact[n-k] % MOD

def nCk(n, k):
  return fact[n] * invfact[k] * invfact[n-k] % MOD

N, M = map(int, input().split())
MOD = 10 ** 9 + 7

fact = [0] * 1100000
invfact = [0] * 1100000

fact[0] = 1

for i in range(1, 1100000):
  fact[i] = fact[i-1] * i % MOD

invfact[-1] = pow(fact[-1], MOD-2, MOD)

for i in range(1100000-1, 0, -1):
  invfact[i-1] = invfact[i] * i % MOD

ans = 0

for k in range(N + 1):
  temp = nCk(N, k) * nPk(M, k) * (nPk(M-k, N-k) ** 2) % MOD
  if k % 2 == 0:
    ans += temp
  else:
    ans -= temp
  ans %= MOD

print(ans)

Solution (Python):

mod = 10**9+7



def inv(a, mod):
    r = [1, 0, a]
    w = [0, 1, mod]
    while w[2] != 1:
        q = r[2]//w[2]
        r_new = [r[0]-q*w[0], r[1]-q*w[1], r[2]-q*w[2]]
        r = w
        w = r_new
    x, y = w[0], w[1]    
    return (mod+x % mod) % mod


max_num = 5*10**5+2  
fact = [0 for _ in range(max_num)]
ifact = [0 for _ in range(max_num)]

fact[0] = fact[1] = 1
ifact[0] = ifact[1] = 1

for i in range(2, max_num):
    fact[i] = fact[i-1] * i % mod

ifact[max_num-1] = inv(fact[max_num-1], mod)

for i in range(2, max_num):
    ifact[max_num-i] = (ifact[max_num-i+1] * (max_num-i+1)) % mod


def comb(x, y):
    if x < 0 or y < 0 or x < y:
        return 0
    else:
        return fact[x] * ifact[y] * ifact[x-y] % mod



def solve():
    tmp = 0
    for i in range(1, n+1):
        tmp2 = comb(n, i) * fact[m-i] * ifact[m-n]
        tmp2 %= mod
        tmp += (-1)**(i-1)*tmp2

    tmp = fact[m]*ifact[m-n] - tmp
    tmp %= mod

    ans = tmp*fact[m]*ifact[m-n]
    ans %= mod

    return ans


if __name__ == "__main__":
    n, m = map(int, input().split())
    print(solve())

Solution (Python):


def cmb(n,r,p):
  if r<0 or n<r:
    return 0
  r=min(r,n-r)
  return facto[n]*factoinv[r]*factoinv[n-r]%p
p=10**9+7
n,m=map(int,input().split())

facto=[1,1]

factoinv=[1,1]

inv=[0,1]

for i in range(2,m+1):
  facto.append((facto[-1]*i)%p)
  inv.append((-inv[p%i]*(p//i))%p)
  factoinv.append((factoinv[-1]*inv[-1])%p)
if n==m:
  kosuub=0
  for i in range(n+1):
    kosuub+=(facto[(n-i)]*cmb(n,i,p)*((-1)**i))%p
  kosuu=(kosuub*facto[(n)])%p
  print(kosuu)
else:
  kosuub=0
  for i in range(n+1):
    kosuub+=(cmb((m-i),(n-i),p)*facto[(n-i)]*cmb(n,i,p)*((-1)**i))%p
  kosuu=(kosuub*cmb(m,n,p)*facto[n])%p
  print(kosuu)

Solution (Python):
import sys
input = sys.stdin.readline

MAX = 5*10**5+100
MOD = 10**9+7
fact = [0]*MAX 
inv = [0]*MAX 
finv = [0]*MAX 
fact[0] = 1
fact[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1
    
for i in range(2, MAX):
    fact[i] = fact[i-1]*i%MOD
    inv[i] = MOD-inv[MOD%i]*(MOD//i)%MOD
    finv[i] = finv[i-1]*inv[i]%MOD

def C(n, r):
    if n<r:
        return 0
    if n<0 or r<0:
        return 0
    return fact[n]*(finv[r]*finv[n-r]%MOD)%MOD

N, M = map(int, input().split())
ans = (C(M, N)*fact[N])**2%MOD

for i in range(1, N+1):
    if i%2==0:
        ans += C(N, i)*C(M, i)*fact[i]*(C(M-i, N-i)*fact[N-i])**2
    else:
        ans -= C(N, i)*C(M, i)*fact[i]*(C(M-i, N-i)*fact[N-i])**2
    
    ans %= MOD

print(ans)

Solution (Python):
mod = 1000000007
fact = []
fact_inv = []
pow_ = []

def pow_ini(nn):
	global pow_
	pow_.append(nn)
	for j in range(62):
		nxt = pow_[j] * pow_[j]
		nxt %= mod
		pow_.append(nxt)
	return

def pow1(k):
    ansk = 1
    for ntmp in range(62):
        if((k >> ntmp) % 2 == 1):
            ansk *= pow_[ntmp]
            ansk %= mod
    return ansk

def fact_ini(n):
	global fact
	global fact_inv
	for i in range(n + 1):
		fact.append(0)
		fact_inv.append(0)
	fact[0] = 1
	for i in range(1, n + 1, 1):
		fact_tmp = fact[i - 1] * i
		fact_tmp %= mod
		fact[i] = fact_tmp
	pow_ini(fact[n])
	fact_inv[n] = pow1(mod - 2)
	for i in range(n - 1, -1, -1):
		fact_inv[i] = fact_inv[i + 1] * (i + 1)
		fact_inv[i] %= mod
	return

def nCm(n, m):
	assert(m >= 0)
	assert(n >= m)
	ans = fact[n] * fact_inv[m]
	ans %= mod
	ans *= fact_inv[n - m]
	ans %= mod
	return ans;

def nPm(n, m):
	assert(m >= 0)
	assert(n >= m)
	ans = fact[n]
	ans *= fact_inv[n - m]
	ans %= mod
	return ans;

fact_ini(5500000)
N, M = (int(x) for x in input().split())
ans = 0
for i in range(0, N+1, 1):
	anstmp = nCm(N, i) * nPm(M-i, N-i)
	anstmp %= mod
	if(i%2 == 1):
		anstmp *= -1
		anstmp += mod
		anstmp %= mod
	ans += anstmp
	ans %= mod
ans *= nPm(M, N)
ans %= mod
print(ans)

Solution (Python):
n, m = map(int, input().split())
mod = 10**9 + 7

factorial = [1]
inverse = [1]
for i in range(1, m+1):
    factorial.append(factorial[-1] * i % mod)
    inverse.append(pow(factorial[-1], mod - 2, mod))
def comb(n, r, mod):
    if n < r or r < 0: return 0
    elif r == 0: return 1
    return factorial[n] * inverse[r] * inverse[n - r] % mod
def perm(n, r, mod):
    if n < r or r < 0: return 0
    elif r == 0: return 1
    return factorial[n] * inverse[n - r] % mod

ans = 0
for i in range(n+1):
    ans += (-1)**i * comb(n, i, mod) * perm(m, i, mod) * (perm(m-i, n-i, mod))**2
    ans %= mod
print(ans)

Solution (Python):
import sys
input=sys.stdin.readline
sys.setrecursionlimit(10 ** 8)
from itertools import accumulate
from itertools import permutations
from itertools import combinations
from collections import defaultdict
from collections import Counter
import fractions
import math
from collections import deque
from bisect import bisect_left
from bisect import bisect_right
from bisect import insort_left
import itertools
from heapq import heapify
from heapq import heappop
from heapq import heappush
import heapq
from copy import deepcopy
from decimal import Decimal
alf = list("abcdefghijklmnopqrstuvwxyz")
ALF = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")

INF = float("inf")


N,M = map(int,input().split())
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
MOD = 10**9+7
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
for i in range(2, M+1):
    fac.append(fac[-1] * i % MOD)
    inv.append(MOD - inv[MOD%i] * (MOD//i) % MOD)
    finv.append(finv[-1] * inv[-1] % MOD) 
def comb_mod(n, r, m):  
    if (n<0 or r<0 or n<r): return 0
    r = min(r, n-r)
    return fac[n] * finv[n-r] * finv[r] % m

def chofuku(n,r): 
    return comb_mod(n+r-1,r)
 
def nPr(n,r,m): 
    return fac[n]*finv[n-r]%m
ans = 0
for i in range(1,N+1):
    if i % 2 == 1:
        flag = 1
    else:
        flag = -1
    s = flag*comb_mod(N,i,MOD)*nPr(M-i,N-i,MOD)
    ans += s
    if ans > MOD:
        ans %= MOD
ans1 = 1
for i in range(N):
    ans1 *= (M-i)
    ans1 %= MOD
ans2 = (ans1*ans1)%MOD
ans3 = (ans*nPr(M,N,MOD)) %MOD
print((ans2-ans3)%MOD)

Solution (Python):
import sys
input=sys.stdin.readline
sys.setrecursionlimit(10 ** 8)
from itertools import accumulate
from itertools import permutations
from itertools import combinations
from collections import defaultdict
from collections import Counter
import fractions
import math
from collections import deque
from bisect import bisect_left
from bisect import bisect_right
from bisect import insort_left
import itertools
from heapq import heapify
from heapq import heappop
from heapq import heappush
import heapq
from copy import deepcopy
from decimal import Decimal
alf = list("abcdefghijklmnopqrstuvwxyz")
ALF = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")

INF = float("inf")


N,M = map(int,input().split())
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
MOD = 10**9+7
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
for i in range(2, M+1):
    fac.append(fac[-1] * i % MOD)
    inv.append(MOD - inv[MOD%i] * (MOD//i) % MOD)
    finv.append(finv[-1] * inv[-1] % MOD) 
def comb_mod(n, r, m):  
    if (n<0 or r<0 or n<r): return 0
    r = min(r, n-r)
    return fac[n] * finv[n-r] * finv[r] % m

def chofuku(n,r): 
    return comb_mod(n+r-1,r)
 
def nPr(n,r,m): 
    return fac[n]*finv[n-r]%m
ans = 0
for i in range(1,N+1):
    if i % 2 == 1:
        flag = 1
    else:
        flag = -1
    s = flag*comb_mod(N,i,MOD)*nPr(M-i,N-i,MOD)
    ans += s
    if ans > MOD:
        ans %= MOD
ans1 = 1
for i in range(N):
    ans1 *= (M-i)
    ans1 %= MOD
ans2 = (ans1*ans1)%MOD
ans3 = (ans*nPr(M,N,MOD)) %MOD
print((ans2-ans3)%MOD)

Solution (Python):
fac=[1]*505050
inv=[1]*505050
finv=[1]*505050
mod=1000000007
for i in range(2,500010):
  fac[i]=fac[i-1]*i%mod
  inv[i]=-inv[mod%i]*(mod//i)%mod
  finv[i]=finv[i-1]*inv[i]%mod

def C(n,r):
  if n<0 or r<0 or n<r: return 0
  return fac[n]*finv[n-r]%mod*finv[r]%mod

N,M=map(int,input().split())
ans=0
for i in range(N+1):
  ans+=C(N,i)*fac[M]%mod*finv[M-i]%mod*((fac[M-i]*finv[M-N]%mod)**2%mod)*(1 if i%2==0 else -1)%mod
  if ans<0: ans+=mod
  if ans>=mod: ans-=mod
print(ans)

Solution (Python):
fac=[1]*505050
inv=[1]*505050
finv=[1]*505050
mod=1000000007
for i in range(2,500010):
  fac[i]=fac[i-1]*i%mod
  inv[i]=-inv[mod%i]*(mod//i)%mod
  finv[i]=finv[i-1]*inv[i]%mod

def C(n,r):
  if n<0 or r<0 or n<r: return 0
  return fac[n]*finv[n-r]%mod*finv[r]%mod

N,M=map(int,input().split())
ans=0
for i in range(N+1):
  ans+=C(N,i)*fac[M]%mod*finv[M-i]%mod*((fac[M-i]*finv[M-N]%mod)**2%mod)*(1 if i%2==0 else -1)%mod
  ans%=mod
print(ans%mod)

Solution (Python):
fac=[1]*505050
inv=[1]*505050
finv=[1]*505050
mod=1000000007
for i in range(2,500010):
  fac[i]=fac[i-1]*i%mod
  inv[i]=-inv[mod%i]*(mod//i)%mod
  finv[i]=finv[i-1]*inv[i]%mod

def C(n,r):
  if n<0 or r<0 or n<r: return 0
  return fac[n]*finv[n-r]%mod*finv[r]%mod

N,M=map(int,input().split())
ans=0
for i in range(N+1):
  ans+=C(N,i)*fac[M]%mod*finv[M-i]%mod*((fac[M-i]*finv[M-N])**2%mod)*(1 if i%2==0 else -1)%mod
  ans%=mod
print(ans%mod)

Solution (Python):
fac=[1]*505050
inv=[1]*505050
finv=[1]*505050
mod=1000000007
for i in range(2,500010):
  fac[i]=fac[i-1]*i%mod
  inv[i]=-inv[mod%i]*(mod//i)%mod
  finv[i]=finv[i-1]*inv[i]%mod

def C(n,r):
  if n<0 or r<0 or n<r: return 0
  return fac[n]*finv[n-r]*finv[r]%mod

N,M=map(int,input().split())
ans=0
for i in range(N+1):
  ans+=C(N,i)*fac[M]%mod*finv[M-i]%mod*((fac[M-i]*finv[M-N])**2%mod)*(1 if i%2==0 else -1)%mod
  ans%=mod
print(ans%mod)

Solution (Python):
fac=[1]*505050
inv=[1]*505050
finv=[1]*505050
mod=1000000007
for i in range(2,500010):
  fac[i]=fac[i-1]*i%mod
  inv[i]=-inv[mod%i]*(mod//i)%mod
  finv[i]=finv[i-1]*inv[i]%mod

def C(n,r):
  if n<0 or r<0 or n<r: return 0
  return fac[n]*finv[n-r]*finv[r]%mod

N,M=map(int,input().split())
ans=0
for i in range(N+1):
  ans+=C(N,i)*fac[M]*finv[M-i]*((fac[M-i]*finv[M-N])**2)*(1 if i%2==0 else -1) % mod
print(ans%mod)

Solution (Python):
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 = 10**9+7 
N = 5*10**5 
g1 = [0]*(N+1) 
g1[0] = g1[1] = 1
g2 = [0]*(N+1) 
g2[0] = g2[1] = 1
inverse = [0]*(N+1) 
inverse[0],inverse[1] = 0,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

N, M = map(int, input().split())
ans = 0
for i in range(N+1):
    ans += cmb(N,i,mod)*((-1)**i)*g1[M]*g2[M-i]*((g1[M-i]*g2[M-N])**2)
    ans %= mod
print(ans)

Solution (Python):
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 = 10**9+7 
N = 5*10**5 
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 
for i in range( 2, N + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )

N, M = map(int, input().split())
ans = 0
for i in range(N+1):
    ans += cmb(N,i,mod)*((-1)**i)*g1[M]*g2[M-i]*((g1[M-i]*g2[M-N])**2)
    ans %= mod
print(ans)

Solution (Python):
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 = 10**9+7 
N = 5*10**5
g1 = [1, 1] 
g2 = [1, 1] 
inverse = [0, 1] 
for i in range( 2, N + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )

N, M = map(int, input().split())
ans = 0
for i in range(N+1):
    ans += cmb(N,i,mod)*((-1)**i)*g1[M]*g2[M-i]*((g1[M-i]*g2[M-N])**2)
    ans %= mod
print(ans)

Solution (Python):
import sys
input=sys.stdin.readline
sys.setrecursionlimit(10 ** 8)
from itertools import accumulate
from itertools import permutations
from itertools import combinations
from collections import defaultdict
from collections import Counter
import fractions
import math
from collections import deque
from bisect import bisect_left
from bisect import bisect_right
from bisect import insort_left
import itertools
from heapq import heapify
from heapq import heappop
from heapq import heappush
import heapq
from copy import deepcopy
from decimal import Decimal
alf = list("abcdefghijklmnopqrstuvwxyz")
ALF = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")

INF = float("inf")


N,M = map(int,input().split())
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
MOD = 10**9+7
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
for i in range(2, M+1):
    fac.append(fac[-1] * i % MOD)
    inv.append(MOD - inv[MOD%i] * (MOD//i) % MOD)
    finv.append(finv[-1] * inv[-1] % MOD) 
def comb_mod(n, r, m):  
    if (n<0 or r<0 or n<r): return 0
    r = min(r, n-r)
    return fac[n] * finv[n-r] * finv[r] % m

def chofuku(n,r): 
    return comb_mod(n+r-1,r)
 
def nPr(n,r,m): 
    return fac[n]*finv[n-r]%m
ans = 0
for i in range(1,N+1):
    if i % 2 == 1:
        flag = 1
    else:
        flag = -1
    s = flag*comb_mod(N,i,MOD)*nPr(M-i,N-i,MOD)
    ans += s
    if ans > MOD:
        ans %= MOD
ans1 = 1
for i in range(N):
    ans1 *= (M-i)
    ans1 %= MOD
ans2 = (ans1*ans1)%MOD
ans3 = (ans*nPr(M,N,MOD)) %MOD
print((ans2-ans3)%MOD)

Solution (Python):
import sys
input = sys.stdin.readline

N,M=map(int,input().split())
mod=10**9+7

FACT=[1]
for i in range(1,5*10**5+1):
    FACT.append(FACT[-1]*i%mod)

FACT_INV=[pow(FACT[-1],mod-2,mod)]
for i in range(5*10**5,0,-1):
    FACT_INV.append(FACT_INV[-1]*i%mod)

FACT_INV.reverse()

def Combi(a,b):
    if 0<=b<=a:
        return FACT[a]*FACT_INV[b]*FACT_INV[a-b]%mod
    else:
        return 0

ANS=0
for i in range(N+1):
    if i%2==0:
        ANS=(ANS+Combi(N,i)*FACT[M-i]*FACT_INV[M-N])%mod
    else:
        ANS=(ANS-Combi(N,i)*FACT[M-i]*FACT_INV[M-N])%mod

print(ANS*FACT[M]*FACT_INV[M-N]%mod)

Solution (Python):
n,m = map(int,input().split())
mod = 10**9+7
rng = 500100
fctr = [1]
finv = [1]
for i in range(1,rng):
  fctr.append(fctr[-1]*i%mod)
for i in range(1,rng):
  finv.append(pow(fctr[i],mod-2,mod))
def cmb(n,k):
  if n<0 or k<0:
    return 0
  else:
    return fctr[n]*finv[n-k]*finv[k]%mod
ans = 0
for i in range(n+1):
  if i%2 == 0:
    ans += fctr[m-i]*finv[m-n]*cmb(n,i)
  else:
    ans -= fctr[m-i]*finv[m-n]*cmb(n,i)
  ans %= mod
print(ans*fctr[m]*finv[m-n]%mod)

Solution (Python):
MOD = 10**9 + 7
MAX = int(5e5)


def div(a, b):
    return a * pow(b, MOD-2, MOD) % MOD


FACT = [1] * (MAX+1)
for i in range(1, MAX+1):
    FACT[i] = (i * FACT[i-1]) % MOD
INV = [1] * (MAX+1)
INV[MAX] = div(1, FACT[MAX])
for i in range(MAX, 0, -1):
    INV[i-1] = (INV[i] * i) % MOD


def combination(m, n):
    return (FACT[m] * INV[n] * INV[m-n]) % MOD


def permutation(m, n):
    return (FACT[m] * INV[m-n]) % MOD


def main():
    N, M = map(int, input().split())
    ans = 0
    for i in range(N+1):
        
        tmp = (combination(N, i) * permutation(M, i)
               * (permutation(M-i, N-i)**2)) % MOD
        ans = (ans + (-1)**i * tmp) % MOD

    print(ans)


if __name__ == "__main__":
    main()

Solution (Python):
def cmb(n, r, p):
    if (r < 0) or (n < r):
        return 0
    r = min(r, n - r)
    return fact[n] * factinv[r] * factinv[n-r] % p

def per(n, r, p):
    if r < 0 or n < r:
        return 0
    return fact[n] * factinv[n-r] % p

p = 10 ** 9 + 7
N = 5*pow(10, 5)  
fact = [1, 1]
factinv = [1, 1]
inv = [0, 1]
 
for i in range(2, N + 1):
    fact.append((fact[-1] * i) % p)
    inv.append((-inv[p % i] * (p // i)) % p)
    factinv.append((factinv[-1] * inv[-1]) % p)

n, m = map(int, input().split())
ans = 0
for i in range(n+1):
    ans += pow(per(m-i, n-i, p), 2) * per(m, i, p) * cmb(n, i, p) * (-1)**i % p
ans = ans % p
print(ans)

Solution (Python):
def cmb(n, r, p):
    if (r < 0) or (n < r):
        return 0
    r = min(r, n - r)
    return fact[n] * factinv[r] * factinv[n-r] % p

def per(n, r, p):
    if r < 0 or n < r:
        return 0
    return fact[n] * factinv[n-r] % p

p = 10 ** 9 + 7
N = 5*pow(10, 5)  
fact = [1, 1]
factinv = [1, 1]
inv = [0, 1]
 
for i in range(2, N + 1):
    fact.append((fact[-1] * i) % p)
    inv.append((-inv[p % i] * (p // i)) % p)
    factinv.append((factinv[-1] * inv[-1]) % p)

n, m = map(int, input().split())
ans = 0
for i in range(n+1):
    ans += pow(per(m-i, n-i, p), 2) * per(m, i, p) * cmb(n, i, p) * (-1)**i % p
ans = ans % p
print(ans)

Solution (Python):
import sys
input=sys.stdin.readline
sys.setrecursionlimit(10 ** 8)
from itertools import accumulate
from itertools import permutations
from itertools import combinations
from collections import defaultdict
from collections import Counter
import fractions
import math
from collections import deque
from bisect import bisect_left
from bisect import bisect_right
from bisect import insort_left
import itertools
from heapq import heapify
from heapq import heappop
from heapq import heappush
import heapq
from copy import deepcopy
from decimal import Decimal
alf = list("abcdefghijklmnopqrstuvwxyz")
ALF = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")

INF = float("inf")


N,M = map(int,input().split())
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
MOD = 10**9+7
for i in range(2, M+2):
    fac.append(fac[-1] * i % MOD)
    inv.append(MOD - inv[MOD%i] * (MOD//i) % MOD)
    finv.append(finv[-1] * inv[-1] % MOD) 
def comb_mod(n, r, m):  
    if (n<0 or r<0 or n<r): return 0
    r = min(r, n-r)
    return fac[n] * finv[n-r] * finv[r] % m

def chofuku(n,r): 
    return comb_mod(n+r-1,r)
 
def nPr(n,r,m): 
    return fac[n]*finv[n-r]%m
ans = 0
for i in range(1,N+1):
    if i % 2 == 1:
        flag = 1
    else:
        flag = -1
    s = flag*comb_mod(N,i,MOD)*nPr(M-i,N-i,MOD)
    ans += s
    if ans > MOD:
        ans %= MOD
ans1 = 1
for i in range(N):
    ans1 *= (M-i)
    ans1 %= MOD
ans2 = (ans1*ans1)%MOD
ans3 = (ans*nPr(M,N,MOD)) %MOD
print((ans2-ans3)%MOD)

Solution (Python):
import sys
input=sys.stdin.readline
sys.setrecursionlimit(10 ** 8)
from itertools import accumulate
from itertools import permutations
from itertools import combinations
from collections import defaultdict
from collections import Counter
import fractions
import math
from collections import deque
from bisect import bisect_left
from bisect import bisect_right
from bisect import insort_left
import itertools
from heapq import heapify
from heapq import heappop
from heapq import heappush
import heapq
from copy import deepcopy
from decimal import Decimal
alf = list("abcdefghijklmnopqrstuvwxyz")
ALF = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ")

INF = float("inf")


N,M = map(int,input().split())
fac = [1, 1]
inv = [0, 1]
finv = [1, 1]
MOD = 10**9+7
for i in range(2, M+2):
    fac.append(fac[-1] * i % MOD)
    inv.append(MOD - inv[MOD%i] * (MOD//i) % MOD)
    finv.append(finv[-1] * inv[-1] % MOD) 
def comb_mod(n, r, m):  
    if (n<0 or r<0 or n<r): return 0
    r = min(r, n-r)
    return fac[n] * finv[n-r] * finv[r] % m

def chofuku(n,r): 
    return comb_mod(n+r-1,r)
 
def nPr(n,r,m): 
    return fac[n]*finv[n-r]%m
ans = 0
for i in range(1,N+1):
    if i % 2 == 1:
        flag = 1
    else:
        flag = -1
    s = flag*comb_mod(N,i,MOD)*nPr(M-i,N-i,MOD)
    ans += s
    if ans > MOD:
        ans %= MOD
ans1 = 1
for i in range(N):
    ans1 *= (M-i)
    ans1 %= MOD
ans2 = (ans1*ans1)%MOD
ans3 = (ans*nPr(M,N,MOD)) %MOD
print((ans2-ans3)%MOD)