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
import time
import math
from collections import deque
import heapq
import itertools
from decimal import Decimal
import bisect
from operator import itemgetter
MAX_INT = int(10e18)
MIN_INT = -MAX_INT
mod = 1000000000+7
sys.setrecursionlimit(1000000)
def IL(): return list(map(int,input().split()))
def SL(): return input().split()
def I(): return int(sys.stdin.readline())
def S(): return input()

import functools
@functools.lru_cache(maxsize=None)
def factrial_memo(n):
  fact = [1, 1]
  for i in range(2, n + 1):
    fact.append((fact[-1] * i) % mod)
  return fact

def permutation(n,r): 
  return fact[n]*pow(fact[n-r],mod-2,mod)%mod
  
def combination(n,r): 
  return permutation(n,r)*pow(fact[r],mod-2,mod)%mod

def homogeneous(n,r): 
  return combination(n+r-1,r)%mod

N,M = IL()

fact = factrial_memo(10**6)

s = 0
for k in range(1,N+1):
  s += combination(N, k) * permutation(M-k, N-k) * (-1)**(k-1)
ans = (permutation(M, N)**2 - s*permutation(M, N)) %mod
print(ans)

Solution (Python):
SIZE=5*10**5+1; MOD=10**9+7 
 
SIZE += 1
inv = [0]*SIZE  
fac = [0]*SIZE  
finv = [0]*SIZE 
inv[1] = 1
fac[0] = fac[1] = 1
finv[0] = finv[1] = 1
for i in range(2,SIZE):
    inv[i] = MOD - (MOD//i)*inv[MOD%i]%MOD
    fac[i] = fac[i-1]*i%MOD
    finv[i]= finv[i-1]*inv[i]%MOD
 
def choose(n,r): 
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0
 
def chofuku(n,r): 
    return choose(n+r-1,r)

def narabekae(n,r): 
    return fac[n]*finv[n-r]%MOD


N,M=map(int,input().split())
ans=0
for k in range(N+1):
    ans+=choose(N,k)*narabekae(M-k,N-k)*(-1 if k%2==1 else 1)
    ans%=MOD
ans*=narabekae(M,N)
print(ans%MOD)

Solution (Python):
import sys
read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

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

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 C(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 P(self, n, r):
        if not 0 <= r <= n: return 0
        return self.fac[n] * self.facinv[n-r] % self.mod
      
f = Combination(m, mod)
ans = 0
for k in range(m+1):
  chk = f.C(n,k)*f.P(m,k)*pow(f.P(m-k,n-k),2,mod)%mod
  ans += pow(-1,k)*chk
  ans %= mod
print(ans)

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

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

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

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 )

ans = 0
for k in range(n+1):
    a = cmb(n,k,mod)
    p1 = perm(m,k)
    p2 = perm(m-k, n-k)
    ans += (a * (-1)**k * p1 * p2**2)

print(ans%mod)

Solution (Python):
N,M = [int(i) for i in input().split()]
P = 10**9+7
nn = M+1
fa = [1] * (nn+1)
fainv = [1] * (nn+1)
for i in range(nn):
    fa[i+1] = fa[i] * (i+1) % P
fainv[-1] = pow(fa[-1], P-2, P)
for i in range(nn)[::-1]:
    fainv[i] = fainv[i+1] * (i+1) % P
C = lambda a, b: fa[a] * fainv[b] % P * fainv[a-b] % P if 0 <= b <= a else 0

ans = ((fa[M]*fainv[M-N]%P)**2)%P
for i in range(1,N+1):
  ans += (-1)**i * C(N,i) * fa[M]%P * fa[M-i]%P * fainv[M-N]%P * fainv[M-N]%P
  ans %= P
print(ans)

Solution (Python):
N,M = [int(i) for i in input().split()]
P = 10**9+7
nn = 500000
fa = [1] * (nn+1)
fainv = [1] * (nn+1)
for i in range(nn):
    fa[i+1] = fa[i] * (i+1) % P
fainv[-1] = pow(fa[-1], P-2, P)
for i in range(nn)[::-1]:
    fainv[i] = fainv[i+1] * (i+1) % P
C = lambda a, b: fa[a] * fainv[b] % P * fainv[a-b] % P if 0 <= b <= a else 0

ans = ((fa[M]*fainv[M-N]%P)**2)%P
for i in range(1,N+1):
  ans += (-1)**i * C(N,i) * fa[M]%P * fa[M-i]%P * fainv[M-N]%P * fainv[M-N]%P
  ans %= P
print(ans)

Solution (Python):

"""
玉：各要素(1~m)(区別あり)
箱：配列の格納場所(1~n)(区別あり)
条件：全射かつ、任意のiにおいてAi != Biとなるような2つの組み合わせ
問題の余事象は
Ai = Biとなるような要素の選び方
となる。
全事象(mCn ^2)からこれを引けば良いので、
包除原理を用いて、
A0 = B0となる場合
A1 = B1となる場合…
を数え上げる。

0~n個までの重複要素について、
-1^i *(n箱からi個選んで重複させる*その時のm個の整数からの選び方はmPi)*残りのn-i個はm-i個から任意に選ぶ(m-iPn-i)^2
を実装する
n箱からi箇所選ぶのは区別しなくて良いが、(その選び方に区別は無いため)
それ以外の、その場所にどの数字を入れるかはちゃんと区別しないといけないため順列になり、
また、残りの箇所もAとBの配列それぞれについて(ここが^2)整数を区別して(順列)自由に選べることに注意。
"""
import sys
readline = sys.stdin.buffer.readline
def even(n): return 1 if n%2==0 else 0
mod = 10**9+7
n,m = map(int,readline().split())
def pow(n,p,mod=mod): 
    res = 1
    while p > 0:
        if p % 2 == 0:
            n = n ** 2 % mod
            p //= 2
        else:
            res = res * n % mod
            p -= 1
    return res % mod

def factrial_memo(n=10**6,mod=mod):
    fact = [1, 1]
    for i in range(2, n + 1):
        fact.append((fact[-1] * i) % mod)
    return fact

fact = factrial_memo()

def permutation(n,r): 
    return fact[n]*pow(fact[n-r],mod-2)%mod
def combination(n,r): 
    return permutation(n,r)*pow(fact[r],mod-2)%mod
    
def homogeneous(n,r): 
    return combination(n+r-1,r)%mod
    

ans = 0
for i in range(n+1):
    
    ans += pow(-1,i)*combination(n,i)*permutation(m,i)*pow(permutation(m-i,n-i),2)%mod
    ans %= mod

print(ans)

Solution (Python):
def combs_mod(n,k,mod):
    
    inv = [1]*(k+1)
    for i in range(1,k+1):
        inv[i] = pow(i,mod-2,mod)
    ans = [1]*(k+1)
    for i in range(1,k+1):
        ans[i] = ans[i-1]*(n+1-i)*inv[i]%mod
    return ans

def solve():
    ans = 0
    mod = 10**9+7
    N, M = map(int, input().split())
    NCi = combs_mod(N,N,mod)
    MPi = [1]*(N+1)
    for i in range(1,N+1):
        MPi[i] = MPi[i-1]*(M+1-i)%mod
    MiPNi = [1]*(N+1)
    for i in range(N-1,-1,-1):
        MiPNi[i] = MiPNi[i+1]*(M-i)%mod
    for i in range(N+1):
        ans += pow(-1,i)*NCi[i]*MPi[i]*(MiPNi[i]**2)%mod
        ans %= mod
    return ans
print(solve())

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

mod = 10**9 + 7
l = 5*10**5 + 3

M = [1]  
mul = 1
for i in range(1, l):
    mul = (mul * i) % mod
    M.append(mul)

def pow(x, y, mod):  
    ans = 1
    while y > 0:
        if y % 2 == 1:
            ans = (ans * x) % mod
        x = (x**2) % mod
        y //= 2
    return ans

def inv(x, mod):  
    return pow(x, mod-2, mod)



MI = [0] * (l-1) +[inv(M[l-1], mod)]  
for i in range(l-2, -1, -1):
    MI[i] = MI[i+1] * (i+1) % mod

def C(x, y):  
    if y < 0 or y > x:
        return 0
    elif x > l:  
        y = min(y, x-y)
        ans = 1
        for i in range(x, x-y, -1):
            ans = (ans * i) % mod
        return (ans * MI[y]) % mod
    else:  
        ans = M[x]
        ans = (ans * MI[y]) % mod
        return (ans * MI[x-y]) % mod

def P(x, y):  
    if y < 0 or y > x:
        return 0
    elif x > l:  
        y = min(y, x-y)
        ans = 1
        for i in range(x, x-y, -1):
            ans = (ans * i) % mod
        return ans % mod
    else:  
        ans = M[x]
        return (ans * MI[x-y]) % mod

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

Solution (Python):
n,m=map(int,input().split())
mod=10**9+7
ans=0
a=1
b=1
k=1
for i in range(m-n+1,m+1):
    b*=i
    b%=mod
for i in range(n+1):
    ans+=k*a*(b**2)
    ans%=mod
    a*=(n-i)*(m-i)*pow(i+1,mod-2,mod)
    a%=mod
    b*=pow(m-i,mod-2,mod)
    b%=mod
    k*=-1
print(ans)

Solution (Python):
N, M = map(int, input().split())
MOD = 10**9+7
n = max(N, M)
fac = [1]*(n+1)
rev = [1]*(n+1)
 
for i in range(1,n+1):
  fac[i] = i*fac[i-1]%MOD
  rev[i] = pow(fac[i], MOD-2, MOD)
 
comb = lambda a,b:(fac[a]*rev[a-b]*rev[b])%MOD
perm = lambda a,b:(fac[a]*rev[a-b])%MOD

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

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

N = m + 100
MOD = 10 ** 9 + 7
fact = [0 for _ in range(N)]
invfact = [0 for _ in range(N)]
fact[0] = 1
for i in range(1, N):
    fact[i] = i * fact[i - 1] % MOD

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

for i in range(N - 2, -1, -1):
    invfact[i] = invfact[i + 1] * (i + 1) % MOD
    
def nCk(n, k):
    if k < 0 or n < k:
        return 0
    else:
        return fact[n] * invfact[k] * invfact[n - k] % MOD
        
def nPk(n, k):
    return fact[n] * invfact[n - k] % MOD
    
ans = 0
sign = 1

for i in range(n + 1):
    ans = (ans + nPk(m, i) * nCk(n, i) * nPk(m - i, n - i) ** 2 * sign) % MOD
    sign *= -1
    
print(ans)

Solution (Python):
MOD = 10 ** 9 + 7
n, m = map(int, input().split())
a = 1
d = [1] * (n + 1)
for i in range(n):
  d[i + 1] = ((m - n + i) * d[i] + i * d[i - 1]) % MOD
  a = a * (m - i) % MOD
print(a * d[n] % MOD)

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

MOD_t_MAX = 5*10**5+10

fac  = [None] * MOD_t_MAX
finv = [None] * MOD_t_MAX
inv  = [None] * MOD_t_MAX
def MOD_COM_init():
    fac[0] = fac[1] = 1
    finv[0] = finv[1] = 1
    inv[1] = 1
    for i in range(2, MOD_t_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 MOD_COM(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 MOD_perm(n, r):
    rtn = 1
    for _ in range(r):
        rtn *= n
        rtn %= MOD
        n -= 1
    return rtn

def main():
    n, m = map(int, input().split())
    ans = 0
    perm = MOD_perm(m, n)
    MOD_COM_init()
    for i in range(n+1):
        tmp = perm * MOD_COM(n, i) % MOD
        ans += (-1)**i * tmp
        ans %= MOD
        if i < n:
            perm *= inv[m-i]
            perm %= MOD
    print(ans*MOD_perm(m, n)%MOD)

if __name__ == "__main__":
    main()

Solution (Python):
import sys


MOD = 10**9 + 7  
N_MAX = 5 * 10**5  

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

for i in range(2, N_MAX + 1):
    fac.append((fac[-1] * i) % MOD)
    inv.append((-inv[MOD % i] * (MOD // i)) % MOD)
    facinv.append((facinv[-1] * inv[-1]) % MOD)


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


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


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

    all = prm(M, N) ** 2 % MOD

    dup = 0

    for i in range(1, N + 1):

        tmp = prm(M, i) * prm(M - i, N - i) ** 2 % MOD
        tmp *= cmb(N, i) * (-1) ** (i - 1)

        dup += tmp % MOD

    print((all - dup) % MOD)

main()

Solution (Python):
n,m=map(int,input().split())
a = 1
Mod=10**9+7
d = [1]*(n+1)
for i in range(n):
  d[i+1]=((m-n+i)*d[i]+i*d[i-1])%Mod
  a=a*(m-i)%Mod
print(a*d[n]%Mod)

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

    
    
    mPn = M
    for i in range(M-N+1,M):
        mPn = mPn * i % P
    

    '''
    
    
    
    
    
    
    

    '''
    num_b = 0
    nCk = 1
    m_kPn_k = mPn

    for k in range(1,N+1):
        nCk = nCk * (N-k+1) * pow(k,P-2,P) % P
        m_kPn_k = m_kPn_k * pow(M-k+1,P-2,P) % P
        num_b = (num_b + ((-1) ** (k-1)) * nCk * m_kPn_k) % P
    print((mPn ** 2 - mPn * num_b)% P)
main()

Solution (Python):
class Factorial:
    def __init__(self, n, mod):
        self.f = [1]
        self.mod = mod
        for j in range(1, n + 1):
            self.f.append(self.f[-1] * j % mod)
        self.i = [pow(self.f[-1], mod - 2, mod)]
        for j in range(n, 0, -1):
            self.i.append(self.i[-1] * j % mod)
        self.i.reverse()
    def factorial(self, j):
        return self.f[j]
    def ifactorial(self, j):
        return self.i[j]
    def comb(self, n, k):
        return self.f[n] * self.i[n - k] % self.mod * self.i[k] % self.mod if n >= k else 0


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

MOD = 10 ** 9 + 7

F = Factorial(M + 1, MOD)

ans = 0
for k in range(N + 1):
    tmp1 = F.comb(N, k)
    tmp2 = (F.factorial(M) * F.ifactorial(M - k)) % MOD
    tmp3 = ((F.factorial(M - k) * F.ifactorial(M - N))) ** 2 % MOD
    tmp = (tmp1 * tmp2 * tmp3) % MOD
    if k % 2 == 0:
        ans += tmp
    else:
        ans -= tmp
    ans %= MOD

print (ans)

Solution (Python):
import sys

readline = sys.stdin.readline
MOD = 10 ** 9 + 7
INF = float('INF')
sys.setrecursionlimit(10 ** 5)


class ModInt:
    MOD = 10 ** 9 + 7

    def __init__(self, x):
        self.x = x % MOD

    def __str__(self):
        return str(self.x)

    def __repr__(self):
        return str(self.x)

    def __add__(self, other):
        return (
            ModInt(self.x + other.x) if isinstance(other, ModInt) else
            ModInt(self.x + other)
        )

    def __sub__(self, other):
        return (
            ModInt(self.x - other.x) if isinstance(other, ModInt) else
            ModInt(self.x - other)
        )

    def __mul__(self, other):
        return (
            ModInt(self.x * other.x) if isinstance(other, ModInt) else
            ModInt(self.x * other)
        )

    def __truediv__(self, other):
        return (
            ModInt(
                self.x * pow(other.x, MOD - 2, MOD)
            ) if isinstance(other, ModInt) else
            ModInt(self.x * pow(other, MOD - 2, MOD))
        )

    def __pow__(self, other):
        return (
            ModInt(pow(self.x, other.x, MOD)) if isinstance(other, ModInt) else
            ModInt(pow(self.x, other, MOD))
        )

    __radd__ = __add__

    def __rsub__(self, other):
        return (
            ModInt(other.x - self.x) if isinstance(other, ModInt) else
            ModInt(other - self.x)
        )

    __rmul__ = __mul__

    def __rtruediv__(self, other):
        return (
            ModInt(
                other.x * pow(self.x, MOD - 2, MOD)
            ) if isinstance(other, ModInt) else
            ModInt(other * pow(self.x, MOD - 2, MOD))
        )

    def __rpow__(self, other):
        return (
            ModInt(pow(other.x, self.x, MOD)) if isinstance(other, ModInt) else
            ModInt(pow(other, self.x, MOD))
        )

    def set_mod_divisor(self, divisor):
        ModInt.MOD = divisor


def comb_mod(n, r):
    res = 1
    r = min(n - r, r)

    for i in range(r):
        res *= (n - i)
        res %= MOD
        res *= pow((r - i), MOD - 2, MOD)

    return res


def main():
    n, m = map(int, readline().split())
    comb = ModInt(0)
    p1 = ModInt(1)  
    p2 = ModInt(1)  
    c1 = ModInt(1)  
    sign = -1

    for i in range(n):
        p2 *= (m - i)

    ans = p2 * p2

    for i in range(1, n + 1):
        p1 *= (m - i + 1)
        p2 *= pow(m - i + 1, MOD - 2, MOD)
        c1 *= (n - i + 1)
        c1 *= pow(i, MOD - 2, MOD)
        sign *= (-1)
        comb += p1 * p2 * p2 * c1 * sign

    ans -= comb

    print(ans)



if __name__ == '__main__':
    main()

Solution (Python):
import sys, re
from collections import deque, defaultdict, Counter
from math import ceil, sqrt, hypot, factorial, pi, sin, cos, radians, gcd
from itertools import accumulate, permutations, combinations, product, groupby, combinations_with_replacement
from operator import itemgetter, mul
from copy import deepcopy
from string import ascii_lowercase, ascii_uppercase, digits
from bisect import bisect, bisect_left
from heapq import heappush, heappop
from functools import reduce
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)))
sys.setrecursionlimit(10 ** 9)
INF = float('inf')
mod = 10 ** 9 + 7

N, M = MAP()

lim = M+10  
fact = [1] * (lim+1)
fact_inv = [1] * (lim+1)
for n in range(1, lim+1):
    fact[n] = (fact[n-1] * n) % mod
fact_inv[lim] = pow(fact[lim], mod-2, mod)
for n in range(lim, 0, -1):
    fact_inv[n-1] = (n * fact_inv[n]) % mod
def C(n, r):
    return (((fact[n] * fact_inv[r]) % mod) * fact_inv[n-r]) % mod

def P(n, r):
    return fact[n]*fact_inv[n-r]%mod

ans = 0
for k in range(N+1):
    ans += (-1)**(k%2) * C(N, k) * P(M-k, N-k)
    ans %= mod

ans *= P(M, N)
print(ans%mod)

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

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

N = 10**6  
fact = [1, 1]  
factinv = [1, 1]  
inv = [0, 1]  

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

"""
1つ当たっている→Freeはn-1個→m-1 C n-1
同様に2つ→m-2 C n-2
数は？→ n C 2
包除原理
"""
mode = -1
ans = 0
for i in range(n+1):
	mode = -mode
	ans = (ans + mode * cmb(n, i, mod) * cmb(m-i, n-i, mod) * fact[n-i]) % mod
print(ans * cmb(m, n, mod) * fact[n] % mod)

Solution (Python):
import sys
sys.setrecursionlimit(300000)
class BinomialCoefficient(object):

    def __init__(self, N=510000, MOD=10**9+7):
        self.fac = [1, 1]
        self.finv = [1, 1]
        self.inv = [0, 1]
        self.MOD = MOD
        for i in range(2, N + 1):
            self.fac.append((self.fac[-1] * i) % MOD)
            self.inv.append((-self.inv[MOD % i] * (MOD // i)) % MOD)
            self.finv.append((self.finv[-1] * self.inv[-1]) % MOD)

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

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


def I(): return int(sys.stdin.readline())
def MI(): return map(int, sys.stdin.readline().split())
def MI0(): return map(lambda s: int(s) - 1, sys.stdin.readline().split())
def LMI(): return list(map(int, sys.stdin.readline().split()))
def LMI0(): return list(map(lambda s: int(s) - 1, sys.stdin.readline().split()))
MOD = 10 ** 9 + 7
INF = float('inf')


N, M = MI()
bi = BinomialCoefficient()

ans = 0
for i in range(0, N + 1):
    sgn = (-1) ** i
    lans = bi.C(N, i) * bi.P(M - i, N - i) % MOD
    lans = lans * sgn
    ans = (ans + lans + MOD) % MOD

ans = ans * bi.P(M, N) % MOD
print(ans)

Solution (Python):
MOD = 10**9 + 7
fac = [1,1] 
inv = [0,1] 
finv = [1,1] 

def finv_init(n, m):
    for i in range(2, n+1):
        fac.append( (fac[-1] * i ) % m )
        inv.append( (-inv[m%i] * (m//i)) % m )
        finv.append( (finv[-1] * inv[-1]) % m )

def nPr(n, r, m): return fac[n] * finv[n-r] % m

n, m = map(int, input().split())
finv_init(m, MOD)
ans = 0
for d in range(0, n+1):
    a = nPr(m-d, n-d, MOD) * nPr(n, d, MOD) * finv[d] % MOD
    if d & 1: a = -a
    ans += a; ans %= MOD
ans *= nPr(m, n, MOD)
print(ans % MOD)

Solution (Python):
MOD = 10**9 + 7
fac = [1,1] 
inv = [0,1] 
finv = [1,1] 

def finv_init(n, m):
    for i in range(2, n+1):
        fac.append( (fac[-1] * i ) % m )
        inv.append( (-inv[m%i] * (m//i)) % m )
        finv.append( (finv[-1] * inv[-1]) % m )

def nPr(n, r, m): return fac[n] * finv[n-r]

n, m = map(int, input().split())
finv_init(m, MOD)
ans = 0
for d in range(0, n+1):
    a = nPr(m-d, n-d, MOD) * nPr(n, d, MOD) * finv[d] % MOD
    if d & 1: a = -a
    ans += a; ans %= MOD
ans *= nPr(m, n, MOD)
print(ans % MOD)

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 perm(n,r,p):
    if (r < 0) or (n < r):
        return 0
    return fact[n] * factinv[n-r] % p 

p = 10 ** 9 + 7
N = 5 * 10 ** 6 + 10 
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):
    dat=perm(m,n,p)*perm(m-i,n-i,p)*cmb(n,i,p)%p
    
    if i%2==0:
        ans+=dat
    else:
        ans-=dat

print(ans%p)

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

def plmi(k):
    if k%2 == 0:
        return 1
    else:
        return -1


def mod_factorials(n, m):
    mod_factorials = [1] * (n + 1)
    for i in range(1, n+1):
        mod_factorials[i] = mod_factorials[i-1]*i
        mod_factorials[i] %= m
    return mod_factorials


def mod_inv_factorials(n, m):
    mod_invs = [1] * (n + 1)
    invs = [1] * (n + 1)
    for i in range(2, n + 1):
        invs[i] = m - invs[m % i] * (m//i) % m

    for i in range(1, n+1):
        mod_invs[i] = invs[i]*mod_invs[i-1]
        mod_invs[i] %= m
    return mod_invs

mod_facts = mod_factorials(M,mod)
mod_inv_facts = mod_inv_factorials(M,mod)

ans = 0

for i in range(0,N+1):
    tmp = 1
    tmp *= plmi(i)*mod_facts[N]*mod_facts[M-i]*mod_inv_facts[i]*mod_inv_facts[N-i]*mod_inv_facts[M-N]
    ans += tmp
    ans %= mod

ans *= mod_facts[M]*mod_inv_facts[M-N]
ans %= mod

print(ans)

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

def plmi(k):
    if k%2 == 0:
        return 1
    else:
        return -1


def mod_factorials(n, m):
    mod_factorials = [1] * (n + 1)
    for i in range(1, n+1):
        mod_factorials[i] = mod_factorials[i-1]*i
        mod_factorials[i] %= m
    return mod_factorials


def mod_inv_factorials(n, m):
    mod_invs = [1] * (n + 1)
    invs = [1] * (n + 1)
    for i in range(2, n + 1):
        invs[i] = m - invs[m % i] * (m//i) % m

    for i in range(1, n+1):
        mod_invs[i] = invs[i]*mod_invs[i-1]
        mod_invs[i] %= m
    return mod_invs

mod_facts = mod_factorials(M,mod)
mod_inv_facts = mod_inv_factorials(M,mod)

ans = 0

for i in range(0,N+1):
    tmp = 1
    pm = plmi(i)
    tmp *= mod_facts[N]*mod_facts[M-i]
    tmp %= mod
    tmp *= mod_inv_facts[i]*mod_inv_facts[N-i]*mod_inv_facts[M-N]
    tmp %= mod
    tmp *= pm
    ans += tmp
    ans %= mod

ans *= mod_facts[M]*mod_inv_facts[M-N]
ans %= mod

print(ans)

Solution (Python):

from collections import defaultdict,deque
from heapq import heappush, heappop
from bisect import bisect_left, bisect_right
import sys, random, itertools, math
sys.setrecursionlimit(10**5)
input = sys.stdin.readline
sqrt = math.sqrt
def LI(): return list(map(int, input().split()))
def LF(): return list(map(float, input().split()))
def LI_(): return list(map(lambda x: int(x)-1, input().split()))
def II(): return int(input())
def IF(): return float(input())
def S(): return input().rstrip()
def LS(): return S().split()
def IR(n): return [II() for _ in range(n)]
def LIR(n): return [LI() for _ in range(n)]
def FR(n): return [IF() for _ in range(n)]
def LFR(n): return [LI() for _ in range(n)]
def LIR_(n): return [LI_() for _ in range(n)]
def SR(n): return [S() for _ in range(n)]
def LSR(n): return [LS() for _ in range(n)]
mod = 1000000007
inf = float("INF")
mod = 1000000007
inf = float('INF')
fact = [i for i in range(5 * 10 ** 5 + 1)]
fact[0] = 1
for i in range(5 * 10 ** 5):
    fact[i + 1] *= fact[i]
    fact[i + 1] %= mod
invfact = fact[:]
invfact[-1] = pow(invfact[-1], mod - 2, mod)
for i in range(5 * 10 ** 5, 0, -1):
    invfact[i - 1] = i * invfact[i]
    invfact[i - 1] %= mod

def combination_mod(n, k, mod=mod):
    """ power_funcを用いて(nCk) mod p を求める """ 
    """ nCk = n!/((n-k)!k!)を使用 """
    from math import factorial
    if n < 0 or k < 0 or n < k: return 0
    if n == 0 or k == 0: return 1
    a = fact[n]
    b = invfact[k]
    c = invfact[n - k]
    return (a * b * c) % mod

def nPk(a, b):
    return (fact[a] * invfact[a - b]) % mod


def solve():
    n, m = LI()
    ans = 0
    a = 1
    for i in range(n + 1):
        ans += nPk(m, i) * nPk(m - i, n - i) * nPk(m - i, n - i) * a * combination_mod(n, i) % mod
        ans %= mod 
        a *= - 1
    print(ans)
            
    return



if __name__ == '__main__':
    solve()

Solution (Python):
N, M = map(int, input().split())
mod = int(1e9) + 7

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

for i in range( 2, N + M + 1 ):
    g1.append( ( g1[-1] * i ) % mod )
    inverse.append( ( -inverse[mod % i] * (mod//i) ) % mod )
    g2.append( (g2[-1] * inverse[-1]) % mod )


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 pmt(n, r, mod):
    if (r < 0 or r > n):
        return 0
    return cmb(n, r, mod) * g1[r] % mod

if __name__ == "__main__":
    ans = 0
    for i in range(N+1):
        if i % 2 == 1:
            ans -= cmb(N, i, mod) * pmt(M, i, mod) * (pmt(M - i, N - i, mod))**2 % mod
        else:
            ans += cmb(N, i, mod) * pmt(M, i, mod) * (pmt(M - i, N - i, mod))**2 % mod
    print(ans % mod)

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


factorial = [1]
for i in range(1, M + 1):
    factorial.append(factorial[-1] * i % MOD)
inv_factorial = [-1] * (M + 1)
inv_factorial[-1] = pow(factorial[-1], MOD - 2, MOD)
for i in reversed(range(M)):
    inv_factorial[i] = inv_factorial[i + 1] * (i + 1) % MOD

ans = 0
for k in range(N + 1):
    ans += (
        factorial[N]
        * inv_factorial[k]
        * inv_factorial[N - k]
        * (-1) ** (k % 2)
        * factorial[M]
        * inv_factorial[M - k]
        * pow(factorial[M - k] * inv_factorial[M - N], 2, MOD)
    )
    ans %= MOD

print(ans)

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

def mod_prepare(N):
    fac = [1, 1]
    finv = [1, 1]
    inv = [0, 1]
    for i in range(2, N + 1):
        fac.append(fac[i - 1] * i % mod)
        inv.append(mod - inv[mod % i] * (mod // i) % mod)
        finv.append(finv[i - 1] * inv[i] % mod)
    return fac, finv

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

fac, inv = mod_prepare(M)

ans = 0
sign = 1
for k in range(N + 1):
    ans += sign * fac[M - k] * inv[k] * inv[N - k]
    sign *= -1
    ans %= mod

print(ans * fac[N] * fac[M] * inv[M - N] * inv[M - N] % mod)

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



max = 5 * 10**5 + 100
fac, finv, inv = [0]*max, [0]*max, [0]*max

def comInit(max):
    fac[0] = fac[1] = 1
    finv[0] = finv[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

comInit(max)


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

cnt = [0] * (n+1)
tot = 0
for i in range(n,-1,-1):
    cnt[i] = com(n,i)*(fac[n-i])*com(m-i ,n-i)
    cnt[i] %= mod
    tot = cnt[i] - tot
    tot %= mod

ans = tot * fac[m] * finv[m-n]
ans %= mod
print(ans)

Solution (Python):
from collections import deque


def permutation(n, x, mod=10**9+7):
    
    
    
    tmp = 1
    for i in range(n, n-x, -1):
        tmp = (tmp * i) % mod
    return tmp

def power_func(a,n,p=10**9+7):
    bi = str(format(n,"b"))
    res = 1
    for i in range(len(bi)):
        res = (res*res) %p
        if bi[i] == "1":
            res = (res*a) %p
    return res


n,v=map(int,input().split())
m=10**9+7
perm1=deque([1])
perm2=deque([permutation(v,n)])
comb=deque([1])
ones=deque([1])
acc=((comb[-1]*perm1[-1]%m)*perm2[-1]**2)%m
for k in range(n):
  
  comb.append((comb.pop()*(n-k)*power_func(k+1,m-2))%m)
  
  perm1.append((perm1.pop()*(v-k))%m)
  
  perm2.append((perm2.pop()*power_func(v-k,m-2))%m)
  
  ones.append(ones.pop()*(-1))
  acc+=ones[-1]*((comb[-1]*perm1[-1]%m)*perm2[-1]**2)%m
print(acc%m)

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

a = 10 ** 7 + 10 ** 6
fact = [0] * a
fact[0] = 1
for i in range(1, a):
    fact[i] = fact[i-1] * i % MOD

inv_fact = [0] * a
inv_fact[a-1] = pow(fact[a-1], MOD-2, MOD)
for i in range(a-1, 0, -1):
    inv_fact[i-1] = inv_fact[i] * i % MOD


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


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


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

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

dp = [0] * n
mod = 10**9 + 7


if n == 1:
    print((m-n) * m % mod)
elif n == 2:
    ans = (m - n)**2 + (m - n + 1)
    for i in range(2):
        ans *= (m - i)
        ans %= mod
    print(ans)
else:
    dp[0] = m - n
    dp[1] = (m - n)**2 + (m - n + 1)


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

    ans = dp[-1]

    for i in range(n):
        ans *= (m - i)
        ans %= mod

    print(ans)

Solution (Python):
class nCrMod():
    def __init__(self, mod):
        self.mod = mod
        self.fac = [1, 1]
        self.finv = [1, 1]
        self.inv = [0, 1]

    def prep(self, n):
        mod = self.mod
        f, fi = self.fac[-1], self.finv[-1]
        for i in range(len(self.fac), n + 1):
            fn = f * i % mod
            v = -self.inv[mod % i] * (mod // i) % mod
            fin = fi * v % mod
            f, fi = fn, fin
            self.fac.append(f)
            self.finv.append(fi)
            self.inv.append(v)

    def __call__(self, n, r):
        if len(self.fac) <= n:
            self.prep(n)
        return self.fac[n] * self.finv[r] * self.finv[n - r] % self.mod


def main():
    N, M = map(int, input().split())
    mod = 10**9 + 7
    nCr = nCrMod(mod)
    nPr = lambda n, r: nCr(n, r) * nCr.fac[r] % mod
    r = nPr(M, N) ** 2 % mod
    for i in range(1, N + 1):
        r = (r + (-1)**i * nCr(N, i) * nPr(M, i) * (nPr(M - i, N - i) ** 2)) % mod
    return r % mod

print(main())

Solution (Python):
n,m=map(int,input().split())
mod=10**9+7
framod=[1]
def framod_calc(n, mod, a=1):
    for i in range(1,n+1):
        a=a * i % mod
        framod.append(a)

def permmod(n, k, mod):
    a=framod[n]
    c=framod[n-k]
    return (a * pow(c, mod-2, mod)) % mod

def combmod(n, k, mod):
    a=framod[n]
    b=framod[k]
    c=framod[n-k]
    return (a * pow(b, mod-2, mod) * pow(c, mod-2, mod)) % mod
framod_calc(m, mod)

a=combmod(m,n,mod)*framod[n]
a%=mod
b=0
for k in range(n+1):
    b+=((((-1)**k)*combmod(n,k,mod)*combmod(m-k,n-k,mod))%mod*framod[n-k])%mod
    b%=mod
print(a*b%mod)