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

Solution (Python):
class Factorial:
    def __init__(self, max_fact, mod):
        
        
        f = [1] * 2
        for idx in range(2, max_fact + 1):
            f += [f[idx - 1] * 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):

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

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


    ｍPn = fact[m] * pow(fact[m-n],MOD-2,MOD) % MOD

    a = pow(mPn,2,MOD)
    b = 0
    for i in range(1,n+1):
        nCr = fact[n] * fact_inv[i] % MOD * fact_inv[n-i] % MOD
        b += (-1)**(i%2) * nCr * fact[m-i] % MOD * fact_inv[m-n] % MOD
        b %= MOD
    ans = a + b * mPn
    print(ans%MOD)


main()

Solution (Python):
class Factorial:
    def __init__(self, max_fact, mod):
        
        
        f = [1] * 2
        for idx in range(2, max_fact + 1):
            f += [f[idx - 1] * 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 = 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):
max_fact = 5 * 10**5
mod = 10**9 + 7

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]

def factorial(n):
    return f[n]

def factorial_inverse(n):
    return fi[n]

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

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

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


comb = combination
perm = 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 Combination():
    def __init__(self, n, mod=10**9+7):
        self.mod = mod
        self.fac = [1]*(n+1)
        for i in range(1,n+1):
            self.fac[i] = self.fac[i-1] * i % self.mod
        self.invfac = [1]*(n+1)
        self.invfac[n] = pow(self.fac[n], self.mod - 2, self.mod)
        for i in range(n-1, 0, -1):
            self.invfac[i] = self.invfac[i+1] * (i+1) % self.mod

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

    def permutation(self, n, r):
        return self.factorial(n) * self.invfactorial(n-r) % self.mod

    def factorial(self, i):
        return self.fac[i]

    def invfactorial(self, i):
        return self.invfac[i]

def main():
    import sys
    def input(): return sys.stdin.readline().rstrip()
    n, m = map(int, input().split())
    c = Combination(m)
    mod = 10**9 + 7
    ans = 0
    for i in range(n+1):
        tmp = c.combination(n,i)*pow(-1,i%2)*c.permutation(m, i)%mod
        tmp *= c.permutation(m-i, n-i)**2
        tmp %= mod
        ans += tmp
        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):
N, M = map(int, input().split())
ans = 0
mod = 10**9+7
def modinv(a, mod=10**9+7):
    return pow(a, mod-2, mod)

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

def combi(n,r,mod=10**9+7):
    return( ( fac[n] * modinv(fac[r]) )%mod * modinv(fac[n-r]) % mod)

for k in range(N+1):
    ans += (fac[M-k]*modinv(fac[M-N])) % mod * combi(N, k) * (-1)**k
    ans %= mod

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

Solution (Python):
N, M = map(int, input().split())
ans = 0
mod = 10**9+7
def modinv(a, mod=10**9+7):
    return pow(a, mod-2, mod)

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

def combi(n,r,mod=10**9+7):
    return( ( fac[n] * modinv(fac[r]) )%mod * modinv(fac[n-r]) % mod)

for k in range(N+1):
    ans += (fac[M-k]*modinv(fac[M-N])) % mod * combi(N, k) * (-1)**k

print(ans*fac[M]*modinv(fac[M-N]) % 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):
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 prepare(n):
    global mod
    modFacts = [0] * (n + 1)
    modFacts[0] = 1
    for i in range(n):
        modFacts[i + 1] = (modFacts[i] * (i + 1)) % mod

    invs = [1] * (n + 1)
    invs[n] = pow(modFacts[n], mod - 2, mod)
    for i in range(n, 1, -1):
        invs[i - 1] = (invs[i] * i) % mod

    return modFacts, invs

