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):
import sys
sys.setrecursionlimit(10**6)

MOD = 10**9+7


def make_cmb_table(n):
    fact = [1]
    for i in range(1, n+1):
        fact.append(fact[i-1] * i % MOD)

    ifact = [0] * (n+1)
    ifact[n] = pow(fact[n], MOD-2, MOD)
    for i in range(n, 0, -1):
        ifact[i-1] = ifact[i] * i % MOD
    return (fact, ifact)


def combination(n, r, table):
    if r > n:
        return 0
    return table[0][n] * table[1][n - r] * table[1][r] % MOD


def permutation(n, k, table):
    return table[0][n] * table[1][n - k] % MOD


def main(input, print):
    n, m = map(int, input().split())
    ans = 0
    table = make_cmb_table(10**6+5)
    for i in range(n+1):
        now = combination(n, i, table)
        now %= MOD
        now *= (-1) ** i
        now *= permutation(m, i, table)
        now %= MOD
        now *= permutation(m-i, n-i, table) ** 2
        now %= MOD
        ans += now
        ans %= MOD
    print(ans)


if __name__ == '__main__':
    main(sys.stdin.readline, print)

Solution (Python):
import sys

readline = sys.stdin.readline
readall = sys.stdin.read
ns = lambda: readline().rstrip()
ni = lambda: int(readline().rstrip())
nm = lambda: map(int, readline().split())
nl = lambda: list(map(int, readline().split()))
prn = lambda x: print(*x, sep='\n')

n_ = 5 * 10**5
mod = 10**9 + 7
fun = [1] * (n_ + 1)
for i in range(1, n_ + 1):
    fun[i] = fun[i - 1] * i % mod
rev = [1] * (n_ + 1)
rev[n_] = pow(fun[n_], mod - 2, mod)
for i in range(n_ - 1, 0, -1):
    rev[i] = rev[i + 1] * (i + 1) % mod

def nCr(n, r):
    if r > n:
        return 0
    return fun[n] * rev[r] % mod * rev[n - r] % mod

def modinv(x, mod):
    a, b = x, mod
    u, v = 1, 0
    while b:
        t = a // b
        a -= t * b; a, b = b, a
        u -= t * v; u, v = v, u
    return u % mod

def solve():
    n, m = nm()
    s = 0
    v = 1
    for i in range(n+1):
        s = (s + v * nCr(n, i) * nCr(m, i) * fun[i] % mod * (nCr(m-i, n-i) * fun[n-i] % mod)**2) % mod
        v *= -1
    print(s)
    return

solve()

Solution (Python):
N, M = map(int, input().split())
mod = 10 ** 9 + 7
fact = [1, 1]
factinv = [1, 1]
inv = [0, 1]


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

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


