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):
N,M = map(int,input().split())
MOD = 10**9+7

MAXN = N+5
fac = [1,1] + [0]*MAXN
finv = [1,1] + [0]*MAXN
inv = [0,1] + [0]*MAXN
for i in range(2,MAXN+2):
    fac[i] = fac[i-1] * i % MOD
    inv[i] = -inv[MOD%i] * (MOD // i) % MOD
    finv[i] = finv[i-1] * inv[i] % MOD
def comb(n,r):
    if n < r: return 0
    if n < 0 or r < 0: return 0
    return fac[n] * (finv[r] * finv[n-r] % MOD) % MOD

ans = 1
for i in range(N):
    ans *= (M-i)
    ans %= MOD

ps = [1]
for i in range(N):
    ps.append((ps[-1] * (M-N+1+i)) % MOD)
ps.reverse()

tmp = 0
for i,p in enumerate(ps):
    t = p * comb(N,i)
    if i%2:
        t *= -1
    tmp += t

ans *= tmp
ans %= MOD
print(ans)

Solution (Python):
size = 500001
mod = 10 ** 9 + 7
fac = [0] * size
finv = [0] * size
inv = [0] * size

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

g1 = [1, 1]
g2 = [1, 1]
inverse = [0, 1]

for i in range( 2, size + 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, mod):
    if ( r<0 or r>n ):
        return 0
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod

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

ans = 0
flag = 1

for k in range(N+1):
    ans += flag * nPr(M - k, N - k, mod) * nCr(N, k, mod) * nPr(M, N, mod)
    flag *= -1
    ans %= mod

print(ans)

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

def I(): return int(input())
def MI(): return map(int, input().split())
def LI(): return list(map(int, input().split()))

def main():
    mod=10**9+7
    N,M=MI()
    
    
    
    def cmb(n, r, mod):
        if (r < 0) or (n < r):
            return 0
        r = min(r, n - r)
        return (fact[n] * factinv[r] * factinv[n-r])%mod

    fact=[1,1]
    factinv=[1,1]
    inv=[0,1]
    
    for i in range(2, max(N,M) + 3):
        fact.append((fact[-1] * i) % mod)
        inv.append((-inv[mod % i] * (mod // i)) % mod)
        factinv.append((factinv[-1] * inv[-1]) % mod)
    
    
    """
    Aは固定．
    
    「1つもない」は数え辛い，「少なくともk個ある」は数えやすい
    i=Piがk個以上ある
        kこ選ぶ：NCK
        残りを適当に並べる：(M-K)PK
        係数：(-1)^k
    
    """
    ans=0
    for k in range(N+1):
        cm=cmb(N,k,mod)
        rem=(fact[M-k]*factinv[M-N])%mod
        if k%2==0:
            c=1
        else:
            c=-1
        ans=(ans+c*cm*rem)%mod
        
    ans=(ans*fact[M]*factinv[M-N])%mod
    print(ans)
    
    


        
    

main()

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


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


ans_a=g1[m]*g2[m-n]


ans_b=0
for i in range(n+1):
  tmp=cmb(n,i,mod)
  tmp*=(g1[m-i]*g2[m-n])%mod
  tmp%=mod
  ans_b+=tmp*(-1)**i
  ans_b%=mod
  

print((ans_a*ans_b)%mod)

Solution (Python):
import sys
def input(): return sys.stdin.readline().strip()
def mapint(): return map(int, input().split())
sys.setrecursionlimit(10**9)

N, M = mapint()
mod = 10**9+7
pos = [1]+[0]*M
neg = [1]+[0]*M
n = 1
for i in range(1, M+1):
    n = n*i%mod
    pos[i] = n
    neg[i] = pow(n, mod-2, mod)
ans = 0

def comb(n, r):
    return pos[n]*neg[r]*neg[n-r]%mod


for i in range(N+1):
    ans += (-1)**i*comb(N, i)*pos[M]*neg[M-i]*(pow(pos[M-i]*neg[M-N], 2, mod))
    ans %= mod
print(ans)

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

n = N + M
fac = [1] * (n + 1) 
inv = [1] * (n + 1)
for j in range(1, n + 1):
    fac[j] = (fac[j - 1] * j) % MOD

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

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

def perm(n, r):
    if r > n or n < 0 or r < 0:
        return 0
    return comb(n, r) * fac[r] % MOD

ans = pow(perm(M, N), 2, MOD)

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

print(ans)

Solution (Python):

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


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

invfact = [1] * (M+1)
invfact[M] = pow(fact[M], MOD-2, MOD)
for i in range(M, 1, -1):
	invfact[i-1] = invfact[i] * i % MOD

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

def choose(n, k):
	return permute(n, k) * invfact[k] % MOD

total = pow(permute(M, N), 2, MOD)
add = -1

for i in range(1, N+1):
	total = (total + add * choose(N, i) * permute(M, i) * pow(permute(M-i, N-i), 2, MOD)) % MOD
	add *= -1

print(total)

Solution (Python):

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


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

invfact = [1] * (M+1)
invfact[M] = pow(fact[M], MOD-2, MOD)
for i in range(M, 1, -1):
	invfact[i-1] = invfact[i] * i % MOD

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

def choose(n, k):
	return permute(n, k) * invfact[k] % MOD

total = pow(permute(M, N), 2, MOD)
add = -1

for i in range(1, N+1):
	total = (total + add * choose(N, i) * permute(M, i) % MOD * pow(permute(M-i, N-i), 2, MOD)) % MOD
	add *= -1

print(total)

Solution (Python):
N,M = map(int,input().split())
mod = 10**9 + 7
leng = 500010
fac = [1,1]
for i in range(2,500010):
    fac.append(fac[-1]*i%mod)

invfac = [0]*500010
invfac[500009] = pow(fac[500009],-1,mod)
for i in range(500009):
    invfac[500008-i] = invfac[500009-i]*(500009-i)%mod

def combi(n,m):
    return(fac[n]*invfac[m]*invfac[n-m]%mod)

def permu(n,m):
    return(fac[n]*invfac[n-m]%mod)

Ans = 0

for i in range(N+1):
    if i%2 == 0:
        Ans += combi(N,i)*permu(M,i)*(permu(M-i,N-i)**2)
    else:
        Ans -= combi(N,i)*permu(M,i)*(permu(M-i,N-i)**2)

print(Ans%mod)

Solution (Python):
'''
完全順列(derangement)
モンモール数(Montmort number)
'''

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


ans = 1
for i in range(N):
    ans *= M-i
    ans %= MOD


d2, d1 = 1, M-N
d0 = d1
for i in range(2, N+1):
    
    
    d0 = (i-1)*d2 % MOD
    
    d0 += (M-N+i-1)*d1
    d0 %= MOD
    d2, d1 = d1, d0

print(ans*d0 % MOD)

Solution (Python):
'''
完全順列(derangement)
モンモール数(Montmort number)
'''

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


ans = 1
for i in range(N):
    ans *= M-i
    ans %= MOD


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

print(ans*d[-1] % MOD)

Solution (Python):
import sys 
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines
mod = 10**9+7
N, M = map(int, readline().split())
 
fact = [1] * (M+1)
fact_inv = [1] * (M+1)
for i in range(1,M+1):
    fact[i] = fact[i-1] * i % mod 
fact_inv[M] = pow(fact[M],mod-2,mod)
for i in range(M,0,-1):
    fact_inv[i-1] = (i * fact_inv[i]) % mod 
 
def comb(n,r):
  return fact[n] * fact_inv[r] * fact_inv[n-r] % mod
total = 0
for k in range(N+1):
    a = fact[M-k] * fact_inv[M-N] * comb(N,k) % mod
    if k % 2 == 0:
        total += a
    else:
        total -= a
ans = (comb(M,N) * fact[N]) % mod * total 
print(ans%mod)

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

fac = [1, 1]
for i in range(2, table_len):
    fac.append(fac[-1] * i % MOD)
    
finv = [0] * table_len
finv[-1] = pow(fac[-1], MOD - 2, MOD)
for i in range(table_len-1, -1, -1):
    finv[i-1] = finv[i] * i % MOD

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

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

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

print(sum( ((-1)**k * comb(N, k) * perm(M, k) * perm(M-k, N-k)**2 % MOD ) for k in range(N+1)) % MOD)

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

fac = [1, 1]
for i in range(2, table_len):
    fac.append(fac[-1] * i % MOD)
    
finv = [0] * table_len
finv[-1] = pow(fac[-1], MOD - 2, MOD)
for i in range(table_len-1, -1, -1):
    finv[i-1] = finv[i] * i % MOD

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

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

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

print(sum( ((-1)**k * comb(N, k) * perm(M, k) * perm(M-k, N-k)**2 % MOD ) for k in range(N+1)) % MOD)

Solution (Python):
n, m = [int(x) for x in input().split(" ")]
MOD = int(1e9 + 7)


fact = [0 for _ in range(m + 1)]
fact[0] = 1
for i in range(1, m + 1):
    fact[i] = fact[i-1] * i % MOD

finv = [0 for _ in range(m + 1)]
finv[m] = pow(fact[m], MOD-2, MOD)
for i in range(m-1, -1, -1):
    finv[i] = finv[i+1] * (i + 1) % MOD

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

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

ans = 0
for k in range(0, n + 1):
    tmp = nCk(n, k) * nPk(m, k) * (nPk(m-k, n-k) ** 2)
    if k % 2 == 0:
        ans += tmp
    else:
        ans -= tmp
    ans = ans % MOD

print(ans)

Solution (Python):

class Combi():
    def __init__(self, N, mod):
        self.power = [1 for _ in range(N+1)]
        self.rev = [1 for _ in range(N+1)]
        self.mod = mod
        for i in range(2, N+1):
            self.power[i] = (self.power[i-1]*i) % self.mod
        self.rev[N] = pow(self.power[N], self.mod-2, self.mod)
        for j in range(N, 0, -1):
            self.rev[j-1] = (self.rev[j]*j) % self.mod

    def com(self, K, R):
        if K < R:
            return 0
        else:
            return ((self.power[K])*(self.rev[K-R])*(self.rev[R])) % self.mod

    def pom(self, K, R):
        if K < R:
            return 0
        else:
            return (self.power[K])*(self.rev[K-R]) % self.mod


def main():
    N, M = map(int, input().split())
    mod = 10**9+7
    c = Combi(5*10**5, mod)
    ans = c.pom(M, N)**2 % mod
    for k in range(1, N+1):
        s = (-1)**((k-1) % 2)
        t = c.com(N, k)
        u = c.pom(M-k, N-k)**2
        v = c.pom(M, k)
        ans -= s*t*u*v
        ans = ans % mod
    print(ans)


if __name__ == "__main__":
    main()

Solution (Python):


import sys
import math
from bisect import bisect_right as br
from bisect import bisect_left as bl
sys.setrecursionlimit(2147483647)
from heapq import heappush, heappop,heappushpop
from collections import defaultdict
from itertools import accumulate
from collections import Counter
from collections import deque
from operator import itemgetter
from itertools import permutations
mod = 10**9 + 7
inf = float('inf')
def I(): return int(sys.stdin.readline())
def LI(): return list(map(int,sys.stdin.readline().split()))

def comb(a, b):
    return fact[a] * pow(fact[a-b], mod-2, mod) * pow(fact[b], mod-2, mod) % mod
n, m = LI()
fact = [1] * (max(n, m) + 1)
for i in range(1, max(n, m)+1):
    fact[i] = i * fact[i-1]
    fact[i] %= mod

ret = 0
for x in range(n+1):
    if x % 2:
        ret -= comb(n, x) * fact[m-x] * pow(fact[m-x-(n-x)], mod-2, mod)
        ret %= mod
    else:
        ret += comb(n, x) * fact[m-x] * pow(fact[m-x-(n-x)], mod-2, mod)
        ret %= mod

ans = fact[m] * pow(fact[m-n], mod-2, mod) * ret
ans %= mod
print(ans)

Solution (Python):

class Combi():
    def __init__(self, N, mod):
        self.power = [1 for _ in range(N+1)]
        self.rev = [1 for _ in range(N+1)]
        self.mod = mod
        for i in range(2, N+1):
            self.power[i] = (self.power[i-1]*i) % self.mod
        self.rev[N] = pow(self.power[N], self.mod-2, self.mod)
        for j in range(N, 0, -1):
            self.rev[j-1] = (self.rev[j]*j) % self.mod

    def com(self, K, R):
        if K < R:
            return 0
        else:
            return ((self.power[K])*(self.rev[K-R])*(self.rev[R])) % self.mod

    def pom(self, K, R):
        if K < R:
            return 0
        else:
            return (self.power[K])*(self.rev[K-R]) % self.mod


def main():
    N, M = map(int, input().split())
    mod = 10**9+7
    c = Combi(5*10**5, mod)
    ans = c.pom(M, N)**2 % mod
    for k in range(1, N+1):
        s = (-1)**((k-1) % 2)
        t = c.com(N, k)
        u = c.pom(M-k, N-k)**2
        v = c.pom(M, k)
        ans -= s*t*u*v
        ans = ans % mod
    print(ans)


if __name__ == "__main__":
    main()

Solution (Python):
class Factorial:
    def __init__(self, max_fact, mod):
        
        
        f = [1] * (max_fact + 1)
        for idx in range(2, max_fact + 1):
            f[idx] = f[idx - 1] * idx
            f[idx] %= mod
        fi = [pow(f[-1], mod - 2, mod)]
        for idx in range(max_fact, 0, -1):
            fi += [fi[-1] * idx % mod]
        fi = fi[::-1]
        self.mod = mod
        self.f = f
        self.fi = fi

    def factorial(self, n):
        return self.f[n]

    def factorial_inverse(self, n):
        return self.fi[n]

    def combination(self, n, r):
        f = self.f
        fi = self.fi
        return f[n] * fi[r] * fi[n - r] % self.mod

    def permutation(self, n, r):
        return self.f[n] * self.fi[n - r] % self.mod

    def homogeneous_product(self, n, r):
        f = self.f
        fi = self.fi
        return f[n + r - 1] * fi[r] * fi[n - 1] % self.mod


max_fact = 5*10**5
mod = 10**9 + 7
fact_instance = Factorial(max_fact, mod)
comb = fact_instance.combination
perm = fact_instance.permutation

N, M = [int(_) for _ in input().split()]
ans = 0
for p in range(N + 1):
    ans += (-1)**p * comb(N, p) * perm(M - p, N - p)
    ans %= mod
ans *= perm(M, N)
ans %= mod
print(ans)

Solution (Python):
class Factorial:
    def __init__(self, max_fact, mod):
        
        
        f = [1] * (max_fact + 1)
        for idx in range(2, max_fact + 1):
            f[idx] = f[idx - 1] * idx
            f[idx] %= mod
        fi = [pow(f[-1], mod - 2, mod)]
        for idx in range(max_fact, 0, -1):
            fi += [fi[-1] * idx % mod]
        fi = fi[::-1]
        self.mod = mod
        self.f = f
        self.fi = fi

    def factorial(self, n):
        return self.f[n]

    def factorial_inverse(self, n):
        return self.fi[n]

    def combination(self, n, r):
        f = self.f
        fi = self.fi
        return f[n] * fi[r] * fi[n - r] % self.mod

    def permutation(self, n, r):
        return self.f[n] * self.fi[n - r] % self.mod

    def homogeneous_product(self, n, r):
        f = self.f
        fi = self.fi
        return f[n + r - 1] * fi[r] * fi[n - 1] % self.mod


max_fact = 10**6
mod = 10**9 + 7
fact_instance = Factorial(max_fact, mod)
comb = fact_instance.combination
perm = fact_instance.permutation

N, M = [int(_) for _ in input().split()]
ans = 0
for p in range(N + 1):
    ans += (-1)**p * comb(N, p) * perm(M - p, N - p)
    ans %= mod
ans *= perm(M, N)
ans %= mod
print(ans)

Solution (Python):
class Factorial:
    def __init__(self, max_fact, mod):
        
        
        f = [1] * (max_fact + 1)
        for idx in range(2, max_fact + 1):
            f[idx] = f[idx - 1] * idx
            f[idx] %= mod
        fi = [pow(f[-1], mod - 2, mod)]
        for idx in range(max_fact, 0, -1):
            fi += [fi[-1] * idx % mod]
        fi = fi[::-1]
        self.mod = mod
        self.f = f
        self.fi = fi

    def factorial(self, n):
        return self.f[n]

    def factorial_inverse(self, n):
        return self.fi[n]

    def combination(self, n, r):
        f = self.f
        fi = self.fi
        return f[n] * fi[r] * fi[n - r] % self.mod

    def permutation(self, n, r):
        return self.f[n] * self.fi[n - r] % self.mod

    def homogeneous_product(self, n, r):
        f = self.f
        fi = self.fi
        return f[n + r - 1] * fi[r] * fi[n - 1] % self.mod


max_fact = 10**6
mod = 10**9 + 7
fact_instance = Factorial(max_fact, mod)
comb = fact_instance.combination
perm = fact_instance.permutation

N, M = [int(_) for _ in input().split()]
ans = 0
for p in range(N + 1):
    ans += (-1)**p * comb(N, p) * perm(M - p, N - p)
    ans %= mod
ans *= perm(M, N)
ans %= mod
print(ans)

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

mod = 1000000007
fact = [1] * (m + 1)
ifact = [1] * (m + 1)

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


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


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


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

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

md=1000000007
fac=[0]*(M+1)
inv=[0]*(M+1)
ifac=[0]*(M+1)
inv[1]=1
for i in range(2,M+1):
	inv[i] = (md - (md//i) * inv[md%i] % md) % md;
fac[0]=1
ifac[0]=1
for i in range(1,M+1):
    fac[i]=fac[i-1]*i%md;
    ifac[i]=ifac[i-1]*inv[i]%md;

def choose(n,k):
    return fac[n]*ifac[k]*ifac[n-k]%md
def perm(m,n):
    return fac[m]*ifac[m-n]%md
ans=0
for k in range(N+1):
    if k%2==0:
        ans+=choose(N,k)*perm(M,k)*perm(M-k,N-k)**2
        ans%=md
    else:
        ans-=choose(N,k)*perm(M,k)*perm(M-k,N-k)**2
        ans%=md
print(ans)

Solution (Python):
import sys

read = sys.stdin.read
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 9)
INF = 1 << 60
MOD = 1000000007


def main():
    N, M = map(int, readline().split())

    COM_MAX = M

    fac, finv, inv = [0] * (COM_MAX + 1), [0] * (COM_MAX + 1), [0] * (COM_MAX + 1)
    fac[0] = fac[1] = finv[0] = finv[1] = inv[1] = 1

    for i in range(2, COM_MAX + 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

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

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

    ans = 0
    for k in range(N + 1):
        tmp = perm(M - k, N - k) * com(N, k) % MOD
        if k & 1:
            tmp *= -1
        ans = (ans + tmp) % MOD

    ans = ans * perm(M, N) % MOD

    print(ans)
    return


if __name__ == '__main__':
    main()

Solution (Python):
import sys

read = sys.stdin.read
readline = sys.stdin.readline
readlines = sys.stdin.readlines
sys.setrecursionlimit(10 ** 9)
INF = 1 << 60
MOD = 1000000007


def main():
    N, M = map(int, readline().split())

    COM_MAX = M

    fac, finv, inv = [0] * (COM_MAX + 1), [0] * (COM_MAX + 1), [0] * (COM_MAX + 1)
    fac[0] = fac[1] = 1
    finv[0] = finv[1] = 1
    inv[1] = 1

    for i in range(2, COM_MAX + 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

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

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

    ans = 0
    for k in range(N + 1):
        tmp = perm(M - k, N - k) * com(N, k) % MOD
        if k & 1:
            tmp *= -1
        ans = (ans + tmp) % MOD

    ans = ans * perm(M, N) % MOD

    print(ans)
    return


if __name__ == '__main__':
    main()

Solution (Python):
import bisect
import collections
import copy
import functools
import heapq
import math
import sys
from collections import deque
from collections import defaultdict
input = sys.stdin.readline
sys.setrecursionlimit(10**9)
MOD = 10**9+7

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

kaijo_memo = []
def kaijo(n):
  if(len(kaijo_memo) > n):
    return kaijo_memo[n]
  if(len(kaijo_memo) == 0):
    kaijo_memo.append(1)
  while(len(kaijo_memo) <= n):
    kaijo_memo.append(kaijo_memo[-1] * len(kaijo_memo) % MOD)
  return kaijo_memo[n]
 
gyaku_kaijo_memo = []
def gyaku_kaijo(n):
  if(len(gyaku_kaijo_memo) > n):
    return gyaku_kaijo_memo[n]
  if(len(gyaku_kaijo_memo) == 0):
    gyaku_kaijo_memo.append(1)
  while(len(gyaku_kaijo_memo) <= n):
    gyaku_kaijo_memo.append(gyaku_kaijo_memo[-1] * pow(len(gyaku_kaijo_memo),MOD-2,MOD) % MOD)
  return gyaku_kaijo_memo[n]
 
def nCr(n,r):
  if(n == r):
    return 1
  if(n < r or r < 0):
    return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(r) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret

def nPr(n,r):
  if(n < r or r < 0):
    return 0
  ret = 1
  ret = ret * kaijo(n) % MOD
  ret = ret * gyaku_kaijo(n-r) % MOD
  return ret

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

Solution (Python):
import sys

def input(): return sys.stdin.readline().strip()
def list2d(a, b, c): return [[c] * b for i in range(a)]
def list3d(a, b, c, d): return [[[d] * c for j in range(b)] for i in range(a)]
def list4d(a, b, c, d, e): return [[[[e] * d for j in range(c)] for j in range(b)] for i in range(a)]
def ceil(x, y=1): return int(-(-x // y))
def INT(): return int(input())
def MAP(): return map(int, input().split())
def LIST(N=None): return list(MAP()) if N is None else [INT() for i in range(N)]
def Yes(): print('Yes')
def No(): print('No')
def YES(): print('YES')
def NO(): print('NO')
sys.setrecursionlimit(10 ** 9)
INF = 10 ** 19
MOD = 10 ** 9 + 7
EPS = 10 ** -10

class ModTools:
    """ 階乗・逆元用のテーブルを構築する """

    def __init__(self, MAX, MOD):

        
        MAX += 1
        self.MAX = MAX
        self.MOD = MOD
        
        factorial = [1] * MAX
        factorial[0] = factorial[1] = 1
        for i in range(2, MAX):
            factorial[i] = factorial[i-1] * i % MOD
        
        inverse = [1] * MAX
        
        inverse[MAX-1] = pow(factorial[MAX-1], MOD-2, MOD)
        for i in range(MAX-2, -1, -1):
            
            inverse[i] = inverse[i+1] * (i+1) % MOD
        self.fact = factorial
        self.inv = inverse

    def nCr(self, n, r):
        """ 組み合わせ """

        if n < r: return 0
        
        r = min(r, n-r)
        
        numerator = self.fact[n]
        
        denominator = self.inv[r] * self.inv[n-r] % self.MOD
        return numerator * denominator % self.MOD

    def nHr(self, n, r):
        """ 重複組み合わせ """

        
        return self.nCr(r+n-1, r)

    def nPr(self, n, r):
        """ 順列 """

        if n < r: return 0
        return self.fact[n] * self.inv[n-r] % self.MOD

N, M = MAP()

mt = ModTools(M+1, MOD)
acnt = mt.nPr(M, N)

bcnt = 0
for i in range(N+1):
    bcnt += mt.nCr(N, i) * mt.nPr(M-i, N-i) * (-1)**i
    bcnt %= MOD
ans = acnt * bcnt % MOD
print(ans)

Solution (Python):


N, M = map(int, input().split())
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)

fixed = (fac[N] * fac[M] * pow(finv[M-N], 2, MOD)) % MOD

def f(k):
    return (fac[M-k] * finv[k] * finv[N-k]) % MOD

ans = 0
for i in range(N+1):
    ans = (ans + f(i) * pow(-1, i)) % MOD

print((fixed * ans) % MOD)

Solution (Python):


N, M = map(int, input().split())
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)

fixed = (fac[N] * fac[M] * pow(finv[M-N], 2, MOD)) % MOD

def f(k):
    return (fac[M-k] * finv[k] * finv[N-k]) % MOD

ans = 0
for i in range(N+1):
    ans = (ans + fixed * f(i) * pow(-1, i)) % MOD

print(ans)

Solution (Python):


N, M = map(int, input().split())
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)

fixed = (fac[N] * fac[M] * pow(finv[M-N], 2, MOD)) % MOD

def f(k):
    return (fac[M-k] * finv[k] * finv[N-k]) % MOD

ans = 0
for i in range(N+1):
    ans = (ans + fixed * f(i) * pow(-1, i)) % MOD

print(ans)

Solution (Python):

import sys

sysread = sys.stdin.readline
read = sys.stdin.read
sys.setrecursionlimit(10 ** 7)





import numpy as np

from collections import deque


def generate_inv(n,mod):
    """
    逆元行列
    n >= 2
    Note: mod must bwe a prime number
    """
    ret = [0, 1]
    for i in range(2,n+1):
        next = -ret[mod%i] * (mod // i)
        next %= mod
        ret.append(next)
    return ret

def generate_comb_m_i(M, N, inv, mod):
    ret = [1, M]
    for i in range(2, N+1):
        tmp = (ret[-1] * (M-(i-1)) * inv[i]) % mod
        ret.append(tmp)
    return ret
def generate_p_N_i(N, mod):
    ret = [1]
    for i in range(1, N+1):
        tmp = (ret[-1] * (N - (i-1))) % mod
        ret.append(tmp)
    return ret

def generate_p_mi_ni(M, N, mod):
    ret = deque([1])
    cache = deque([1])
    for i in range(1, N+1):
        tmp = (cache[0] * (M-N+i)) % mod
        cache.appendleft(tmp)
        ret.appendleft((tmp ** 2) % mod)
    return ret


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

    inv = generate_inv(M, mod)

    comb_m_i = generate_comb_m_i(M, N, inv, mod)
    
    p_N_i = generate_p_N_i(N, mod)
    
    p_mi_ni = generate_p_mi_ni(M, N, mod)
    

    ret = p_mi_ni[0]

    arr01 = [-1, 1] * (N//2 + 1)
    arr01 = arr01[:N+1]

    comb_m_i = np.array(comb_m_i, dtype=np.int64)
    p_N_i = np.array(p_N_i, dtype=np.int64)
    p_mi_ni = np.array(p_mi_ni, dtype=np.int64)
    arr01 = np.array(arr01, dtype=np.int64)

    tmp = (comb_m_i * p_N_i) % mod
    tmp = (tmp * p_mi_ni) % mod
    tmp = (tmp * arr01) % mod
    tmp = tmp[1:]
    sub = 0
    for val in tmp:
        sub += val
        sub %= mod

    ret -= sub
    ret %= mod

    print(ret)

if __name__ == "__main__":
    run()

Solution (Python):
n, m = map(int, input().split())
mod = 10**9+7
nums = 5*(10**5) 
 
g1, g2, inverse = [1, 1] , [1, 1], [0, 1]
 
for num in range(2, nums + 1):
    g1.append((g1[-1] * num) % mod)
    inverse.append((-inverse[mod % num] * (mod//num)) % mod)
    g2.append((g2[-1] * inverse[-1]) % mod)
    
def nPr(n, r):
    return g1[n] * g2[n-r] % mod
 
def nCr(n, r):
    r = min(r, n-r)
    return g1[n] * g2[r] * g2[n-r] % mod
 
ans = 0
for r in range(n+1):
    tmp = nCr(n, r)*nPr(m, r)*(nPr(m-r, n-r)**2)%mod
    ans = (ans+(-1)**r*tmp)%mod
print(ans)

Solution (Python):
def makefac(n):
    global faclist
    faclist = [1, 1]
    for i in range(2, n + 1):
        faclist.append(faclist[-1] * i % MOD)
    return None
def mcomb(n, r):
    return faclist[n] * pow(faclist[r],MOD - 2, MOD) * pow(faclist[n-r], MOD - 2, MOD)
def mperm(n, r):
    return faclist[n] * pow(faclist[n-r],MOD - 2, MOD)

N, M = map(int, input().split())
MOD = 10 ** 9 + 7
makefac(M)
ans = mperm(M, N)
for k in range(1, N + 1):
    ans += ((-1)**k * mcomb(N, k) * mperm(M - k, N - k)) % MOD
ans = ans * mperm(M, N) % MOD
print(ans)

Solution (Python):
def makefac(n):
    global faclist
    faclist = [1, 1]
    for i in range(2, n + 1):
        faclist.append(faclist[-1] * i % MOD)
    return None
def mcomb(n, r):
    return faclist[n] * pow(faclist[r],MOD - 2, MOD) * pow(faclist[n-r], MOD - 2, MOD)
def mperm(n, r):
    return faclist[n] * pow(faclist[n-r],MOD - 2, MOD)

N, M = map(int, input().split())
MOD = 10 ** 9 + 7
makefac(M)
ans = mperm(M, N)
for k in range(1, N + 1):
    ans += ((-1)**k * mcomb(N, k) * mperm(M - k, N - k)) % MOD
ans = ans * mperm(M, N) % MOD
print(ans)

Solution (Python):

class Combination:
    def __init__(self, n_max, mod=10 ** 9 + 7):
        
        self.mod = mod
        f = 1
        self.fac = fac = [f]
        for i in range(1, n_max + 1):
            f = f * i % mod
            fac.append(f)
        f = pow(f, mod - 2, mod)
        self.facinv = facinv = [f]
        for i in range(n_max, 0, -1):
            f = f * i % mod
            facinv.append(f)
        facinv.reverse()

    
    

    def __call__(self, n, r):  
        return self.fac[n] * self.facinv[r] % self.mod * self.facinv[n - r] % self.mod

    def nCr(self, n, r):
        if not 0 <= r <= n:
            return 0
        return self.fac[n] * self.facinv[r] % self.mod * self.facinv[n - r] % self.mod

    def nPr(self, n, r):
        if not 0 <= r <= n:
            return 0
        return self.fac[n] * self.facinv[n - r] % self.mod

MOD = 10**9+7
def resolve():
    N, M = map(int, input().split())
    cmb = Combination(5 * 10**5)

    ans = 0
    for p in range(N+1):
        ans += (-1)**p * cmb.nCr(N, p) * cmb.nPr(M-p, N-p)
        ans %= MOD

    ans *= cmb.nPr(M, N)
    ans %= MOD

    print(ans)


if __name__ == "__main__":
    resolve()

Solution (Python):
class Factorial:
    def __init__(self, max_fact, mod):
        
        
        f = [1] * 2 + [None] * (max_fact - 1)
        for idx in range(2, max_fact + 1):
            f[idx] = f[idx - 1] * idx % mod
        fi = [None] * max_fact + [pow(f[-1], mod - 2, mod)]
        for idx in range(max_fact, 0, -1):
            fi[idx - 1] = fi[idx] * idx % mod
        self.mod = mod
        self.f = f
        self.fi = fi

    def factorial(self, n):
        return self.f[n]

    def factorial_inverse(self, n):
        return self.fi[n]

    def combination(self, n, r):
        f = self.f
        fi = self.fi
        return f[n] * fi[r] * fi[n - r] % self.mod

    def permutation(self, n, r):
        return self.f[n] * self.fi[n - r] % self.mod

    def homogeneous_product(self, n, r):
        f = self.f
        fi = self.fi
        return f[n + r - 1] * fi[r] * fi[n - 1] % self.mod


max_fact = 5 * 10**5
mod = 10**9 + 7
fact_instance = Factorial(max_fact, mod)
comb = fact_instance.combination
perm = fact_instance.permutation

N, M = [int(_) for _ in input().split()]
ans = 0
for p in range(N + 1):
    ans += (-1)**p * comb(N, p) * perm(M - p, N - p)
    ans %= mod
ans *= perm(M, N)
ans %= mod
print(ans)