n,m = map(int,input().split())
mod = 10**9+7
modFacts,invs = prepare(m)
ans = 0
nci = 0
for i in range(n+1):
    nci = (modFacts[n]*invs[n-i]*invs[i]) % mod
    mpi = (modFacts[m]*invs[m-i]) % mod
    mipni = (modFacts[m-i]*invs[m-n]) % mod
    ans += ((nci * mpi * mipni**2) % mod )* (-1)**i
    ans %= mod
    
print(ans)

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


MOD = 10**9+7
size = 500001
frac = [1]*size
finv = [1]*size
for i in range(size-1):
    frac[i+1] = (frac[i]*(i+1))%MOD
finv[-1] = pow(frac[-1],MOD-2,MOD)
for i in range(size-2,-1,-1):
    finv[i] = (finv[i+1]*(i+1))%MOD
def comb(n,r):
    if n < r or r < 0:
    	return 0
    elif n == 0 or r == 0 or n == r:
        return 1
    else:
        return (frac[n]*finv[n-r]*finv[r])%MOD
def nPk(n,k):
    return frac[n]*finv[n-k]%MOD
ans = 0
for i in range(n+1):
    res = comb(n,i) * nPk(m,i) * (nPk(m-i,n-i)**2)
    res %=MOD
    if i%2: res*=-1
    ans += res
    ans %= MOD

print(ans%MOD)

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

f_list = [1] * list_size
f_r_list = [1] * list_size

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

f_r_list[-1] = pow(f_list[-1], MOD - 2, MOD)

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

def comb(n, r):
	if n < r or r < 0:
		return 0
	elif n == 0 or r == 0 or n == r:
		return 1
	else:
		return (f_list[n] * f_r_list[n-r] * f_r_list[r]) % MOD 

n, m = map(int, input().split())
ans = 0
for i in range(n+1):
	res = comb(n, i) * f_list[m-i] * f_r_list[m-n]
	res %= MOD
	if i%2:
		res *= -1
	ans += res
	ans %= MOD

ans *= (f_list[m] * f_r_list[m-n])
print(ans%MOD)

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


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


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


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


ans = pow(nCr(M, N) * factorial[N], 2, MOD)
for r in range(1, N + 1):
    p = nCr(M, r) * nPr(N, r) * pow(nPr(M - r, N - r), 2, MOD) % MOD
    p *= ((-1) ** ((r % 2) ^ 1))
    ans -= p
    ans %= MOD

print(ans)

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


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


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


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


ans = pow(nCr(M, N) * factorial[N], 2, MOD)
for r in range(1, N + 1):
    p = nCr(M, r) * pow(nCr(M - r, N - r), 2, MOD) * nCr(N, r) * factorial[r] * factorial[N - r] * factorial[N - r] % MOD
    p *= ((-1) ** ((r % 2) ^ 1))
    ans -= p
    ans %= MOD

print(ans)

Solution (Python):
'''
包除原理 |A∪B∪C|=|A|+|B|+|C|−|A∩B|−|A∩C|−|B∩C|+|A∩B∩C|
悪魔の証明のような否定の個数を求めるのは難しい
具体的には、A[i], B[i], 1 <= i <= N のとき
A[i] != B[i] となる組み合わせの個数は求めにくい。
しかし、A[i] == B[i] となる項が少なくともk個ある組み合わせは求めやすい(nCk * (n-k)!)
そこで包除原理を使えば、
Σ[k=0,n](nCk * (n-k)! * (-1)^k) で求まる (k=0がここでは全事象を表す)
'''

mod = pow(10,9)+7
n,m = map(int,input().split()) 
inv = [0,1]
fact = [1,1]
ifact = [1,1]
for i in range(2,m+1):
    inv.append(inv[mod%i]*(mod-int(mod/i))%mod)
    fact.append(i*fact[-1]%mod)
    ifact.append(ifact[-1]*inv[i]%mod)

def combi(n,r):
    r = min(r, n - r)
    return fact[n]*ifact[r]*ifact[n-r]%mod
def permi(n,r):
    return fact[n]*ifact[n-r]%mod