for i in range(2, M + 1):
    fact.append((fact[-1] * i) % mod)
    inv.append((-inv[mod % i] * (mod // i)) % mod)
    factinv.append((factinv[-1] * inv[-1]) % mod)


ans = 0
for i in range(N+1):
    ans += cmb(N, i) * perm(M - i, N - i) * (-1) ** i
    ans %= mod
ans *= perm(M, N)
ans %= mod

print(ans)

Solution (Python):

from sys import stdin
input = stdin.readline
printn = lambda x: print(x,end='')
inn = lambda : int(input())
inl   = lambda: list(map(int, input().split()))
inm   = lambda:      map(int, input().split())
ins = lambda : input().strip()
DBG = True 
BIG = 10**18
R = 10**9 + 7

def ddprint(x):
  if DBG:
    print(x)

def modinv2(x,r):
  return pow(x,r-2,r)

def perm(n,a,R):
    return (fact[n]*modinv2(fact[n-a],R))

def comb(n,a,R):
    return (fact[n]*modinv2(fact[n-a]*fact[a],R))

def factinit(n,R):
    global fact
    fact = [1]*(n+1)
    for i in range(2,n+1):
        fact[i] = (fact[i-1]*i)%R

n,m = inm()
factinit(m,R)
x = 0
for k in range(n+1):
    inc = comb(n,k,R)*((-1)**k)*perm(m,k,R)* \
          pow(perm(m-k,n-k,R),2,R)
    inc %= R
    x = (x+inc)%R
    
print(x%R)

Solution (Python):

import sys
sys.setrecursionlimit(10**6)
MOD = 10**9+7
mon_l = dict({})

def mon(m,n):
    if (m,n) in mon_l:
        return mon_l[(m,n)]
    if n==1:
        mon_l[(m,n)] = m-1
        return mon_l[(m,n)]
    if n==2:
        mon_l[(m,n)] = ((m-1)+(m-2)**2)%MOD
        return mon_l[(m,n)]
    else:
        mon_l[(m,n)] = ((n-1)*(mon(m-1,n-1)+mon(m-2,n-2)) + (m-n)*mon(m-1,n-1))%MOD
        return mon_l[(m,n)]
    
def permu(m,n):
    ans = 1
    for _ in range(n):
        ans = m*ans%MOD
        m -= 1
    return ans

N,M = map(int,input().split())
ans = permu(M,N)*mon(M,N)%MOD
print(ans)

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

fact = [1,1]
finv = [1,1]
inv = [0,1]
 
for i in range(2,max(n,m)+5):
    fact.append((fact[-1]*i)%mod)
    inv.append((inv[mod%i]*(mod-mod//i))%mod)
    finv.append((finv[-1]*inv[-1])%mod)
 
def nCr(n,r,mod):
    if r > n:
        return 0
    else: 
        return fact[n]*finv[r]*finv[n-r]%mod

def nPr(n,r,mod):
    if r > n:
        return 0
    else: 
        return fact[n]*finv[n-r]%mod


ans = 0
for i in range(n+1):
    now = nCr(n,i,mod)
    now *= nPr(m-i,n-i,mod)
    now %= mod
    if i%2:
        now = -now
    ans += now
    ans % mod
ans *= nPr(m,n,mod)
ans %= mod
print(ans)

Solution (Python):

MAX_NUM = 5* 10**5 + 1
MOD = 10**9+7

fac  = [0 for _ in range(MAX_NUM)]
finv = [0 for _ in range(MAX_NUM)]
inv  = [0 for _ in range(MAX_NUM)]

fac[0]  = fac[1] = 1
finv[0] = finv[1] = 1
inv[1] = 1

for i in range(2,MAX_NUM):
  fac[i] = fac[i-1] * i % MOD
  inv[i] = MOD - inv[MOD%i] * (MOD // i) % MOD
  finv[i] = finv[i-1] * inv[i] % MOD

def comb(n,k):
  if n < k:
    return 0
  if n < 0 or k < 0:
    return 0
  return fac[n] * (finv[k] * finv[n-k] % MOD) % MOD

def perm(n,k):
  if n < k:
    return 0
  if n < 0 or k < 0:
    return 0
  return fac[n] * finv[n-k] % MOD


N, M = map(int,input().split())

ans = 0
for K in range(N+1):
  term = comb(N, K) * perm(M, K) 
  term *= perm(M - K, N - K)**2
  term %= MOD
  if K % 2 == 0:
    ans = (ans + term) % MOD
  else:
    ans = (ans - term) % MOD
    
print(ans)

Solution (Python):
MOD = 10**9 + 7

class Combinatorics:
    def __init__(self, N, mod):
        '''
        Preprocess for calculating binomial coefficients nCr (0 <= r <= n, 0 <= n <= N)
        over the finite field Z/(mod)Z.
        Input:
            N (int): maximum n
            mod (int): a prime number. The order of the field Z/(mod)Z over which nCr is calculated.
        '''
        self.mod = mod
        self.fact = [0] * (N+1)           
        self.inverse = [None] + [0] * N   
        self.fact_inverse = [0] * (N+1)   
        
        
        self.fact[0] = self.fact[1] = 1
        self.fact_inverse[0] = self.fact_inverse[1] = 1
        self.inverse[1] = 1
        for i in range(2, N+1):
            self.fact[i] = i * self.fact[i-1] % self.mod
            q, r = divmod(self.mod, i)
            self.inverse[i] = (- (q % self.mod) * self.inverse[r]) % self.mod
            self.fact_inverse[i] = self.inverse[i] * self.fact_inverse[i-1] % self.mod
    
    def perm(self, n, r):
        '''
        Calculate nPr = n! / (n-r)! % mod
        '''
        if n < r or n < 0 or r < 0: return 0
        else: return (self.fact[n] * self.fact_inverse[n-r]) % self.mod
    
    def binom(self, n, r):
        '''
        Calculate nCr = n! /(r! (n-r)!) % mod
        '''
        if n < r or n < 0 or r < 0: return 0
        else: return self.fact[n] * (self.fact_inverse[r] * self.fact_inverse[n-r] % self.mod) % self.mod
        
    def hom(self, n, r):
        '''
        Calculate nHr = {n+r-1}Cr % mod.
        Assign r objects to one of n classes.
        Arrangement of r circles and n-1 partitions:
            o o o | o o | | | o | | | o o | | o
        '''
        if n == 0 and r > 0: return 0
        if n >= 0 and r == 0: return 1
        return self.binom(n + r - 1, r)

N, M = map(int, input().split())
com = Combinatorics(M, MOD)
ans = 0
for i in range(N+1):
    temp = com.perm(M, N)
    temp *= com.binom(N, i); temp %= MOD
    temp *= com.perm(M-i, N-i); temp %= MOD
    ans += ((-1)**i) * temp; ans %= MOD
print(ans)

Solution (Python):


MOD = 10**9 + 7
m = 5 * 10**5 + 500
fac = [0] * m
finv = [0] * m
inv = [0] * m


def COMBinitialize(m):
    fac[0] = 1
    finv[0] = 1
    if m > 1:
        fac[1] = 1
        finv[1] = 1
        inv[1] = 1
        for i in range(2, m):
            fac[i] = fac[i-1] * i % MOD
            inv[i] = MOD - inv[MOD % i] * (MOD // i) % MOD
            finv[i] = finv[i - 1] * inv[i] % MOD


def COMB(n, k):
    if n < k:
        return 0
    if n < 0 or k < 0:
        return 0
    return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD


def PERM(n, k):
    if n < k:
        return 0
    if n < 0 or k < 0:
        return 0
    return fac[n] * (finv[n - k] % MOD) % MOD


COMBinitialize(m)


def main():
    N, M = (int(i) for i in input().split())
    ans = 0
    for k in range(N+1):
        ans += PERM(M-k, N-k) * COMB(N, k) * (-1)**(k % 2)
        ans %= MOD
    print(ans * PERM(M, N) % MOD)


if __name__ == '__main__':
    main()

Solution (Python):


MOD = 10**9 + 7
m = 5 * 10**5 + 500
fac = [0] * m
finv = [0] * m
inv = [0] * m


def COMBinitialize(m):
    fac[0] = 1
    finv[0] = 1
    if m > 1:
        fac[1] = 1
        finv[1] = 1
        inv[1] = 1
        for i in range(2, m):
            fac[i] = fac[i-1] * i % MOD
            inv[i] = MOD - inv[MOD % i] * (MOD // i) % MOD
            finv[i] = finv[i - 1] * inv[i] % MOD


def COMB(n, k):
    if n < k:
        return 0
    if n < 0 or k < 0:
        return 0
    return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD


def PERM(n, k):
    if n < k:
        return 0
    if n < 0 or k < 0:
        return 0
    return fac[n] * (finv[n - k] % MOD) % MOD


COMBinitialize(m)


def main():
    N, M = (int(i) for i in input().split())
    if N == M == 1:
        return print(0)
    ans = 0
    for k in range(N+1):
        ans += PERM(M-k, N-k) * COMB(N, k) * (-1)**(k % 2)
        ans %= MOD
    print(ans * PERM(M, N) % MOD)


if __name__ == '__main__':
    main()

Solution (Python):
class Factorials:
    def __init__(self, n, mod):
        self.mod = mod

        
        self.fac = [1]
        num = 1
        for i in range(1, n+1):
            num *= i
            num %= mod
            self.fac.append(num)

        
        num = pow(num, mod-2, mod)
        self.rec = [1 for i in range(n+1)]
        self.rec[n] = num
        for i in range(n-1, 0, -1):
            num *= i + 1
            num %= mod
            self.rec[i] = num

    
    def comb(self, n, r):
        return self.fac[n] * self.rec[r] * self.rec[n - r] % self.mod
    
    
    def perm(self, n, r):
        return self.fac[n] * self.rec[n-r] % self.mod

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

    f = Factorials(m, mod)

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

    print(ans)

main()

Solution (Python):
class Factorials:
    def __init__(self, n, mod):
        self.mod = mod

        
        self.fac = [1]
        num = 1
        for i in range(1, n+1):
            num *= i
            num %= mod
            self.fac.append(num)

        
        num = pow(num, mod-2, mod)
        self.rec = [1 for i in range(n+1)]
        self.rec[n] = num
        for i in range(n-1, 0, -1):
            num *= i + 1
            num %= mod
            self.rec[i] = num

    
    def comb(self, n, r):
        return self.fac[n] * self.rec[r] * self.rec[n - r] % self.mod
    
    
    def perm(self, n, r):
        return self.fac[n] * self.rec[n-r] % self.mod

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

    f = Factorials(m, mod)

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

    print(ans)

main()

Solution (Python):
class Factorials:
    def __init__(self, n, mod):
        self.mod = mod

        
        self.fac = [1]
        num = 1
        for i in range(1, n+1):
            num *= i
            num %= mod
            self.fac.append(num)

        
        num = pow(num, mod-2, mod)
        self.rec = [1 for i in range(n+1)]
        self.rec[n] = num
        for i in range(n-1, 0, -1):
            num *= i + 1
            num %= mod
            self.rec[i] = num

    
    def comb(self, n, r):
        return self.fac[n] * self.rec[r] * self.rec[n - r] % self.mod
    
    
    def perm(self, n, r):
        return self.fac[n] * self.rec[n-r] % self.mod

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

f = Factorials(m, mod)

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

print(ans)

Solution (Python):

import sys
sys.setrecursionlimit(10**8)
def input(): return sys.stdin.readline().strip()
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(): return list(map(int, input().split()))
def ZIP(n): return zip(*(MAP() for _ in range(n)))


INF = float("inf")
MOD = 10**9 + 7

class Combination(object):
    def __init__(self, N, mod=MOD):
        fac, finv, inv = [0] * (N + 1), [0] * (N + 1), [0] * (N + 1)
        fac[:2] = 1, 1
        finv[:2] = 1, 1
        inv[1] = 1
        for i in range(2, N + 1):
            fac[i] = fac[i - 1] * i % mod
            inv[i] = -inv[mod % i] * (mod // i) % mod
            finv[i] = finv[i - 1] * inv[i] % mod
        self.N = N
        self.MOD = mod
        self.fac = fac
        self.finv = finv
        self.inv = inv

    def __call__(self, n, k):
        return self.C(n, k)

    def C(self, n, k):
        if n < k:
            return 0
        if n < 0 or k < 0:
            return 0
        fac = self.fac
        finv = self.finv
        mod = self.MOD
        return fac[n] * (finv[k] * finv[n - k] % mod) % mod

    def P(self, n, k):
        if n < k:
            return 0
        if n < 0 or k < 0:
            return 0
        return self.fac[n] * self.finv[n - k]


def main():
    N, M = MAP()
    cmb = Combination(M)

    ans = 0
    for i in range(N + 1):
        buf = (-1) ** i * cmb(N, i)
        buf %= MOD
        buf *= cmb.P(M, i)
        buf %= MOD
        buf *= pow(cmb.P(M - i, N - i), 2, MOD)
        buf %= MOD
        ans += buf
        ans %= MOD
    print(ans)

    return


if __name__ == '__main__':
    main()

Solution (Python):

import sys
sys.setrecursionlimit(10**8)
def input(): return sys.stdin.readline().strip()
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(): return list(map(int, input().split()))
def ZIP(n): return zip(*(MAP() for _ in range(n)))


INF = float("inf")
MOD = 10**9 + 7

class Combination(object):
    def __init__(self, N, mod=MOD):
        fac, finv, inv = [0] * (N + 1), [0] * (N + 1), [0] * (N + 1)
        fac[:2] = 1, 1
        finv[:2] = 1, 1
        inv[1] = 1
        for i in range(2, N + 1):
            fac[i] = fac[i - 1] * i % mod
            inv[i] = -inv[mod % i] * (mod // i) % mod
            finv[i] = finv[i - 1] * inv[i] % mod
        self.N = N
        self.MOD = mod
        self.fac = fac
        self.finv = finv
        self.inv = inv

    def __call__(self, n, k):
        return self.C(n, k)

    def C(self, n, k):
        if n < k:
            return 0
        if n < 0 or k < 0:
            return 0
        fac = self.fac
        finv = self.finv
        mod = self.MOD
        return fac[n] * (finv[k] * finv[n - k] % mod) % mod

    def P(self, n, k):
        if n < k:
            return 0
        if n < 0 or k < 0:
            return 0
        return self.fac[n] * self.finv[n - k]


def main():
    N, M = MAP()
    cmb = Combination(M)

    ans = 0
    for i in range(N + 1):
        buf = (-1) ** i * cmb(N, i)
        buf %= MOD
        buf *= cmb.P(M, i)
        buf %= MOD
        buf *= pow(cmb.P(M - i, N - i), 2, MOD)
        buf %= MOD
        ans += buf
        ans %= MOD
    print(ans)

    return


if __name__ == '__main__':
    main()

Solution (Python):
import sys,queue,math,copy,itertools,bisect,collections,heapq

def main():
    MOD = 10**9 + 7
    LI = lambda : [int(x) for x in sys.stdin.readline().split()]

    N,M = LI()

    f = [1]
    r = [1]
    c = 1
    for i in range(1, M + 1):
        c = (c * i) % MOD
        f.append(c)
        r.append(pow(c, MOD - 2, MOD))

    def comb(n, k):
        return (f[n] * r[k] * r[n - k]) % MOD
    def perm(n, k):
        return (f[n] * r[n-k]) % MOD


    ans = 0
    for i in range(N+1):
        ans = (ans + ((-1)**i) * comb(N,i) * perm(M,i) * (perm(M-i,N-i) ** 2)) % MOD
    print(ans)



if __name__ == '__main__':
    main()

Solution (Python):
MOD = 10**9+7
MAXN = 5*(10**5)+5
fact = [1]
ifact = [1]
for i in range(MAXN):
    fact.append(fact[-1] * (i+1) % MOD)
    ifact.append(pow(fact[-1], MOD-2, MOD))
def nCr(n, r):
    return fact[n] * ifact[r] * ifact[n-r] % MOD

def nPr(n, r):
    return fact[n] * ifact[n-r] % MOD

n, m = map(int, input().split())
ans = 0
for k in range(n+1):
    
    
    tmp = nPr(m-k, n-k) * nCr(n,k) % MOD
    if k%2:
        ans -= tmp
    else:
        ans += tmp
ans = (ans * nPr(m, n) + MOD) % MOD
print(ans)

Solution (Python):
MOD = 10**9+7
MAXN = 5*(10**5)+5
f = [1]
ifact = [1]
for i in range(MAXN):
    f.append(f[-1] * (i+1) % MOD)
    ifact.append(pow(f[-1], MOD-2, MOD))
def nCr(n, r):
    return f[n] * ifact[r] * ifact[n-r] % MOD

def nPr(n, r):
    return f[n] * ifact[n-r] % MOD

n, m = map(int, input().split())
ans = 0
for k in range(n+1):
    
    tmp = nPr(m-k, n-k) * nPr(m, n) * nCr(n,k) % MOD
    if k%2:
        ans -= tmp
    else:
        ans += tmp
ans = (ans + MOD) % MOD
print(ans)

Solution (Python):
MOD = 10**9+7
MAXN = 5*(10**5)+5
f = [1]
for i in range(MAXN):
    f.append(f[-1] * (i+1) % MOD)
def nCr(n, r, mod=MOD):
    return f[n] * pow(f[r], mod-2, mod) * pow(f[n-r], mod-2, mod) % mod

def nPr(n, r, mod=MOD):
    return f[n] * pow(f[n-r], mod-2, mod) % mod

n, m = map(int, input().split())
ans = 0
for k in range(n+1):
    
    tmp = nPr(m-k, n-k) * nPr(m, n) * nCr(n,k) % MOD
    if k%2:
        ans -= tmp
    else:
        ans += tmp
ans = (ans + MOD) % MOD
print(ans)

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

MAX = 5 * 10 ** 5 + 1
MOD = 10 ** 9 + 7


fac = [0] * (MAX + 1)
fac[0] = 1
fac[1] = 1
for i in range(2, MAX + 1):
    fac[i] = fac[i - 1] * i % MOD


finv = [0] * (MAX + 1)
finv[MAX] = pow(fac[MAX], MOD - 2, MOD)
for i in reversed(range(1, MAX + 1)):
    finv[i - 1] = finv[i] * i % MOD

def comb(n, k):
    if n < k or k < 0:
        return 0
    return (fac[n] * finv[k] * finv[n - k]) % MOD

def permute(n, k):
    return comb(n, k) * fac[k]

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

print(ans)

Solution (Python):
MOD = 10**9+7
MAXN = 5*(10**5)+5
f = [1]
for i in range(MAXN):
    f.append(f[-1] * (i+1) % MOD)
def nCr(n, r, mod=MOD):
    return f[n] * pow(f[r], mod-2, mod) * pow(f[n-r], mod-2, mod) % mod

def nPr(n, r, mod=MOD):
    return f[n] * pow(f[n-r], mod-2, mod) % mod

n, m = map(int, input().split())
ans = 0
for k in range(n+1):
    tmp = pow(nPr(m-k, n-k), 2, MOD) * nPr(m, k) * nCr(n,k) % MOD
    if k%2:
        ans -= tmp
    else:
        ans += tmp
ans = (ans + MOD) % MOD
print(ans)

Solution (Python):
MOD = 10**9+7
MAXN = 5*(10**5)+5
f = [1]
for i in range(MAXN):
    f.append(f[-1] * (i+1) % MOD)
def nCr(n, r, mod=MOD):
    return f[n] * pow(f[r], mod-2, mod) * pow(f[n-r], mod-2, mod) % mod

def nPr(n, r, mod=MOD):
    return f[n] * pow(f[n-r], mod-2, mod) % mod

n, m = map(int, input().split())
ans = 0
for k in range(n+1):
    tmp = nPr(m-k, n-k)
    tmp *= tmp
    tmp %= MOD
    tmp *= nPr(m, k)
    tmp %= MOD
    tmp *= nCr(n,k)
    tmp %= MOD
    if k%2:
        ans -= tmp
    else:
        ans += tmp
ans += MOD
ans %= MOD
print(ans)

Solution (Python):

import sys
MOD = 10**9 + 7
def main(N, M):
    if N > M:
        return 0
    fu = FactorialUtils(M)
    a = fu.permutation(M, N)
    t = Mint(a)
    for i in range(1, N + 1):
        t += (-1) ** (i & 1) * fu.choose(N, i) * fu.permutation(M - i, N - i)
    return t * a

class FactorialUtils:
    def __init__(self, n):
        self.fac = [1] * (n + 1)
        self.ifac = [1] * (n + 1)
        for i in range(2, n + 1): self.fac[i] = self.fac[i - 1] * i % MOD
        self.ifac[n] = pow(self.fac[n], MOD - 2, MOD)
        for i in range(n, 1, -1): self.ifac[i - 1] = self.ifac[i] * i % MOD

    def choose(self, n, r):
        if r < 0 or r > n: return 0
        return (self.fac[n] * self.ifac[n - r] % MOD) * self.ifac[r] % MOD

    def multichoose(self, u, k):
        return (self.fac[u + k - 1] * self.ifac[u - 1] % MOD) * self.ifac[k] % MOD

    def permutation(self, n, r):
        if r < 0 or r > n: return 0
        return self.fac[n] * self.ifac[n - r] % MOD

class Mint:
    def __init__(self, value=0):
        self.value = value % MOD
        if self.value < 0: self.value += MOD

    @staticmethod
    def get_value(x): return x.value if isinstance(x, Mint) else x

    def inverse(self):
        a, b = self.value, MOD
        u, v = 1, 0
        while b:
            t = a // b
            b, a = a - t * b, b
            v, u = u - t * v, v
        if u < 0: u += MOD
        return u

    def __repr__(self): return str(self.value)
    def __eq__(self, other): return self.value == other.value
    def __neg__(self): return Mint(-self.value)
    def __hash__(self): return hash(self.value)
    def __bool__(self): return self.value != 0

    def __iadd__(self, other):
        self.value = (self.value + Mint.get_value(other)) % MOD
        return self

    def __add__(self, other):
        new_obj = Mint(self.value)
        new_obj += other
        return new_obj
    __radd__ = __add__

    def __isub__(self, other):
        self.value = (self.value - Mint.get_value(other)) % MOD
        if self.value < 0: self.value += MOD
        return self

    def __sub__(self, other):
        new_obj = Mint(self.value)
        new_obj -= other
        return new_obj

    def __rsub__(self, other):
        new_obj = Mint(Mint.get_value(other))
        new_obj -= self
        return new_obj

    def __imul__(self, other):
        self.value = self.value * Mint.get_value(other) % MOD
        return self

    def __mul__(self, other):
        new_obj = Mint(self.value)
        new_obj *= other
        return new_obj
    __rmul__ = __mul__

    def __ifloordiv__(self, other):
        other = other if isinstance(other, Mint) else Mint(other)
        self *= other.inverse()
        return self

    def __floordiv__(self, other):
        new_obj = Mint(self.value)
        new_obj //= other
        return new_obj

    def __rfloordiv__(self, other):
        new_obj = Mint(Mint.get_value(other))
        new_obj //= self
        return new_obj

if __name__ == '__main__':
    input = sys.stdin.readline
    N, M = map(int, input().split())
    print(main(N, M))

Solution (Python):
import sys
from collections import deque
def input(): return sys.stdin.readline().strip()
mod = 10**9+7

class Combination:
    """
    O(n)の前計算を1回行うことで，O(1)でnCr mod mを求められる
    n_max = 10**6のとき前処理は約950ms (PyPyなら約340ms, 10**7で約1800ms)
    使用例：
    comb = Combination(1000000)
    print(comb(5, 3))  
    """
    def __init__(self, n_max, mod=10**9+7):
        self.mod = mod
        self.modinv = self.make_modinv_list(n_max)
        self.fac, self.facinv = self.make_factorial_list(n_max)
 
    def __call__(self, n, r):
        return self.fac[n] * self.facinv[r] % self.mod * self.facinv[n-r] % self.mod
 
    def make_factorial_list(self, n):
        
        
        fac = [1]
        facinv = [1]
        for i in range(1, n+1):
            fac.append(fac[i-1] * i % self.mod)
            facinv.append(facinv[i-1] * self.modinv[i] % self.mod)
        return fac, facinv
 
    def make_modinv_list(self, n):
        
        modinv = [0] * (n+1)
        modinv[1] = 1
        for i in range(2, n+1):
            modinv[i] = self.mod - self.mod//i * modinv[self.mod%i] % self.mod
        return modinv


comb = Combination(5 * 10**5)


def main():
    """
    「異なる」はネガティブだから逆を考えるべきなのか。。。
    そうすると包除原理が浮かんできそう
    A=[1,2,...,N]は仮定してよい。後で{}_M P_Nをかければよい。
    """
    N, M = map(int, input().split())
    ans = 0
    for i in range(N + 1):
        """
        B[j] = jなるものがi箇所ある場合。
        まずB[j] = jなるi箇所を選ぶためにcomb(N, i)必要。
        そのあとそれ以外の場所(N - i)箇所に埋める数を(M - i)個から選ぶための{}_{M-i} P_{N - i}が必要。
        """
        val = comb(N, i) * comb(M - i, N - i) * comb.fac[N - i] * (-1) ** i
        ans += val % mod
        
    ans *= comb(M, N) * comb.fac[N]
    print(ans % mod)
    

if __name__ == "__main__":
    main()

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

class Facts():
   def __init__(self, max_num=10**5, p=10**9+7):
      self.p = p
      self.max_num = max_num
      self.fact = [1] * (self.max_num + 1)
      self.rev = [1] * (self.max_num + 1)
      for i in range(1, self.max_num + 1):
         self.fact[i] = (self.fact[i-1] * i) % self.p
      self.rev[self.max_num] = self.mod_pow(self.fact[self.max_num], self.p-2)
      for i in range(self.max_num-1, 0, -1):
         self.rev[i] = self.rev[i+1] * (i + 1) % self.p

   def get_fact(self, n):
      return self.fact[n]

   def comb(self, n, k):
      if n < 0 or k < 0 or n < k:
         return 0
      if n == 0 or k == 0:
         return 1
      res = ((self.fact[n] * self.rev[k] % self.p) * self.rev[n - k]) % self.p
      return res

   def perm(self, n, k):
      if n < 0 or k < 0 or n < k:
         return 0
      if n == 0 or k == 0:
         return 1
      a = self.fact[n]
      b = self.fact[n-k]
      return (a * self.mod_pow(b, self.p-2)) % self.p

   def mod_pow(self, a, b):
      ans = 1
      while b > 0:
         if b & 1:
            ans = ans * a % self.p
         a = a * a % self.p; b >>= 1
      return ans

facts = Facts(max_num=5*10**5)

ans = 0
for k in range(n+1):
   tmp = (-1) ** k 
   tmp *= facts.comb(n, k) * facts.perm(m, k)
   tmp %= MOD
   tmp *= (facts.perm(m-k, n-k)) ** 2
   tmp %= MOD
   ans += tmp
   ans %= MOD
print(ans)

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

mod = 10 ** 9 + 7

fac = [1, 1]
finv = [1, 1]
inv = [0, 1]

def comb(n, r):
  return fac[n] * ( finv[r] * finv[n-r] % mod ) % mod

def perm(n, r):
  return  fac[n] * finv[n-r] % mod
 
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 )

ans = 0

for K in range(N+1):
  ans += comb(N, K) * (-1)**(K%2) % mod * perm(M, K) % mod * pow(perm(M-K, N-K), 2, mod)
  ans %= mod

print(ans)

Solution (Python):
def modInverse(a, p):  
    assert a % p != 0
    return pow(a, p - 2, p)


MOD = 10 ** 9 + 7
N, M = [int(x) for x in input().split()]


MAX_N = M

fact = [1]
for i in range(1, MAX_N + 1):
    fact.append((fact[-1] * i) % MOD)


invFact = [0] * (MAX_N + 1)
invFact[MAX_N] = modInverse(fact[MAX_N], MOD)
for i in range(MAX_N - 1, -1, -1):
    invFact[i] = (invFact[i + 1] * (i + 1)) % MOD
    assert fact[i] * invFact[i] % MOD == 1


def nCr(n, r):  
    
    return (fact[n] * invFact[r] * invFact[n - r]) % MOD



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

print(ans)

Solution (Python):
import sys

sys.setrecursionlimit(10 ** 7)
rl = sys.stdin.readline


class Combination:
    def __init__(self, n: int, mod: int):
        self.mod = mod
        self.fact = [0] * (n + 1)
        self.factinv = [0] * (n + 1)
        self.inv = [0] * (n + 1)
        
        self.fact[0] = self.fact[1] = 1
        self.factinv[0] = self.factinv[1] = 1
        self.inv[1] = 1
        for i in range(2, n + 1):
            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 ncr(self, n: int, r: int):
        if r < 0 or n < r:
            return 0
        r = min(r, n - r)
        return self.fact[n] * self.factinv[r] % self.mod * self.factinv[n - r] % self.mod
    
    def nhr(self, n: int, r: int):
        return self.ncr(n + r - 1, r)
    
    def npr(self, n: int, r: int):
        if r < 0 or n < r:
            return 0
        return self.fact[n] * self.factinv[n - r] % self.mod


def solve():
    N, M = map(int, rl().split())
    MOD = 10 ** 9 + 7
    
    com = Combination(5 * 10 ** 5, MOD)
    ans = 0
    for i in range(N + 1):
        exp = -1 if i % 2 == 1 else 1
        ans += com.ncr(N, i) * exp * com.npr(M, i) * com.npr(M - i, N - i) ** 2
        ans %= MOD
    print(ans)


if __name__ == '__main__':
    solve()

Solution (Python):
N,M=map(int,input().split())
mod=10**9+7
Fact=[0 for i in range(M+1)]
Finv=[0 for i in range(M+1)]
Fact[0]=1
def inv(x):
    return pow(x,mod-2,mod)
for i in range(M):
    Fact[i+1]=Fact[i]*(i+1)
    Fact[i+1]%=mod
Finv[M]=inv(Fact[M])
for i in range(M-1,-1,-1):
    Finv[i]=(i+1)*Finv[i+1]
    Finv[i]%=mod
def P(n,k):
    return (Fact[n]*Finv[n-k])%mod
def C(n,k):
    return (Fact[n]*Finv[n-k]*Finv[k])%mod
ans=0
for k in range(N+1):
    if k%2==0:
        ans+=C(N,k)*P(M-k,N-k)
    else:
        ans-=C(N,k)*P(M-k,N-k)
    ans%=mod
print((ans*P(M,N))%mod)

Solution (Python):
n, m = map(int, input().split())
mod = pow(10, 9) + 7
x = 0
y = m - n + 1
u, v = n, 1
c, p = 1, 1
for i in range(n, 0, -1):
    if i % 2 == 1:
        x += ((c * p) % mod)
    else:
        x -= ((c * p) % mod)
    x %= mod
    p *= y
    p %= mod
    y += 1
    c = (c * u * pow(v, mod - 2, mod)) % mod
    u -= 1
    v += 1
p = 1
for i in range(n):
    p *= (m - i)
    p %= mod
ans = (pow(p, 2, mod) - (p * x)) % mod
print(ans)

Solution (Python):

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


def table():
    k = 5 * 10**5 + 5  
    fac = [1] * k
    finv = [1] * k
    inv = [0] * k
    inv[1] = 1
    for i in range(2, k):
        fac[i] = fac[i - 1] * i % mod
        inv[i] = mod - inv[mod % i] * (mod // i) % mod
        finv[i] = finv[i - 1] * inv[i] % mod
    return fac, finv


fac, finv = table()


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


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


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