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
input=sys.stdin.readline
MOD=10**9+7

fac=[1]*(5*10**5+1)
for i in range(1,5*10**5+1):
    fac[i]=(fac[i-1]*i)%MOD
    
def inv(x):
    return pow(x,MOD-2,MOD)

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

def main():
    N,M=map(int,input().split())
    hojo=0
    for i in range(1,N+1):
        hojo+=(-1)**(i)*nCr(N,i)*nPr(M,i)*nPr(M-i,N-i)**2
        hojo%=MOD
    print((nPr(M,N)**2+hojo)%MOD)
    
    
    
if __name__=="__main__":
    main()

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

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

fac = [1]*(m+1)
fac_inv =  [0]*(m+1)

for i in range(m):
    fac[i+1] = (fac[i] * (i+1)) %MOD
fac_inv[-1] = pow(fac[-1],MOD-2,MOD)
for i in range(m,0,-1):
    fac_inv[i-1] = (fac_inv[i] * i) % MOD
    
def comb(n,k):
    if k<0 or k>n:
        return 0
    x = (fac[n] * (fac_inv[k] * fac_inv[n-k])%MOD)%MOD
    return x


ans = 0
for i in range(n+1):
    ans += (-1)**(i%2) * comb(n,i) * (fac[m-i] * fac_inv[m-n] %MOD) % MOD
    ans %= MOD

ans *= comb(m,n)*fac[n] %MOD
ans %= MOD
    
print(ans)

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