ans = 0
sign = 1
for k in range(n+1):
    ans += combi(n,k) * permi(m-k,n-k) * sign % mod
    ans %= mod
    sign *= -1
ans = ans * permi(m,n) % mod
print(ans)

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):
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):
    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):
MOD = 10**9 + 7
fact = [0] * 1100000
factinv = [0] * 1100000
fact[0] = 1
for i in range(1,1100000):
    fact[i] = fact[i-1] * i% MOD
factinv[-1] = pow(fact[-1], MOD-2, MOD)
for i in range(1100000 - 1, 0, -1):
    factinv[i-1] = factinv[i] * i % MOD

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

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

ans = 0
n, m = map(int,input().split())
for i in range(n+1):
    temp = nCk(n,i)
    temp = (nCk(n,i) * nPk(m,i) * nPk(m-i, n-i)**2)% MOD
    if i&1:
        temp *= -1
    ans += temp
    ans %= MOD
print(ans)

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


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


ans = 0
t = 0
cnt = 0
a = 1
for j in range(n + 1):
    if j % 2 == 0:
        ans += comb(m - j, n - j) * fac[n - j] * comb(n, j)
    else:
        ans -= comb(m - j, n - j) * fac[n - j] * comb(n, j)
print(ans * comb(m, n) * fac[n] % 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<0 or r>n:
    return 0
  return (g1[n]*g2[n-r])%mod

g1=[1,1]
g2=[1,1]
inv=[0,1]
for i in range(2,N+M+10):
  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=(P+pow(-1,i&1)*cmb(N,i)*perm(M,i)*perm(M-i,N-i)**2)%mod
print(P)

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

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


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


def nCr(n, r):
    if r < 0 or n < r:
        return 0
    return fact[n]*invfact[r]*invfact[n-r]


def nPr(n, r):
    return fact[n]*invfact[n-r]


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

Solution (Python):
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
mod = 10**9 + 7

fact = [1, 1]
factinv = [1, 1]
inv = [0,1]

n , m = map(int,input().split())
for i in range(2,m+1):
    fact.append((fact[-1] * i) % mod)
    inv.append((inv[mod % i] * (mod - (mod // i))) % mod)
    factinv.append((factinv[-1] * inv[-1]) % mod)

ans = 0

for i in range(0,n+1):
    ans += cmb(n,i,mod)*((-1)**i)*fact[m]*factinv[m-i]*(fact[m-i]*factinv[m-n])**2
    ans %= mod
print(ans)

Solution (Python):
class PrimeTool():
    def __init__(self, max=5 * 10 ** 5 + 500, mod=10 ** 9 + 7):
        self.mod = mod
        self.modm2 = self.mod - 2
        self.max = max
        self.fac = [1] * max
        self.inv = [1] * max
        if self.max > 1:
            for i in range(2, self.max):
                self.fac[i] = (self.fac[i - 1] * i) % self.mod
                self.inv[i] = pow(self.fac[i], self.modm2, self.mod)

    def perm(self, n, r):
        if r == 0:
            return 1
        return self.fac[n] * self.inv[n - r] % self.mod

    def cmb(self, n, r):
        r = max(r, n - r)
        if n == r:
            return 1
        return self.fac[n] * self.inv[n - r] * self.inv[r] % self.mod


def resolve():
    N, M = map(int, input().split(' '))
    res = 0
    tool = PrimeTool()

    for k in range(N + 1):
        res += (-1) ** k * tool.cmb(N, k) * tool.perm(M, k) * tool.perm(M - k, N - k) ** 2
        res %= tool.mod

    print(res)


resolve()

Solution (Python):
n,m = map(int,input().split())
MOD = 10**9 + 7
FAC = [1]
INV = [1]
for i in range(1,m+1):
  FAC.append((FAC[i-1]*i) % MOD)
  INV.append(pow(FAC[-1],MOD-2,MOD))

def nCr(n,r):
  return (FAC[n]*INV[n-r]*INV[r])%MOD

def nPr(n,r):
  return (FAC[n]*INV[n-r])%MOD














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

Solution (Python):
n,m=map(int,input().split())
i=b=0
a=1
while i<n:a,b,i=((m-i)*((m-n+i)*a+i*(m-i+1)*b))%(10**9+7),a,i+1
print(a)

Solution (Python):










import sys

def input(): return sys.stdin.readline()
def iIn(): return int(input())
def iInM(): return map(int, input().split())
def iInM1(): return map(int1, input().split())
def iInLH(): return list(map(int, input().split()))
def iInLH1(): return list(map(int1, input().split()))
def iInLV(n): return [iIn() for _ in range(n)]
def iInLV1(n): return [iIn()-1 for _ in range(n)]
def iInLD(n): return [iInLH() for _ in range(n)]
def iInLD1(n): return [iInLH1() for _ in range(n)]
def sInLH(): return list(input().split())
def sInLV(n): return [input().rstrip("\n") for _ in range(n)]
def sInLD(n): return [sInLH() for _ in range(n)]

def OutH(lst): print(*lst)
def OutV(lst): print(*lst, sep="\n")

sys.setrecursionlimit(10 ** 6)

int1 = lambda x: int(x) - 1

INF = int(1e9)
MOD = 1000000007
dx = (-1, 0, 1, 0)
dy = (0, -1, 0, 1)






class Factorial:

    def __init__(self, N):
        self.fact = [0] * (N+1)
        self.fact_inv = [0] * (N+1)
        self.calc_factorial(N)
        self.calc_inv_factorial(N)

    def calc_factorial(self, N):
        self.fact[0] = 1
        for i in range(N):
            self.fact[i+1] = self.fact[i] * (i+1) % MOD

    def calc_inv_factorial(self, N):
        self.fact_inv[N] = pow(self.fact[N], MOD-2, MOD)
        for i in range(N, 0, -1):
            self.fact_inv[i-1] = self.fact_inv[i] * i % MOD

    def comb(self, n, r):
        if (n < 0 or r < 0 or n < r): return 0
        res = self.fact[n]
        res = res * self.fact_inv[r] % MOD
        res = res * self.fact_inv[n-r] % MOD
        return res

    def perm(self, n, r):
        if (n < 0 or r < 0 or n < r): return 0
        res = self.fact[n]
        res = res * self.fact_inv[n-r] % MOD
        return res



N, M = iInM()

fac = Factorial(500000)
ans = 0

for i in range(N+1):
    res = (-1)**(i) * fac.comb(N, i)
    res = res * fac.perm(M, i) % MOD
    res = res * fac.perm(M-i, N-i) % MOD
    res = res * fac.perm(M-i, N-i) % MOD
    ans = (ans + res) % MOD

print(ans)

Solution (Python):
n,m=map(int,input().split())
i=b=0
a=1
while i<n:a,b=((m-i)*((m-n+i)*a+i*(m-i+1)*b))%(10**9+7),a;i+=1
print(a)

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

kai = [1 for _ in range(N+M+1)]
for i in range(2,N+M+1):
    kai[i] = (kai[i-1]*i)%mod

ap = 0
for i in range(N+1):
    sp = (kai[M]*pow(kai[M-i], mod-2, mod))%mod
    restp = (kai[M-i]*pow(kai[M-N], mod-2, mod))%mod
    chose = ((kai[N]*pow(kai[N-i], mod-2, mod))%mod*pow(kai[i], mod-2,mod))%mod
    sign = (-1)**i
    ap = (ap + (sign*sp)%mod*(chose*restp**2)%mod)%mod
print(ap%mod)

Solution (Python):
n,m=map(int,input().split())
i=b=0
a=1
while i<n:c=((m-i)*((m-n+i)*a+(m-i+1)*i*b))%(10**9+7);i+=1;a,b=c,a
print(a)

Solution (Python):
n,m=map(int,input().split())
i,a,b=0,1,0
while i<n:c=((m-i)*((m-n+i)*a+i*b*(m-i+1)))%(10**9+7);i+=1;a,b=c,a
print(a)

Solution (Python):
MOD = 10 ** 9 + 7
def init(n=5*10**5):
    fact = [1] * (n + 1)
    for i in range(2, n+1):
        fact[i] = i * fact[i-1] % MOD
    inv = [1] * (n + 1)
    inv[n] = pow(fact[n], MOD-2, MOD)
    for i in range(n, 0, -1):
        inv[i-1] = i * inv[i] % MOD
    return fact, inv

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

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

N, M = map(int, input().split())
fact, inv = init()
res = 0
sgn = 1
for i in range(N+1):
    res += sgn * nPr(M - i, N - i) * nCr(N, i)
    res %= MOD
    sgn *= -1
print(nPr(M, N) * res % MOD)

Solution (Python):
n,m=map(int,input().split())
i,d=0,[1]*(n+1)
while i<n:d[i+1]=((m-i)*((m-n+i)*d[i]+i*d[i-1]*(m-i+1)))%(10**9+7);i+=1
print(d[-1])

Solution (Python):


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

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 pmt(n, r, p):
    if (r < 0) or (n < r):
        return 0
    
    return fact[n] * factinv[n - r] % p
  
p = 10 ** 9 + 7
N = m + 1
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)



ans = 0
for i in range(n + 1):
  a = (-1) ** i
  b = cmb(n, i, p)
  c = pmt(m, i, p)
  d = pmt(m - i, n - i, p)
  
  ans += (a * b * c * (d ** 2)) % p
  
print(ans % p)

Solution (Python):
n,m=map(int,input().split())
i,d=0,[1]*(n+1)
while i<n:
  d[i+1]=((m-i)*((m-n+i)*d[i]+i*d[i-1]*(m-i+1)))%(10**9+7);i+=1
print(d[-1])

Solution (Python):
n,m=map(int,input().split())
mod=10**9+7
fact=[1]*(m+1)
inv=[1]*(m+1)
for i in range(2,m+1):
    fact[i]=i*fact[i-1]%mod
inv[-1]=pow(fact[-1],mod-2,mod)
for i in range(m,1,-1):
    inv[i-1]=inv[i]*i%mod
def comb(x,y):return fact[x]*inv[y]%mod*inv[x-y]%mod if x>=y>=0 else 0
def p(x,y):return fact[x]*inv[x-y]%mod if x>=y>=0 else 0
ans=p(m,n)**2%mod
a=-1
for i in range(1,n+1):
    ans=(ans+a*(p(m,i)*comb(n,i)*p(m-i,n-i)**2%mod))%mod
    a*=-1
print(ans)

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

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

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

Solution (Python):
"""
nCk := f[n] * g[n - k] * g[k](mod)
"""
def pre_c(n, mod):
    f = [0] * (n + 1)
    g = [0] * (n + 1)
    for i in range(n + 1):
        if i == 0:
            f[i] = 1
            g[i] = 1
        else:
            f[i] = (f[i - 1] * i) % mod
            g[i] = pow(f[i], mod - 2, mod)

    return f, g

n, m = map(int,input().split())
mod = 10 ** 9 + 7
f, g = pre_c(m, mod)
res = 0
for i in range(n + 1):
    C = f[n] * g[i] * g[n - i]
    P1 = f[m] * g[m - i]
    P2 = f[m - i] * g[m - n]
    temp = C * P1 * (P2 ** 2) * pow(-1, i)
    res = (res + temp%mod) %mod
print(res%mod)

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

Solution (Python):
def prepare(n, MOD):
 
    
    f = 1
    fn=[1]*(n+1)
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        fn[m]=f
     
    
    inv = pow(f, MOD - 2, MOD)
    
    invs = [1] * (n + 1)
    invs[n] = inv
    for m in range(n, 1, -1):
        inv *= m
        inv %= MOD
        invs[m - 1] = inv
     
    return fn, invs

def main():
    n,m=map(int,input().split())
    MOD=10**9+7
    fn,invs=prepare(m,MOD)
    ans=0
    for i in range(n+1):
        x=(fn[m-i]*invs[m-n])%MOD
        x=(x*fn[n]*invs[i]*invs[n-i])%MOD
        x*=(-1)**i
        ans+=x
    ans=(ans*fn[m]*invs[m-n])%MOD
    print(ans)

if __name__ == '__main__':
    main()

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

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

Solution (Python):
n,m=map(int,input().split())
M,F=10**9+7,[1]
for i in range(m): F+=[-~i*F[i]%M]
c=lambda n,r: F[n]*pow(F[r]*F[n-r],M-2,M)
p=lambda n,r: F[n]*pow(F[n-r],M-2,M)
print(sum((-1)**k*c(n,k)*p(m,k)*p(m-k,n-k)**2 for k in range(n+1))%M)

Solution (Python):
M=10**9+7
n,m=map(int,input().split())
F=[1]
for i in range(1,m+1): F+=[i*F[-1]%M]
c=lambda n,r: F[n]*pow(F[r]*F[n-r],M-2,M)%M
p=lambda n,r: F[n]*pow(F[n-r],M-2,M)%M
a=0
for k in range(n+1):
  a+=(-1)**k*c(n,k)*p(m,k)*p(m-k,n-k)**2%M
print(a%M)

Solution (Python):
import sys
import math

sys.setrecursionlimit(10 ** 8)

input = sys.stdin.readline


def main():
    N, M = [int(x) for x in input().split()]

    MOD = 10 ** 9 + 7

    fact = [1, 1]
    factinv = [1, 1]
    inv = [0, 1]

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

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

    X = fact[M] * factinv[M - N]

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

        ans %= MOD

    print(ans % MOD)


if __name__ == '__main__':
    main()

Solution (Python):
mod = pow(10,9)+7
n,m = map(int,input().split())
ans = 0
inv = [0,1]
fact = [1,1]
ifact = [1,1]
for i in range(2,m+1):
    inv.append(inv[mod%i]*(mod-int(mod/i))%mod)
    fact.append(i*fact[-1]%mod)
    ifact.append(ifact[-1]*inv[i]%mod)

def combi(a,b):
    return fact[a]*ifact[b]*ifact[a-b]%mod
def permi(a,b):
    return fact[a]*ifact[a-b]%mod

for i in range(n+1):
    now = combi(n,i)
    now *= permi(m-i,n-i)
    now %= mod
    ans += now*(-1)**(i&1)
    ans %= mod
ans *= permi(m,n)
ans %= mod
print(ans)

Solution (Python):
n,m=map(int,input().split())
M,F=10**9+7,[1]
for i in range(m): F+=[-~i*F[i]%M]
print(sum((-1)**k*F[m-k]*pow(F[k]*F[n-k]*F[m-n]**2,-1,M) for k in range(n+1))*F[n]*F[m]%M)

Solution (Python):
mod = pow(10,9)+7
n,m = map(int,input().split())
ans = 0
inv = [0,1]
fact = [1,1]
ifact = [1,1]
for i in range(2,m+1):
    inv.append(inv[mod%i]*(mod-int(mod/i))%mod)
    fact.append(i*fact[-1]%mod)
    ifact.append(ifact[-1]*inv[i]%mod)

def combi(a,b):
    return fact[a]*ifact[b]*ifact[a-b]%mod
def permi(a,b):
    return fact[a]*ifact[a-b]%mod

for i in range(n+1):
    now = combi(n,i)
    now *= permi(m-i,n-i)
    now %= mod
    ans += now*(-1)**(i&1)
    ans %= mod
ans *= permi(m,n)
ans %= mod
print(ans)

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