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):

from collections import defaultdict, deque, Counter, OrderedDict
from bisect import bisect_left, bisect_right
from functools import reduce, lru_cache
from heapq import heappush, heappop, heapify

import itertools, bisect
import math, fractions
import sys, copy

def L(): return sys.stdin.readline().split()
def I(): return int(sys.stdin.readline().rstrip())
def S(): return list(sys.stdin.readline().rstrip())
def LI(): return [int(x) for x in sys.stdin.readline().split()]
def LI1(): return [int(x) - 1 for x in sys.stdin.readline().split()]
def LS(): return [list(x) for x in sys.stdin.readline().split()]
def IR(n): return [I() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def LIR1(n): return [LI1() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
def LR(n): return [L() for _ in range(n)]

alphabets = "abcdefghijklmnopqrstuvwxyz"
ALPHABETS = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"

sys.setrecursionlimit(1000000)
dire = [[1, 0], [0, 1], [-1, 0], [0, -1]]
dire8 = [[1, 0], [1, 1], [0, 1], [-1, 1], [-1, 0], [-1, -1], [0, -1], [1, -1]]
MOD = 1000000007


class FactInv:
    def __init__(self, N, MOD=1000000007):
        fact = [1]*(N+1)
        for i in range(1, N+1):
            fact[i] = fact[i - 1] * i % MOD

        inv = [None]*(N+1)
        inv[N] = pow(fact[N], MOD - 2, MOD)

        for i in range(N)[::-1]:
            inv[i] = inv[i + 1] * (i + 1) % MOD

        self.N = N
        self.MOD = MOD
        self.fact = fact
        self.inv = inv

    def perm(self, a, b):
        if a > self.N or b > self.N:
            raise ValueError(
                "\nPermutaion arguments are bigger than N\n N = {}, a = {}, b = {}".format(self.N, a, b))
        return self.fact[a] * self.inv[a-b] % self.MOD

    def comb(self, a, b):
        if a > self.N or b > self.N:
            raise ValueError(
                "\nCombination arguments are bigger than N\n N = {}, a = {}, b = {}".format(self.N, a, b))
        return self.fact[a] * self.inv[b] * self.inv[a-b] % self.MOD


def main():
    N, M = LI()
    factinv = FactInv(M)

    a = factinv.perm(M, N)

    tmp = a
    for k in range(1, N + 1):
        tmp = (tmp + factinv.comb(N, k) * factinv.perm(M - k, N - k) * pow(-1, k)) % MOD

    print(a * tmp % MOD)





if __name__ == '__main__':
    main()

Solution (Python):
from sys import stdin
import sys
import math
from functools import reduce
import functools
import itertools
from collections import deque,Counter
from operator import mul
from functools import reduce


N = 10**6
mod = 10**9+7

inv_t = [0] + [1]
for i in range(2, N):
    inv_t += [inv_t[mod % i] * (mod - int(mod / i)) % mod]

kai = [1, 1]
rev_kai = [1, inv_t[1]]
for i in range(2, N):
    kai.append(kai[-1] * i % mod)
    rev_kai.append(rev_kai[-1] * inv_t[i] % mod)

def cmb(n, r):
    return kai[n] * rev_kai[r] * rev_kai[n - r] % mod

def prm(n, r):
    return kai[n] * rev_kai[n - r] % mod

n,m = list(map(int, input().split()))

a = 0
for k in range(n+1):
    a += prm(m, k)*prm(m-k,n-k)**2*(-1)**k*cmb(n,k) % mod

print(a % mod)

Solution (Python):
import sys
INF = 1 << 60
MOD = 10**9 + 7 
sys.setrecursionlimit(2147483647)
input = lambda:sys.stdin.readline().rstrip()

class modfact(object):
    def __init__(self, n):
        fact, invfact = [1] * (n + 1), [1] * (n + 1)
        for i in range(1, n + 1): fact[i] = i * fact[i - 1] % MOD
        invfact[n] = pow(fact[n], MOD - 2, MOD)
        for i in range(n - 1, -1, -1): invfact[i] = invfact[i + 1] * (i + 1) % MOD
        self._fact, self._invfact = fact, invfact

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

    def invfact(self, n):
        return self._invfact[n]

    def comb(self, n, k):
        return self._fact[n] * self._invfact[k] % MOD * self._invfact[n - k] % MOD if 0 <= k <= n else 0

    def perm(self, n, k):
        return self._fact[n] * self._invfact[n - k] % MOD if 0 <= k <= n else 0

def resolve():
    n, m = map(int, input().split())
    mf = modfact(m)
    ans = 0
    for k in range(n + 1):
        if k & 1:
            ans -= mf.comb(n, k) * mf.perm(m - k, n - k)
        else:
            ans += mf.comb(n, k) * mf.perm(m - k, n - k)
        ans %= MOD
    ans *= mf.perm(m, n)
    ans %= MOD
    print(ans)
resolve()

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

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

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

def cmb(n, r, mod=MOD):
    return (pmu(n, r) * pow(FCT[r], mod-2, mod)) % mod

def solve():
    
    
    ans = 0
    for i in range(N+1):
        cp = (cmb(N, i) * pmu(M-i, N-i))%MOD
        ans = (ans + (1 if i%2==0 else -1) * cp) % MOD
    return (ans * pmu(M, N)) % MOD

if __name__ == "__main__":
    print(solve())

Solution (Python):
import sys
import itertools

import time
import math
import heapq
from collections import defaultdict
sys.setrecursionlimit(10 ** 7)
 
INF = 10 ** 18
MOD = 10 ** 9 + 7
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines


MAX = 10 ** 6
MOD = 10 ** 9 + 7

fac = [0 for i in range(MAX)]
finv = [0 for i in range(MAX)]
inv = [0 for i in range(MAX)]

def comInit(mod):
    fac[0], fac[1] = 1, 1
    finv[0], finv[1] = 1, 1
    inv[1] = 1
    for i in range(2, MAX):
        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 com(n, r, mod):
    if n < r:
        return 0
    if n < 0 or r < 0:
        return 0
    return fac[n] * (finv[r] * finv[n - r] % mod) % mod

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

comInit(MOD)

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

ans = 0
for i in range(N + 1):
    now = com(N, i, MOD)
    now *= p(M - i, N - i, MOD)
    if i & 1 > 0:
        now = -now
    ans += now

ans *= p(M, N, MOD)
print(ans % MOD)

Solution (Python):
n, m = map(int, input().split())

mod = 10 ** 9 + 7
N = 10 ** 6


inv_t = [0]+[1]
for i in range(2, N):
    inv_t += [inv_t[mod % i] * (mod - int(mod / i)) % mod]


kai = [1, 1]
rev_kai = [1, inv_t[1]]
for i in range(2, N):
    kai.append(kai[-1] * i % mod)
    rev_kai.append(rev_kai[-1] * inv_t[i] % mod)


def cmb(n, r):
    return kai[n] * rev_kai[r] * rev_kai[n-r] % mod

def prm(n, r):
    return kai[n] * rev_kai[n-r] % mod

ans = 0
tmp = 1


for k in range(n+1):
    
    
    ans += tmp * cmb(n, k) * prm(m-k, n-k) % mod 
    ans %= mod
    tmp *= -1


ans *= prm(m, n)
ans %= mod

print(ans)

Solution (Python):
import sys
N,M=map(int, sys.stdin.readline().split())

mod=10**9+7
fact, inv_fact=[1], [1]
fact_tmp=1
for i in range(1, M+1):
    fact_tmp*=i
    fact_tmp%=mod
    fact.append(fact_tmp)
    inv_fact.append(pow(fact_tmp, mod-2, mod))

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

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

ans=0
for k in range(N+1):
    if k%2==0: ans+=ncr(N,k)*ncr(M,k)*npr(k,k)*(npr(M-k,N-k)**2)
    else:      ans-=ncr(N,k)*ncr(M,k)*npr(k,k)*(npr(M-k,N-k)**2)
    ans%=mod
print ans

Solution (Python):

import sys
N,M=map(int, sys.stdin.readline().split())

mod=10**9+7
fact, inv_fact=[1], [1]
fact_tmp=1
for i in range(1, M+1):
    fact_tmp*=i
    fact_tmp%=mod
    fact.append(fact_tmp)
    inv_fact.append(pow(fact_tmp, mod-2, mod))

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

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

ans=0
for k in range(N+1):
    if k%2==0: ans+=ncr(N,k)*ncr(M,k)*npr(k,k)*(npr(M-k,N-k)**2)
    else:      ans-=ncr(N,k)*ncr(M,k)*npr(k,k)*(npr(M-k,N-k)**2)
    ans%=mod
print ans

Solution (Python):
n, m = map(int, input().split())
mod = 10**9 + 7
def setfact(l, mod):
    global fact
    global ifact
    fact = [1]
    for i in range(1, l+1):
        fact.append(fact[-1]*i%mod)
    ifact = [0] * (l+1)
    ifact[-1] = pow(fact[-1], mod-2, mod)
    for i in range(l, 0, -1):
        ifact[i-1] = ifact[i] * i %mod
def comb(n, m, mod):
    global fact
    global ifact
    return fact[n]*ifact[n-m]*ifact[m]%mod
def p(n, m, mod):
    global fact
    global ifact
    return fact[n]*ifact[n-m]%mod
ans = 0
setfact(m, mod)
for k in range(n+1):
    ans += p(m-k, n-k, mod)*comb(n, k, mod) * (-1)**k
    ans %= mod
    if ans < 0:
        ans += mod
ans *= p(m, n, mod)
ans %= mod
print(ans)

Solution (Python):
n,m=map(int,input().split())
mod=10**9+7
f=[1]
for i in range(1,max(m,n)+1):f.append(f[i-1]*i%mod)
def comb(n,r):return f[n]*pow(f[r],mod-2,mod)*pow(f[n-r],mod-2,mod)%mod
def pern(n,r):return f[n]*pow(f[n-r],mod-2,mod)%mod
a=[0]*(n+1)
a[n]=1
for i in range(n-1,-1,-1):
    a[i]=(comb(n,n-i)*pern(m-i,n-i)-a[i+1])%mod
print(a[0]*pern(m,n)%mod)

Solution (Python):
class Factorial():
    def __init__(self, n, mod):
        self.mod = mod
        self.fct = [0 for _ in range(n + 1)]
        self.inv = [0 for _ in range(n + 1)]
        self.fct[0] = 1
        self.inv[0] = 1
        for i in range(n):
            self.fct[i + 1] = self.fct[i] * (i + 1) % mod
        self.inv[n] = pow(self.fct[n], mod - 2, mod)
        for i in range(n)[::-1]:
            self.inv[i] = self.inv[i + 1] * (i + 1) % mod

    def fact(self, m):
        return self.fct[m]

    def invf(self, m):
        return self.inv[m]

    def perm(self, m, k):
        if m < k: return 0
        return self.fct[m] * self.inv[m - k] % self.mod

    def comb(self, m, k):
        if m < k: return 0
        return self.fct[m] * self.inv[k] * self.inv[m - k] % self.mod

MOD = 1000000007

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

f = Factorial(M, MOD)

res = f.perm(M, N)

for k in range(1, N + 1):
    res -= (-1)**(k - 1) * f.comb(N, k) * f.perm(M - k, N - k)

res *= f.perm(M, N)
res %= MOD

print(res)

Solution (Python):

import sys
sys.setrecursionlimit(10**9)
INF=10**18
MOD=10**9+7
input=lambda: sys.stdin.readline().rstrip()
YesNo=lambda b: bool([print('Yes')] if b else print('No'))
YESNO=lambda b: bool([print('YES')] if b else print('NO'))
int1=lambda x:int(x)-1

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

def COMinit(n,MOD):
    fac,finv,inv=[0]*max(2,n+1),[0]*max(2,n+1),[0]*max(2,n+1)
    fac[0]=fac[1]=1
    finv[0]=finv[1]=1
    inv[1]=1
    for i in range(2,(n+1)):
        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,inv

fac,finv,inv=COMinit(M,MOD)

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

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

ans=0
k=1

for i in range(N+1):
    ans+=k*COM(N,i)*PER(M,i)*(PER(M-i,N-i)**2)
    ans%=MOD
    k*=-1
print(ans)

Solution (Python):



MOD = 10**9 + 7

def resolve():
    N,M = map(int,input().split())
    a = 1
    d = [1]
    for i in range(N):
        d.append(((M-N+i)*d[i]+i*d[i-1])%MOD)
        a=a*(M-i)%MOD
    print(a*d[-1]%MOD)
resolve()

Solution (Python):
'''
Let g(x) = 序列內數字唯一，且最多 x 個位置不同。
Let f(x) = 序列內數字唯一，且給定 x 個位置時，使這些位置不同的方法數。

題目所求是 C(N, N) f(N)。

其中 g(x) 代表有至少 N - x 個位置是相同。可以分三個部份填數字，注意順序會影響，所以要乘上全排列：
(1) N - x 個位置（數字必需相同）：C(M, N - x) * (N - x)! = P(M, N - x)
(2) 序列 A 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
(3) 序列 B 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
g(x) = (1) * (2) * (3)

g(x) 是「選定 0 個位置，使這些位置不同的方法數」、「選定 1 個位置，使這些位置不同的方法數」、…、「選定 x 個位置，使這些位置不同的方法數」，即
g(x) = C(x, 0) f(0) + C(x, 1) f(1) + ... + C(x, x) f(x)
     = sum( C(x, i) * f(i) for i in [0, x] )
根據二項式反演，得到
f(x) = sum( (-1) ** (x - i) * C(x, i) * g(i) for i in [0, x] )
'''

class CombMod:
    def __init__(self, V, p):
        self.fact = [1] * V
        self.finv = [1] * V
        for i in range(1, V):
            self.fact[i] = self.fact[i - 1] * i % p
        self.finv[-1] = pow(self.fact[-1], p - 2, p)
        for i in range(V - 2, 0, -1):
            self.finv[i] = self.finv[i + 1] * (i + 1) % p
        self.p = p
        
    def fact(self, a):
        return self.fact[a]
    
    def finv(self, a):
        return self.finv[a]

    def comb(self, a, b):
        return self.fact[a] * self.finv[b] % self.p * self.finv[a - b] % self.p
        
    def perm(self, a, b):
        return self.fact[a] * self.finv[a - b] % self.p
    
    def hcomb(self, a, b):
        return self.comb(a + b - 1, b)


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

f_N = 0
for i in range(0, N + 1):
    p1 = tool.perm(M, N - i)
    p2 = tool.perm(M - (N - i), i)
    g_i = p1 * p2 % mod * p2 % mod
    val = (-1)**((N - i) % 2) * tool.comb(N, i) % mod * g_i % mod
    f_N = (f_N + val) % mod
ans = tool.comb(N, N) * f_N % mod
print(ans)

Solution (Python):
'''
Let g(x) = 序列內數字唯一，且最多 x 個位置不同。
Let f(x) = 序列內數字唯一，且給定 x 個位置時，使這些位置不同的方法數。

題目所求是 C(N, N) f(N)。

其中 g(x) 代表有至少 N - x 個位置是相同。可以分三個部份填數字，注意順序會影響，所以要乘上全排序：
(1) N - x 個位置（數字必需相同）：C(M, N - x) * (N - x)! = P(M, N - x)
(2) 序列 A 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
(3) 序列 B 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
g(x) = (1) * (2) * (3)

g(x) 是「選定 0 個位置，使這些位置不同的方法數」、「選定 1 個位置，使這些位置不同的方法數」、…、「選定 x 個位置，使這些位置不同的方法數」，即
g(x) = C(x, 0) f(0) + C(x, 1) f(1) + ... + C(x, x) f(x)
     = sum( C(x, i) * f(i) for i in [0, x] )
根據二項式反演，得到
f(x) = sum( (-1) ** (x - i) * C(x, i) * g(i) for i in [0, x] )
'''

class CombMod:
    def __init__(self, V, p):
        self.fact = [1] * V
        self.finv = [1] * V
        for i in range(1, V):
            self.fact[i] = self.fact[i - 1] * i % p
        self.finv[-1] = pow(self.fact[-1], p - 2, p)
        for i in range(V - 2, 0, -1):
            self.finv[i] = self.finv[i + 1] * (i + 1) % p
        self.p = p
        
    def fact(self, a):
        return self.fact[a]
    
    def finv(self, a):
        return self.finv[a]

    def comb(self, a, b):
        return self.fact[a] * self.finv[b] % self.p * self.finv[a - b] % self.p
        
    def perm(self, a, b):
        return self.fact[a] * self.finv[a - b] % self.p
    
    def hcomb(self, a, b):
        return self.comb(a + b - 1, b)


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

f_N = 0
for i in range(0, N + 1):
    p1 = tool.perm(M, N - i)
    p2 = tool.perm(M - (N - i), i)
    g_i = p1 * p2 % mod * p2 % mod
    val = (-1)**((N - i) % 2) * tool.comb(N, i) % mod * g_i % mod
    f_N = (f_N + val) % mod
ans = tool.comb(N, N) * f_N % mod
print(ans)

Solution (Python):
'''
Let g(x) = 序列內數字唯一，且最多 x 個位置不同。
Let f(x) = 序列內數字唯一，且給定 x 個位置時，使這些位置不同的方法數。

題目所求是 C(N, N) f(N)。

其中 g(x) 代表有至少 N - x 個位置是相同。可以分三個部份填數字，注意順序會影響，所以要乘上全排序：
(1) N - x 個位置（數字必需相同）：C(M, N - x) * (N - x)! = P(M, N - x)
(2) 序列 A 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
(3) 序列 B 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
g(x) = (1) * (2) * (3)

g(x) 是「選定 0 個位置，使這些位置不同的方法數」、「選定 1 個位置，使這些位置不同的方法數」、…、「選定 N 個位置，使這些位置不同的方法數」，即

g(x) = C(x, 0) f(0) + C(x, 1) f(1) + ... + C(x, x) f(x)
     = sum( C(x, i) * f(i) for i in [0, x] )
根據二項式反演，得到
f(x) = sum( (-1) ** (x - i) * C(x, i) * g(i) for i in [0, x] )
'''

class CombMod:
    def __init__(self, V, p):
        self.fact = [1] * V
        self.finv = [1] * V
        for i in range(1, V):
            self.fact[i] = self.fact[i - 1] * i % p
        self.finv[-1] = pow(self.fact[-1], p - 2, p)
        for i in range(V - 2, 0, -1):
            self.finv[i] = self.finv[i + 1] * (i + 1) % p
        self.p = p
        
    def fact(self, a):
        return self.fact[a]
    
    def finv(self, a):
        return self.finv[a]

    def comb(self, a, b):
        return self.fact[a] * self.finv[b] % self.p * self.finv[a - b] % self.p
        
    def perm(self, a, b):
        return self.fact[a] * self.finv[a - b] % self.p
    
    def hcomb(self, a, b):
        return self.comb(a + b - 1, b)


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

f_N = 0
for i in range(0, N + 1):
    p1 = tool.perm(M, N - i)
    p2 = tool.perm(M - (N - i), i)
    g_i = p1 * p2 % mod * p2 % mod
    val = (-1)**((N - i) % 2) * tool.comb(N, i) % mod * g_i % mod
    f_N = (f_N + val) % mod
ans = tool.comb(N, N) * f_N % mod
print(ans)

Solution (Python):
'''
Let g(x) = 序列內數字唯一，且最多 x 個位置不同。
Let f(x) = 序列內數字唯一，且給定 x 個位置時，使這些位置不同的方法數。

題目所求是 C(N, N) f(N)。

其中 g(x) 代表有至少 N - x 個位置是相同。可以分三個部份填數字，注意順序會影響，所以要乘上全排序：
(1) N - x 個位置（數字必需相同）：C(M, N - x) * (N - x)! = P(M, N - x)
(2) 序列 A 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
(3) 序列 B 其餘位置（可以隨便填，但不能用到 (1) 的數字）：C(M - (N - x), x) * x! = P(M - (N - x), x)
g(x) = (1) * (2) * (3)

g(x) 是「選定 0 個位置，使這些位置不同的方法數」、「選定 1 個位置，使這些位置不同的方法數」、…、「選定 N 個位置，使這些位置不同的方法數」，即

g(x) = C(x, 0) f(0) + C(x, 1) f(1) + ... + C(x, x) f(x)
     = sum( C(x, i) * f(i) for i in [0, x] )
根據二項式反演，得到
f(x) = sum( (-1) ** (x - i) * C(x, i) * g(i) for i in [0, x] )
'''

class CombMod:
    def __init__(self, V, p):
        self.fact = [1] * V
        self.finv = [1] * V
        for i in range(1, V):
            self.fact[i] = self.fact[i - 1] * i % p
        self.finv[-1] = pow(self.fact[-1], p - 2, p)
        for i in range(V - 2, 0, -1):
            self.finv[i] = self.finv[i + 1] * (i + 1) % p
        self.p = p
        
    def fact(self, a):
        return self.fact[a]
    
    def finv(self, a):
        return self.finv[a]

    def comb(self, a, b):
        return self.fact[a] * self.finv[b] % self.p * self.finv[a - b] % self.p
        
    def perm(self, a, b):
        return self.fact[a] * self.finv[a - b] % self.p
    
    def hcomb(self, a, b):
        return self.comb(a + b - 1, b)


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

ans = 0
for i in range(0, N + 1):
    p1 = tool.perm(M, N - i)
    p2 = tool.perm(M - (N - i), i)
    g_i = p1 * p2 % mod * p2 % mod
    val = pow(-1, N - i, mod) * tool.comb(N, i) % mod * g_i % mod
    ans = (ans + val) % mod
print(ans)

Solution (Python):
from math import factorial
mod = 10**9+7 
N = 10**6
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 )

def nCr(n, r):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

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


N,M=map(int,input().split())
ans=nPr(M,N)**2%mod
for i in range(1,N+1):
    p=nPr(M-i,N-i)**2%mod
    p=p*nPr(N,i)%mod
    p=p*nCr(M,i)%mod
    ans-=((-1) ** ((i % 2) ^ 1))*p
    ans%=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

def per(n, r, mod):
    if (r < 0 or r > n):
        return 0

    return g1[n] * g2[n-r]  % mod

mod = 10 ** 9 + 7
N = 10 ** 6 
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=per(m,n,mod)


cnt=0
for i in range(n+1):
    if i%2==0:
        cnt+=per(m-i,n-i,mod)*cmb(n,i,mod)
    else:
        cnt-=per(m-i,n-i,mod)*cmb(n,i,mod)

    cnt%=mod
print((ans*cnt)%mod)

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

def cmb(n,r):
  global mod
  if r<0 or r>n:
    return 0
  return (g1[n]*g2[r]*g2[n-r])%mod

def perm(n,r):
  global mod
  if r>n:
    return g1[n]
  if r<0:
    return 0
  return (g1[n]*g2[n-r])%mod

g1=[1,1]
g2=[1,1]
inv=[0,1]
for i in range(2,M+M+5):
  g1.append((g1[-1]*i)%mod)
  inv.append((-inv[mod%i]*(mod//i))%mod)
  g2.append((g2[-1]*inv[-1])%mod)
P=0
for i in range(N+1):
  P+=pow(-1,i)*cmb(N,i)*perm(M,i)*perm(M-i,N-i)**2
print(P%mod)

Solution (Python):
class Modulo_Error(Exception):
    pass

class Modulo():
    def __init__(self,a,n):
        self.a=a%n
        self.n=n

    def __str__(self):
        return "{} (mod {})".format(self.a,self.n)

    
    def __pos__(self):
        return self

    def __neg__(self):
        return  Modulo(-self.a,self.n)

    
    def __eq__(self,other):
        if isinstance(other,Modulo):
            return (self.a==other.a) and (self.n==other.n)
        elif isinstance(other,int):
            return (self-other).a==0

    def __neq__(self,other):
        return not(self==other)
    
    
    def __add__(self,other):
        if isinstance(other,Modulo):
            if self.n!=other.n:
                raise Modulo_Error("異なる法同士の演算です.")
            return Modulo(self.a+other.a,self.n)
        elif isinstance(other,int):
            return Modulo(self.a+other,self.n)

    def __radd__(self,other):
        if isinstance(other,int):
            return Modulo(self.a+other,self.n)
        
    
    def __sub__(self,other):
        return self+(-other)

    def __rsub__(self,other):
        if isinstance(other,int):
            return -self+other
        
    
    def __mul__(self,other):
        if isinstance(other,Modulo):
            if self.n!=other.n:
                raise Modulo_Error("異なる法同士の演算です.")
            return Modulo(self.a*other.a,self.n)
        elif isinstance(other,int):
            return Modulo(self.a*other,self.n)
        
    def __rmul__(self,other):
        if isinstance(other,int):
            return Modulo(self.a*other,self.n)
        
    
    def Modulo_Inverse(self):
        x0, y0, x1, y1 = 1, 0, 0, 1
        a,b=self.a,self.n
        while b != 0:
            q, a, b = a // b, b, a % b
            x0, x1 = x1, x0 - q * x1
            y0, y1 = y1, y0 - q * y1

        if a!=1:
            raise Modulo_Error("{}の逆数が存在しません".format(self))
        else:
            return Modulo(x0,self.n)
        
    
    def __truediv__(self,other):
        return self*(other.Modulo_Inverse())
    
    def __rtruediv__(self,other):
        return other*(self.Modulo_Inverse())

    
    def __pow__(self,m):
        u=abs(m)

        r=Modulo(1,self.n)

        while u>0:
            if u%2==1:
                r*=self
            self*=self
            u=u>>1

        if m>=0:
            return r
        else:
            return r.Modulo_Inverse()


def Factorial_Modulo_List(N,M,Inverse=False):
    X=[0]*(N+1)
    X[0]=Modulo(1,M)
    
    for i in range(1,N+1):
        X[i]=X[i-1]*i

    if Inverse:
        Y=[0]*(N+1)
        Y[-1]=1/X[-1]
        
        for j in range(N-1,-1,-1):
            Y[j]=Y[j+1]*(j+1)
        return X,Y
    else:
        return X

def nPr(n,r):
    return F[n]*G[n-r]

def nCr(n,r):
    return F[n]*G[n-r]*G[r]

A=10**9+7
N,M=map(int,input().split())
F,G=Factorial_Modulo_List(M,A,True)

Q=0
sgn=1
for k in range(N+1):
    Q+=sgn*nCr(N,k)*nPr(M,k)*nPr(M-k,N-k)**2
    sgn*=-1
print(Q.a)

Solution (Python):
n, m = map(int, input().split(' '))
mod = int(1e9 + 7)

l = [1] * (m+1)
for i in range(1, m+1):
	l[i] = (l[i-1] * i) % mod

l2 = [1] * (m+1)
l2[m] = pow(l[m], mod-2, mod)
for i in range(m-1, 0, -1):
	l2[i] = (l2[i+1] * (i+1)) % mod

s = 0
for k in range(n+1):
	n_choose_k = (l[n] * l2[k] * l2[n-k]) % mod
	m_perm_k = (l[m] * l2[m-k]) % mod
	mk_perm_nk = (l[m-k] * l2[m-n]) % mod
	s += ((-1) ** k * n_choose_k * m_perm_k * mk_perm_nk ** 2 ) % mod 

print(s % mod)

Solution (Python):

def main():
    import sys
    imput = sys.stdin.readline
    if sys.argv[-1] == 'ONLINE_JUDGE':
        import numba
        from numba.pycc import CC
        i8 = numba.int64
        cc = CC('my_module')
 
        def cc_export(f, signature):
            cc.export(f.__name__, signature)(f)
            return numba.njit(f)
 
            fact_table = cc_export(fact_table, (i8, i8))
            main = cc_export(main, (i8, i8, i8))
            cc.compile()

    
    
    
    n, m = map(int, input().split())

    
    mod = pow(10, 9) + 7

    n_ = 5 * pow(10, 5) + 5
    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 cmb(n,r):
        if n < 0 or r < 0 or r > n: return 0
        return fun[n] * rev[r] % mod * rev[n-r] % mod
    def perm(n, r):
        if n < 0 or r < 0 or r > n: return 0
        return fun[n] * rev[n-r] % mod

    
    
    

    import math
    answer = 0
    for i in range(n+1):
        temp = perm(n, i) * cmb(m, i) * pow(perm(m-i, n-i), 2)
        temp %= mod
        if i % 2 == 0:
            answer += temp
        else:
            answer -= temp
    print(answer % mod)

    
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()
    ans = 0

    cmb = Combination(M)

    sign = +1
    for i in range(N+1):
        ans += sign * cmb.P(M-i, N-i) * cmb(N, i)
        ans %= MOD
        sign *= -1
    ans *= cmb.P(M, N)
    ans %= MOD
    print(ans)
    return


if __name__ == '__main__':
    main()

Solution (Python):

def main():
    import sys
    imput = sys.stdin.readline
    
    
    
    n, m = map(int, input().split())

    
    mod = pow(10, 9) + 7

    n_ = 5 * pow(10, 5) + 5
    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 cmb(n,r):
        if n < 0 or r < 0 or r > n: return 0
        return fun[n] * rev[r] % mod * rev[n-r] % mod
    def perm(n, r):
        if n < 0 or r < 0 or r > n: return 0
        return fun[n] * rev[n-r] % mod

    
    
    

    import math
    answer = 0
    for i in range(n+1):
        temp = perm(n, i) * cmb(m, i) * pow(perm(m-i, n-i), 2)
        temp %= mod
        if i % 2 == 0:
            answer += temp
        else:
            answer -= temp
    print(answer % mod)

    
if __name__ == "__main__":
    main()

Solution (Python):

def main():
    import sys
    imput = sys.stdin.readline
    
    
    
    n, m = map(int, input().split())

    
    mod = pow(10, 9) + 7

    n_ = 5 * pow(10, 5) + 5
    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 cmb(n,r):
        if n < 0 or r < 0 or r > n: return 0
        return fun[n] * rev[r] % mod * rev[n-r] % mod
    def perm(n, r):
        if n < 0 or r < 0 or r > n: return 0
        return fun[n] * rev[n-r] % mod

    
    
    

    import math
    answer = 0
    for i in range(n+1):
        temp = perm(n, i) * cmb(m, i) * pow(perm(m-i, n-i), 2)
        temp %= mod
        if i % 2 == 0:
            answer += temp
        else:
            answer -= temp
    print(answer % mod)

    
if __name__ == "__main__":
    main()