def power(x, y):
    if   y == 0     : return 1
    elif y == 1     : return x % mod
    elif y % 2 == 0 : return power(x, y//2)**2 % mod
    else            : return power(x, y//2)**2 * x % mod

inverseFactorial = [0] * (m)
inverseFactorial[m-1] = power(factorial[m-1], mod-2)
for i in range(m-2, -1, -1):
    inverseFactorial[i] = inverseFactorial[i+1] * (i+1) % mod

def combi(a,b):
    return factorial[a] * inverseFactorial[b] * inverseFactorial[a-b] % mod

m-=1


sin=-1

ans=0
for i in range(n+1):
    sin*=-1
    nck=combi(n,i)
    npk=(combi(m-i,n-i)*factorial[n-i]) % mod
    tasu=(sin*nck*npk)%mod
    ans+=tasu
    ans%=mod

kake=(combi(m,n)*factorial[n]) % mod

print((ans*kake)%mod)

Solution (Python):
fact = [1]*500001

mod = 10**9 + 7

for i in range(1, 500001):
	fact[i] = fact[i-1]*i
	fact[i] %= mod


def nCr(n, r):
	global fact, mod
	if n < r:
		return 0
	return fact[n] * pow(fact[r]*fact[n-r], mod-2, mod)


n, m = map(int, input().split())
mult = fact[m] * pow(fact[m-n], mod-2, mod) % mod

a = 0
for i in range(0, n+1):
	add = nCr(n, i) * nCr(m-i, n-i) * fact[n-i] % mod
	if i%2 == 0:
		a += add
	else:
		a -= add

print(mult*a % mod)

Solution (Python):
fact = [1]
inv_fact = [1]
MOD = 10**9+7
temp = 1
for i in range(1,5*10**5+2):
    temp *= i
    temp %= MOD
    fact += [temp]
    inv_fact += [pow(temp,MOD-2,MOD)]

def bino(a,b):
    MOD = 10**9+7
    up = fact[a]
    down = (inv_fact[a-b]*inv_fact[b])%MOD
    return (up*down)%MOD

def find(N,M):
    MOD = 10**9+7
    temp1 = fact[M]*pow(inv_fact[M-N],2,MOD)
    temp2 = fact[M]
    for k in range(1,N+1):
        temp2 += (-1)**k*bino(N,k)*fact[M-k]
        temp2 %= MOD
    
    ans = (temp2*temp1)%MOD
    return ans

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

print(find(N,M))

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

class Factorial:
  def __init__(self, n, mod):
    self.f = [1]
    for i in range(1, n + 1):
      self.f.append(self.f[-1] * i % mod)
    self.i = [pow(self.f[-1], mod - 2, mod)]
    for i in range(1, n + 1)[: : -1]:
      self.i.append(self.i[-1] * i % mod)
    self.i.reverse()
  def factorial(self, i):
    return self.f[i]
  def ifactorial(self, i):
    return self.i[i]
  def combi(self, n, k):
    return self.f[n] * self.i[n - k] % mod * self.i[k] % mod
  def permi(self, n, k):
    return self.f[n] * self.i[n - k] % mod

f = Factorial(max(N, M), mod)

res = f.permi(M, N) ** 2 % mod
for x in range(1, N + 1):
  t = f.combi(N, x) * f.permi(M - x, N - x) % mod
  if x % 2: res -= t * f.permi(M, N)
  else: res += t * f.permi(M, N)
  res %= mod
print(res)

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

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
N = 510000  
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)

targ = m-n
"""
m = 10, n = 8とする.
a = 1, 2, 3, 4, 5, 6, 7, 8
b = ?????????????????????
dp[i] = i番目まで選んだとき差がm-n
"""

dp = [0 for _ in range(n+1)]
dp[0] = 1
for i in range(1, n+1):
	if i == 1:
		dp[1] = targ
	elif i == 2:
		dp[i] = dp[0] + dp[1]*(targ+1)
	else:
		dp[i] = (dp[i-2]*(i-1) + dp[i-1]*(targ+i-1)) % mod

print((dp[n]*cmb(m, n, mod)*fact[n]) % mod)

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

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

maxn = 5*(10**5) 
maxn+=10
fac = [0 for i in range(maxn)]
fac[0] = 1
for i in range(1,maxn):
    fac[i] = (fac[i-1]*i)%mod

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

def func(n,r):
    ans = fac[n]
    ans = (ans*inv(fac[n-r]))%mod
    ans = (ans*inv(fac[r]))%mod
    return ans







ans = 0
for i in range(n+1):
    if i%2==0:
        t = 1
    else:
        t = -1
    p = func(n,i)
    q = (fac[m]*inv(fac[m-i]))%mod
    r = func(m-i,n-i)
    s = fac[n-i]
    temp = (r*s)%mod
    temp = pow(temp,2,mod)

    h = (((p*q)%mod)*temp)%mod
    h = (h*t)
    ans = (ans+h)%mod
print(ans)

Solution (Python):
n,m=map(int,input().split())
mod=10**9+7
def c(a,b):
	r=fac[b]*fac[a-b]%mod
	return fac[a]*pow(r,mod-2,mod)%mod
def p(a,b):
	return c(a,b)*fac[b]%mod
fac=[1]
for i in range(m+2):
	fac.append(fac[-1]*(i+1)%mod)
ans=p(m,n)
for i in range(1,n+1):
	r=c(n,i)*p(m-i,n-i)%mod
	if i%2==1:
		ans-=r
	else:
		ans+=r
	ans%=mod
print(ans*p(m,n)%mod)

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

g2 = [0,1]

for i in range(2,m + 1):
    g2.append((-g2[mod % i] * (mod//i))% mod)
    g1.append(g1[-1] * g2[-1] % mod)
g3 = [1,m - n + 1]
for i in range(n - 1):
    g3.append((g3[-1] * (m - n + i + 2))%mod)
x = 1
for i in range(2,n + 1):
    x = (x*i)%mod
k = 0
for i in range(1,n + 1):
    k = (k + (-1)**i * x * g1[i] * g1[n - i] * g3[n - i])%mod
k = (k + g3[n])%mod
k = (k * g3[n])%mod
print(k)

Solution (Python):
N, M = map(int, input().split())
nn = 500500
P = 10 ** 9 + 7
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
PP = lambda a, b: fa[a] * fainv[a-b] % P if 0 <= b <= a else 0

ans = 0
for i in range(N + 1):
    ans = (ans + C(N, i) * PP(M - i, N - i) * (-1) ** i) % P
    
print(ans * PP(M, N) % P)

Solution (Python):
n, m = map(int, input().split())
p = 10**9 + 7
a = 1
fct = [a]
for i in range(1, m + 1):
    a = (a * i) % p
    fct.append(a)
b = pow(fct[-1], p - 2, p)
inv = [b]
for i in range(1, m + 1)[::-1]:
    b = (b * i) % p
    inv.append(b)
inv.reverse()
ans = (fct[m] * inv[m - n]) % p

tmp = 0
for i in range(1, n + 1):
    tmp += ((-1)**(i - 1) * fct[n] * inv[n - i] *
            inv[i] * fct[m - i] * inv[m - n]) % p
    tmp %= p
tmp = (fct[m] * inv[m - n] - tmp) % p
print((ans * tmp) % p)

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



class Data():
    def __init__(self):
        self.power = 1
        self.rev = 1


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

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


c = Combi(10**6, mod)
ans = (c.lists[M].power*c.lists[M-N].rev)**2 % mod
for k in range(1, N+1):
    s = (-1)**((k-1) % 2)
    t = c.combi(N, k)
    u=((c.lists[M-k].power*c.lists[M-N].rev)**2) % mod
    v=(c.lists[M].power*c.lists[M-k].rev) % mod
    ans -= s*t*u*v
print(ans % mod)

Solution (Python):
MOD = 10 ** 9 + 7
N,M = map(int,input().split())
ans1 = 1
for i in range(N):
    ans1 *= M - i
    ans1 %= MOD
ans2 = [1,M-N]
for i in range(2,N+1):
    ans2.append((ans2[-2] * (i -1) + ans2[-1] * (M - N + i - 1))% MOD)
print(ans1 * ans2[N] % MOD)

Solution (Python):
mod=10**9+7
fact=[1]
for i in range(1,5*10**5+1):
  fact.append((fact[-1]*i)%mod)
revfact=[1]
for i in range(1,5*10**5+1):
  revfact.append(pow(fact[i],mod-2,mod))
n,m=map(int,input().split())
patternA=(fact[m]*revfact[n]*revfact[m-n]*fact[n])%mod
patternB=patternA
for i in range(1,n+1):
  tmp=(fact[n]*revfact[i]*revfact[n-i])*(fact[m-i]*revfact[n-i]*revfact[m-n])*fact[n-i]
  if i%2==0:
    patternB+=tmp
    patternB%=mod
  else:
    patternB-=tmp
    patternB%=mod
print((patternA*patternB)%mod)

Solution (Python):
import sys
sys.setrecursionlimit(10000000)
MOD = 10 ** 9 + 7
INF = 10 ** 15


def main():
    N,M = map(int,input().split())
    MAXN = M + 2
    fact = [1]
    for i in range(1,MAXN + 1):
        fact.append(fact[-1]*i%MOD)
    inv_fact = [-1] * (MAXN + 1)
    inv_fact[-1] = pow(fact[-1],MOD - 2,MOD)
    for i in range(MAXN - 1,-1,-1):
        inv_fact[i] = inv_fact[i + 1]*(i + 1)%MOD
    
    nck = lambda N,K: 0 if K > N or K < 0 else fact[N]*inv_fact[N - K]*inv_fact[K]%MOD
    npk = lambda N,K: 0 if  K > N or K< 0 else fact[N]*inv_fact[N - K]%MOD

    ans = npk(M,N)*npk(M,N)%MOD
    cum = 0
    for i in range(N,0,-1):
        ret = nck(N,i)*npk(M,i)%MOD*npk(M-i,N-i)*npk(M-i,N-i)%MOD
        ret -= 2*cum
        ans -= ret
        ans %= MOD
        cum += ret
        cum %= MOD
    print(ans)
if __name__ == '__main__':
    main()

Solution (Python):
import sys
readline = sys.stdin.readline
 
MOD = 10**9+7
def frac(limit):
    frac = [1]*limit
    for i in range(2,limit):
        frac[i] = i * frac[i-1]%MOD
    fraci = [None]*limit
    fraci[-1] = pow(frac[-1], MOD -2, MOD)
    for i in range(-2, -limit-1, -1):
        fraci[i] = fraci[i+1] * (limit + i + 1) % MOD
    return frac, fraci
frac, fraci = frac(1341398)
def comb(a, b):
    if not a >= b >= 0:
        return 0
    return frac[a]*fraci[b]*fraci[a-b]%MOD
N, M = map(int, input().split())
X = 0
for i in range(N+1):
    t = comb(N, i)*frac[M-i]*fraci[M-N]%MOD
    if i%2:
        X -= t
    else:
        X += t
    X%=MOD
print(X*frac[M]*fraci[M-N]%MOD)

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

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

def modinv(a):
    b = mod
    u = 1
    v = 0
    while b:
        t = a // b
        a -= t * b
        a, b = b, a
        u -= t * v
        u, v = v, u
    u %= mod
    if u < 0:
        u += m
    return u

ans = fac[m] * fac[n] % mod
tmp = 0
for k in range(n + 1):
    tmp2 = fac[m - k]
    tmp2 = tmp2 * modinv(fac[n - k]) % mod
    tmp2 = tmp2 * modinv(fac[k]) % mod
    if k % 2 == 0:
        tmp = (tmp + tmp2) % mod
    else:
        tmp = (tmp - tmp2) % mod
ans = ans * tmp % mod
inv = modinv(fac[m - n])
ans = ans * inv % mod
ans = ans * inv % mod
print(ans)

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

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

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

f=[0,1]
for i in range(2,n+2):
  f.append(f[-1]*i%mod)

x=[1]
for i in range(m-n+1,m+1):
  x.append(x[-1]*i%mod)
x.reverse()
ans=0
for i in range(len(x)):
  if i%2==1:
    ans-=cmb(n,i,mod)*x[i]
  else:
    ans+=cmb(n,i,mod)*x[i]
  ans%=mod
print((ans*f[n]*cmb(m,n,mod))%mod)

Solution (Python):


LI = lambda: list(map(int, input().split()))


N, M = LI()

MOD = 10 ** 9 + 7
MAX = (10 ** 5) * 5 + 5
fac, finv, inv = [None] * MAX, [None] * MAX, [None] * MAX


def comb_init():
    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] = -inv[MOD%i] * int(MOD / i) % MOD
        finv[i] = finv[i - 1] * inv[i] % MOD


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


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


def main():
    comb_init()
    ans = 0
    v = perm(M, N)
    for i in range(N + 1):
        p = 1 if i % 2 == 0 else -1
        x = comb(N, i) * perm(M - i, N - i) % MOD
        ans = (ans + p * v * x) % MOD
    print(ans)
    

if __name__ == "__main__":
    main()

Solution (Python):
mod = 10**9+7
fac = [1]
fac_in = [1]
for i in range(1,5*10**5+1):
    fac.append((fac[-1]*i)%mod)
    fac_in.append((fac_in[-1]*pow(i,mod-2,mod))%mod)

def comb(a,b):
    return (fac[a]*fac_in[b]*fac_in[a-b])%mod

def perm(a,b):
    return (fac[a]*fac_in[a-b])%mod

n,m = map(int,input().split())
ans = perm(m,n)
ans1 = 0
for i in range(n+1):
    ans1 += (comb(n,i)*perm(m-i,n-i))*(-1)**i
    ans1 %= mod
ans *= ans1
print(ans%mod)

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

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

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

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

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

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

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

N,M = map(int,input().split())
ans = 0
for S in range(N+1):
    hoge = P(M,S)*(P(M-S,N-S)**2)*comb(N,S)    
    if S%2 == 1:
        hoge *= -1
    ans += hoge

print(ans%MOD)

Solution (Python):
import sys; input = sys.stdin.buffer.readline
sys.setrecursionlimit(10**7)
from collections import defaultdict
con = 10 ** 9 + 7; INF = float("inf")

def getlist():
	return list(map(int, input().split()))

class Combination(object):
	def __init__(self, N, con):
		self.fac = [0] * (N + 1)
		self.inv = [0] * (N + 1)
		self.finv = [0] * (N + 1)
		self.fac[0], self.fac[1] = 1, 1
		self.inv[1] = 1
		self.finv[0], self.finv[1] = 1, 1

		
		for i in range(2, N + 1):
			self.fac[i] = self.fac[i - 1] * i % con
			self.inv[i] = self.inv[con % i] * (con - (con // i)) % con
			self.finv[i] = self.finv[i - 1] * self.inv[i] % con

	def com(self, N, k):
		return (self.fac[N] * self.finv[k] * self.finv[N - k]) % con



def main():
	N, M = getlist()
	ComM = Combination(M, con)
	ComN = Combination(N, con)
	anspre = (ComM.fac[M] * ComM.finv[M - N]) % con

	anspre2 = 0
	for i in range(N + 1):
		val = (ComN.com(N, i) * ComM.fac[M - i] * ComM.finv[M - N]) % con
		if i % 2 == 0:
			anspre2 += val
		else:
			anspre2 -= val

	ans = anspre * anspre2 % con
	print(ans)

if __name__ == '__main__':
	main()

Solution (Python):
"""

手の付け方がわからない
ある条件を決める中での余事象？
A != Bでおいていく
毎回横で被る場合を引く

←にk個置かれてて
* MC2　したとする
上側に置いたのがかぶってる確率は (k / M)
下側もそう
両側は (k/M)**2

dp[k+1] = dp[k] * MC2 * (1 - 2*(k/M) + (k/M)**2)

Aを置いた
Bを置く　→　ほう助原理

"""

def inverse(a,mod): 
    return pow(a,mod-2,mod)

def modfac(n, MOD):
 
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(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 factorials, invs


def modnCr(n,r): 

    return fac[n] * inv[n-r] * inv[r] % mod


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

fac,inv = modfac(5*10**5+10,mod)

las = fac[N] * modnCr(M,N)

now = 0
for i in range(N+1):

    np = modnCr(N,i) * fac[N-i] * modnCr(M-i,N-i)
    if i % 2 == 0:
        now += np
    else:
        now -= np
    
    now %= mod


print (las * now % mod)

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


def inv(a, mod):
    r = [1, 0, a]
    w = [0, 1, mod]
    while w[2] != 1:
        q = r[2]//w[2]
        r_new = [r[0]-q*w[0], r[1]-q*w[1], r[2]-q*w[2]]
        r = w
        w = r_new
    x, y = w[0], w[1]    
    return (mod+x % mod) % mod


max_num = 5*10**5+2  
fact = [0 for _ in range(max_num)]
ifact = [0 for _ in range(max_num)]

fact[0] = fact[1] = 1
ifact[0] = ifact[1] = 1

for i in range(2, max_num):
    fact[i] = fact[i-1] * i % mod

ifact[max_num-1] = inv(fact[max_num-1], mod)

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


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


tmp = 0
for i in range(1, n+1):
    tmp2 = comb(n, i) * fact[m-i] * ifact[m-n]
    tmp2 %= mod
    tmp += (-1)**(i-1)*tmp2

tmp = fact[m]*ifact[m-n] - tmp
tmp %= mod

ans = tmp*fact[m]*ifact[m-n]
ans %= mod

print(ans)

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

mod = 10 ** 9 + 7

fac = [1]
for i in range(1, N + M + 10) :
  fac.append(fac[-1] * i % mod)
inv = [pow(fac[-1], mod - 2, mod)]
for i in range(N + M + 9, 0, -1) :
  inv.append(inv[-1] * i % mod)
inv = inv[::-1]

def comb(n, r) :
  if n >= r and r >= 0 :
    return fac[n] * inv[n - r] * inv[r] % mod
  return 1
  
ret = 0
for i in range(N + 1) :
  if i % 2 == 1 :
    sign = -1
  else :
    sign = 1
  x = sign * comb(M, i) * comb(N, i) * fac[i] % mod
  x *= (fac[M - i] ** 2) * (inv[M - N] ** 2) % mod
  ret += x % mod
  ret %= mod
print(ret)

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


mod = 1000000007
F = [1]
f = 1
for i in range(1,M+1):
    f *= i
    f %= mod
    F.append(f)

ans = 0

for k in range(N+1):
    a =  F[M-k]
    b = F[M-N]*F[N-k]*F[k]
   
    a *= pow(b,mod-2,mod)
    if k%2==0:
        ans += a
    else:
        ans += mod-a
    ans %= mod

ans *= F[N]
ans %= mod

for i in range(N):
    ans *= M-i
    ans %= mod
print(ans)

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

def make_fact(n, k):
    tmp = n
    perm = [i for i in range(k)]
    L = [0 for _ in range(k)]
    for i in range(k):
        L[i] = tmp % (i + 1)
        tmp //= i + 1
    LL = [0 for _ in range(k)]
    for i in range(k):
        LL[i] = perm[L[-i-1]]
        for j in range(L[-i-1]+1, k):
            perm[j-1] = perm[j]
    return LL

def doubling(n, m, modulo=mod):
    y = 1
    base = n
    tmp = m
    while tmp != 0:
        if tmp % 2 == 1:
            y *= base
            if modulo > 0:
                y %= modulo
        base *= base
        if modulo > 0:
            base %= modulo
        tmp //= 2
    return y

def inved(a, modulo=mod):
    x, y, u, v, k, l = 1, 0, 0, 1, a, modulo
    while l != 0:
        x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
        k, l = l, k % l
    return x % modulo

fact = [1 for _ in range(maxf+1)]
invf = [1 for _ in range(maxf+1)]

for i in range(maxf):
    fact[i+1] = (fact[i] * (i+1)) % mod
invf[-1] = inved(fact[-1])
for i in range(maxf, 0, -1):
    invf[i-1] = (invf[i] * i) % mod

S = 0
sign = 1
for i in range(N+1):
    S += (sign*fact[M-i]*invf[i]*invf[N-i])%mod
    S %= mod
    sign *= -1
    sign %= mod
print(S*fact[N]*fact[M]*invf[M-N]*invf[M-N]%mod)

Solution (Python):
from collections import defaultdict, deque, Counter
from heapq import heappush, heappop, heapify
import bisect
import random
from itertools import permutations, accumulate, combinations, product
import sys
import string
from bisect import bisect_left, bisect_right
from math import factorial, ceil, floor, gcd, sqrt
from operator import mul
from functools import reduce
from operator import mul
from pprint import pprint



sys.setrecursionlimit(2147483647)
INF = 10 ** 20
def LI(): return list(map(int, sys.stdin.buffer.readline().split()))
def I(): return int(sys.stdin.buffer.readline())
def LS(): return sys.stdin.buffer.readline().rstrip().decode('utf-8').split()
def S(): return sys.stdin.buffer.readline().rstrip().decode('utf-8')
def IR(n): return [I() for i in range(n)]
def LIR(n): return [LI() for i in range(n)]
def SR(n): return [S() for i in range(n)]
def LSR(n): return [LS() for i in range(n)]
def SRL(n): return [list(S()) for i in range(n)]
def MSRL(n): return [[int(j) for j in list(S())] for i in range(n)]
mod = 1000000007


n, m = LI()

fac = [1] * (m + 1)
inv = [1] * (m + 1)
for j in range(1, m + 1):
    fac[j] = fac[j-1] * j % mod


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


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

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

print(ret)

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

MOD=10**9+7
MAXN=500000+10
g1=[1, 1]
g2=[1, 1]
inverse=[0, 1]

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

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

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

t = 0
sign = 1
for i in range(1, n+1):
    t += sign*comb(n, i)*p(m-i, n-i)
    t %= MOD
    sign *= -1
print(p(m, n)*(p(m, n)-t)%MOD)

Solution (Python):
import sys
sys.setrecursionlimit(10**7) 
import math
from copy import copy, deepcopy
from copy import deepcopy as dcp
from operator import itemgetter
from bisect import bisect_left, bisect, bisect_right

from collections import deque


from collections import Counter

from itertools import accumulate,combinations,permutations

from heapq import heapify,heappop,heappush


from functools import lru_cache

from decimal import Decimal

def input(): 
    x=sys.stdin.readline()
    return x[:-1] if x[-1]=="\n" else x
def printl(li): _=print(*li, sep="\n") if li else None
def argsort(s, return_sorted=False): 
    inds=sorted(range(len(s)), key=lambda k: s[k])
    if return_sorted: return inds, [s[i] for i in inds]
    return inds
def alp2num(c,cap=False): return ord(c)-97 if not cap else ord(c)-65
def num2alp(i,cap=False): return chr(i+97) if not cap else chr(i+65)
def matmat(A,B):
    K,N,M=len(B),len(A),len(B[0])
    return [[sum([(A[i][k]*B[k][j]) for k in range(K)]) for j in range(M)] for i in range(N)]
def matvec(M,v):
    N,size=len(v),len(M)
    return [sum([M[i][j]*v[j] for j in range(N)]) for i in range(size)]
def T(M):
    n,m=len(M),len(M[0])
    return [[M[j][i] for j in range(n)] for i in range(m)]
def main():
    mod = 1000000007
    

    
    N, M = map(int, input().split())
    
    
    
    def extgcd1(a0,b0):
        u,v,a,b=1,0,a0,b0
        while b: t=a//b; a-=t*b; a,b=b,a; u,v=v,u-t*v
        if a!=1: return -1
        return u%b0
    
    maxn=M
    fact=[1]*(maxn+1)
    ifact=[1]*(maxn+1)
    x=1
    for i in range(2,maxn+1):
        x=(x*i)%mod
        fact[i]=x; ifact[i]=extgcd1(x,mod)
        
    def comb(n,r): return (fact[n]*ifact[r]%mod)*ifact[n-r]%mod
    def P(n,r): return (fact[n]*ifact[n-r]%mod)
    ans=0
    sin=-1
    for i in range(0,N+1):
        sin*=-1
        ans+=sin*comb(M,i)*P(N,i)%mod*(P(M-i,N-i)**2%mod)
        ans%=mod
        
        
    
    print(ans)




if __name__ == "__main__":
    main()

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

fact = [1]*(n+1)
rfact = [1]*(n+1)
r = 1
for i in range(1, n+1):
  fact[i] = r = r * i % mod
rfact[n] = r = pow(fact[n], mod-2, mod)
for i in range(n, 0, -1):
  rfact[i-1] = r = r * i % mod


def perm(n, k):
  return fact[n] * rfact[n-k] % mod


def comb(n, k):
  return fact[n] * rfact[k] * rfact[n-k] % mod

N,M = map(int, input().split())
a = perm(M,N)
b = 0
for i in range(N+1):
    if i % 2:
        b -= comb(N,i) * perm(M-i,N-i) % mod
    else:
        b += comb(N,i) * perm(M-i,N-i) % mod
b %= mod
ans = a * b % mod
print(ans)

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

fact = [0] * 1100000
invfact = [0] * 1100000

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

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

for i in range(1100000 - 1, 0, -1):
    invfact[i-1] = invfact[i] * i % MOD
    

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

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



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

ans = 0

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

print(ans)

Solution (Python):
import sys
import numpy as np

read = sys.stdin.buffer.read
readline = sys.stdin.buffer.readline
readlines = sys.stdin.buffer.readlines

MOD = 10**9 + 7

def fact_table(N, MOD):
    inv = np.empty(N, np.int64)
    inv[0] = 0
    inv[1] = 1
    for n in range(2, N):
        q, r = divmod(MOD, n)
        inv[n] = inv[r] * (-q) % MOD
    fact = np.empty(N, np.int64)
    fact[0] = 1
    for n in range(1, N):
        fact[n] = n * fact[n - 1] % MOD
    fact_inv = np.empty(N, np.int64)
    fact_inv[0] = 1
    for n in range(1, N):
        fact_inv[n] = fact_inv[n - 1] * inv[n] % MOD
    return fact, fact_inv, inv

def main(N, M, MOD):
    fact, fact_inv, inv = fact_table(N + 1, MOD)
    x = fact[N] * fact_inv % MOD * fact_inv[::-1] % MOD
    x[1::2] *= (-1)
    x = x[::-1]
    p = 1  
    ret = 0
    for k in range(N + 1):
        ret += x[k] * p % MOD
        if k != N:
            p = p * (M - N + 1 + k) % MOD
    ret %= MOD
    return ret * p % MOD

if sys.argv[-1] == 'ONLINE_JUDGE':
    import numba
    from numba.pycc import CC
    i8 = numba.int64
    cc = CC('my_module')

    def cc_export(f, signature):
        cc.export(f.__name__, signature)(f)
        return numba.njit(f)

    fact_table = cc_export(fact_table, (i8, i8))
    main = cc_export(main, (i8, i8, i8))
    cc.compile()

from my_module import main

N, M = map(int, read().split())
print(main(N, M, MOD))

Solution (Python):
import sys

int1 = lambda x: int(x) - 1
p2D = lambda x: print(*x, sep="\n")
def II(): return int(sys.stdin.readline())
def MI(): return map(int, sys.stdin.readline().split())
def LI(): return list(map(int, sys.stdin.readline().split()))
def LLI(rows_number): return [LI() for _ in range(rows_number)]
def SI(): return sys.stdin.readline()[:-1]

def nCr(com_n, com_r):
    if com_n < com_r: return 0
    return fac[com_n] * ifac[com_r] % md * ifac[com_n - com_r] % md

def nPr(com_n, com_r):
    if com_n < com_r: return 0
    return fac[com_n] * ifac[com_n - com_r] % md

md = 10 ** 9 + 7
n_max = 500005
fac = [1]
for i in range(1, n_max + 1): fac.append(fac[-1] * i % md)
ifac = [1] * (n_max + 1)
ifac[n_max] = pow(fac[n_max], md - 2, md)
for i in range(n_max - 1, 1, -1): ifac[i] = ifac[i + 1] * (i + 1) % md

def main():
    n,m=MI()
    ans=0
    for i in range(n+1):
        if i&1:ans-=nCr(n,i)*nPr(m,i)*nPr(m-i,n-i)**2
        else:ans+=nCr(n,i)*nPr(m,i)*nPr(m-i,n-i)**2
        ans%=md
    print(ans)

main()

Solution (Python):
import sys
input=lambda: sys.stdin.readline().rstrip()
n,m=map(int,input().split())
mod=10**9+7
n_max=(10**6+1)
F,FI=[0]*(n_max+1),[0]*(n_max+1)
F[0],FI[0]=1,1
for i in range(n_max):
  F[i+1]=(F[i]*(i+1))%mod
FI[n_max-1]=pow(F[n_max-1],mod-2,mod)
for i in reversed(range(n_max-1)):
  FI[i]=(FI[i+1]*(i+1))%mod
def comb(x,y):
  return (F[x]*FI[x-y]*FI[y])%mod

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

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

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

def choose(n,r): 
    if 0 <= r <= n:
        return (fac[n]*finv[r]%MOD)*finv[n-r]%MOD
    else:
        return 0

def f(x,y):
    return fac[x]*finv[x-y]%MOD
    
    


import sys
read = sys.stdin.read
readline = sys.stdin.readline

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

ans = 0
r = 1
for u in range(n+1):
    ans += r*f(m-u,n-u)*choose(n,u)%MOD
    ans %= MOD
    r *= -1



print(ans*f(m,n)%MOD)

Solution (Python):
mod = 1000000007
eps = 10**-9


def main():
    import sys
    input = sys.stdin.readline

    
    
    nmax = 10 ** 6 + 10  
    fac = [0] * nmax
    finv = [0] * nmax
    inv = [0] * nmax
    fac[0] = 1
    fac[1] = 1
    finv[0] = 1
    finv[1] = 1
    inv[1] = 1
    for i in range(2, nmax):
        fac[i] = fac[i - 1] * i % mod
        inv[i] = mod - inv[mod % i] * (mod // i) % mod
        finv[i] = finv[i - 1] * inv[i] % mod

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

    N, M = map(int, input().split())
    ans = 1
    for i in range(N):
        ans = (ans * (M-i))%mod
    ans = (ans * ans)%mod

    mpn = fac[M] * finv[M-N]
    for k in range(1, N+1):
        tmp = ((comb(N, k) * mpn)%mod * (fac[M-k] * finv[M-N])%mod)%mod
        ans = (ans + tmp * (-1)**k) % mod
    print(ans)


if __name__ == '__main__':
    main()

Solution (Python):
def prepare(n, MOD):
    f = 1
    factorials = [1]
    for m in range(1, n + 1):
        f *= m
        f %= MOD
        factorials.append(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 factorials, invs


n, m = map(int, input().split())
MOD = 10 ** 9 + 7
facts, finvs = prepare(m, MOD)
c = 1
ans = 0
base = facts[n] * facts[m] * finvs[m - n] * finvs[m - n] % MOD
for k in range(n + 1):
    tmp = base * facts[m - k] * finvs[k] * finvs[n - k] % MOD
    ans = (ans + c * tmp) % MOD
    c *= -1
print(ans)