Problem:
 Score : 500 points    Problem Statement Count the pairs of length- N sequences consisting of integers between 1 and M (inclusive), A_1, A_2, \cdots, A_{N} and B_1, B_2, \cdots, B_{N} , that satisfy all of the following conditions:   A_i \neq B_i , for every i such that 1\leq i\leq N .  A_i \neq A_j and B_i \neq B_j , for every (i, j) such that 1\leq i < j\leq N .   Since the count can be enormous, print it modulo (10^9+7) .      Constraints  1\leq N \leq M \leq 5\times10^5  All values in input are integers.         Input Input is given from Standard Input in the following format:  N  M       Output Print the count modulo (10^9+7) .        Sample Input 1 2 2      Sample Output 1 2  A_1=1,A_2=2,B_1=2,B_2=1 and A_1=2,A_2=1,B_1=1,B_2=2 satisfy the conditions.       Sample Input 2 2 3      Sample Output 2 18       Sample Input 3 141421 356237      Sample Output 3 881613484

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

class mod_pack:

    def __init__(self, MOD, MAX):
        self.MOD = MOD
        self.MAX = MAX
        self.fac, self.ifac = [1]*(MAX+1), [1]*(MAX+1)
        for i in range(1, MAX+1):
            self.fac[i] = self.fac[i-1]*i % MOD
            self.ifac[i] = self.ifac[i-1] * self.mpow(i, MOD-2) % MOD
    
    def mpow(self, x, n):
        ans = 1
        while n!=0:
            if n&1: ans = ans*x % self.MOD
            x = x*x % self.MOD
            n = n >> 1
        return ans

    def mcomb(self, m, n):
        if (m==0 and n==0): return 1
        return (self.ifac[m-n] * self.ifac[n] % self.MOD) * self.fac[m] % self.MOD

    def mperm(self, m, n):
        if (m < n): return 0
        return self.ifac[m-n] * self.fac[m] % MOD
    

mp = mod_pack(10**9+7, M)

tot = mp.mperm(M, N)
ans = 0
for k in range(1, N+1):

    ans = (ans + (-1)**(k-1) * mp.mcomb(N, k) * mp.mperm(M-k, N-k)) % MOD


print((tot * (tot - ans))%MOD)

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

MOD = 10**9 + 7





def solve():
    N,M = map(int,input().split())
    
    fac = [1,1] + [0] * M
    inv = [1,1] + [0] * M
    for i in range(2,M+1):
        fac[i] = (fac[i-1]*i)%MOD
    inv[M] = pow(fac[M],MOD-2,MOD)
    for i in range(M-1,1,-1):
        inv[i] = inv[i+1]*(i+1)%MOD
    def p(a,b):
        return (fac[b]*inv[b-a])%MOD
    def c(a,b):
        return (fac[b]*inv[a]*inv[b-a]) % MOD
    res = p(N,M)
    for i in range(1,N+1):
        res += (-1)**(i&1)*c(i,N)*p(N-i,M-i)
        res %= MOD
    print(res* p(N,M)%MOD)




if __name__ == '__main__':
    solve()

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

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

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

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

ans = 0
for k in range(N+1):
    ans += cmb(N, k, mod)*(-1)**k*cmb(M, k, mod)*g1[k]*(cmb(M-k, N-k, mod)*g1[N-k])**2
    ans %= mod
print(ans)

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

MOD = 10**9+7
NUM = 10**6
g1 = [1, 1]
g2 = [1, 1]
inv = [0, 1]

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

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

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

ans = perm(m,n)**2
sub=0

for k in range(1,n+1):
    cnt = comb(n,k)*perm(m-k,n-k) % MOD
    if k%2==1:
        sub+=cnt
    else:
        sub-=cnt

sub*=perm(m,n)%MOD
ans-=sub
print(ans%MOD)

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

for i in range(2, m+1):
    inv[i] = (-inv[mod % i]*(mod//i)) % mod
    invf[i] = invf[i-1]*inv[i] % mod

ans = 0


def chs(x, y):
    return fct[x]*invf[y]*invf[x-y]
def chs2(x,y):
    return fct[x]*invf[x-y]

m_inv = pow(m, mod-2, mod)
mn = pow(m, n, mod)
for i in range(n+1):
    ans += ((-1)**(i % 2))* chs(n, i)*chs2(m-i,n-i)
    ans %= mod
    mn *= (m_inv)
    mn %= mod
    
print(ans*chs2(m,n)%mod)

Solution (Python):

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

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

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

alphabets = "abcdefghijklmnopqrstuvwxyz"
ALPHABETS = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"

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


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

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

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

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

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

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


class ModInt:
    def __init__(self, x, MOD=1000000007): self.x, self.MOD = x % MOD, MOD
    def __str__(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 * other.inverse()) if isinstance(other, ModInt) else ModInt(self.x * pow(other, self.MOD - 2, self.MOD))
    def __pow__(self, other): return ModInt(pow(self.x, other.x, self.MOD)) if isinstance(other, ModInt) else ModInt(pow(self.x, other, self.MOD))
    def __rsub__(self, other): return ModInt(other.x - self.x) if isinstance(other, ModInt) else ModInt(other - self.x)
    def __rtruediv__(self, other): return ModInt(other.x * other.inverse()) if isinstance(other, ModInt) else ModInt(other * pow(self.x, self.MOD - 2, self.MOD))
    def __rpow__(self, other): return ModInt(pow(other.x, self.x, self.MOD)) if isinstance(other, ModInt) else ModInt(pow(other, self.x, self.MOD))
    __repr__, __radd__, __rmul__ = __str__, __add__, __mul__
    def inverse(self): return pow(self.x, self.MOD - 2, self.MOD)

def main():
    N, M = LI()
    factinv = FactInv(M)
    a = ModInt(factinv.perm(M, N))
    tmp = ModInt(0)

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

    print(a * tmp)





if __name__ == '__main__':
    main()

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

def invmod(a):
  return pow(a,MOD-2,MOD)

fact_list=[1]
fact_inv_list=[1]
for i in range(1,M+1):
  fact_mod=(fact_list[-1]*i)%MOD
  fact_list.append(fact_mod)
  fact_inv_list.append(invmod(fact_mod))
def comb_mod_table(n,r):
  if 0<=r<=n:
    return fact_list[n]*fact_inv_list[r]*fact_inv_list[n-r]
  else:
    return 0
def perm_mod_table(n,r):
  if 0<=r<=n:
    return fact_list[n]*fact_inv_list[n-r]
  else:
    return 0
  
answer=0
for i in range(N+1):  
  term=comb_mod_table(N,i)*pow(-1,i)*perm_mod_table(M,i)*perm_mod_table(M-i,N-i)*perm_mod_table(M-i,N-i)
  answer=(answer+term)%MOD
  
print(answer)

Solution (Python):

def factorials_with_inv(k,mod):
    """
    0! ... k! とそれらのmod逆元を求めて返す
    """
    fac = [1]*(k+1)
    inv = [1]*(k+1)
    t = 1
    for i in range(1,k+1):
        t *= i
        t %= mod
        fac[i] = t

    t = pow(t,-1,mod)
    for i in reversed(range(1,k+1)):
        inv[i] = t
        t *= i
        t %= mod

    return fac, inv

MOD = 10**9+7

def solve(N,M):
    fac,inv = factorials_with_inv(max(N,M), MOD)

    

    res = 0
    sign = True
    for k in range(N+1):
        t = fac[M-k]*inv[k]*inv[N-k]
        if sign:
            res += t
        else:
            res -= t
        res %= MOD
        sign = not sign

    res *= fac[N]*fac[M]*inv[M-N]**2
    res %= MOD
    return res

from itertools import permutations

def naive(N,M):
    res = 0
    for A in permutations(range(M),N):
        for B in permutations(range(M),N):
            res += all(a != b for a,b in zip(A,B))
    return res%MOD

if __name__ == '__main__':
    N,M = map(int,input().split())

    
    print(solve(N,M))

Solution (Python):
class Combination:
    def __init__(self, n_max=10**6, mod=10**9+7):
        
        self._fac, self._finv, self._inv = [0]*n_max, [0]*n_max, [0]*n_max
        self._fac[0], self._fac[1] = 1, 1
        self._finv[0], self._finv[1] = 1, 1
        self._inv[1] = 1
        self._mod = mod
        for i in range(2, n_max):
            self._fac[i] = self._fac[i - 1] * i % self._mod
            self._inv[i] = self._mod - self._inv[self._mod%i] * (self._mod // i) % self._mod
            self._finv[i] = self._finv[i - 1] * self._inv[i] % self._mod
    def com(self, n, r):
        if n < r: return 0
        if n < 0 or r < 0: return 0
        return self._fac[n] * (self._finv[r] * self._finv[n - r] % self._mod) % self._mod

    def perm(self,n,r):
        if n < r: return 0
        if n < 0 or r < 0: return 0
        return self._fac[n] * (self._finv[n-r] % self._mod) % self._mod


MOD = 10**9+7
comb = Combination(10**6, MOD)

pls = [-1]*(10**6)
n,m = map(int, input().split())

ans = 0
for i in range(1,n+1):
    i_com = comb.com(n,i)
    pattern = comb.perm(m-i,n-i)
    if i%2 == 1:
        ans += i_com*pattern
    else:
        ans -= i_com*pattern
    ans%=MOD

ans = comb.perm(m,n)-ans
ans%=MOD

ans *= comb.perm(m,n)
print(ans%MOD)

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

fact = [0 for _ in range(MAX)]
factinv = [0 for _ in range(MAX)]

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

factinv[MAX - 1] = pow(fact[MAX - 1], MOD - 2, MOD)
for k in range(MAX - 1, 0, -1):
    factinv[k - 1] = factinv[k]*k
    factinv[k - 1] %= MOD

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

def nPk(n, k):
    return fact[n]*factinv[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)*nPk(m - k, n - k) % MOD
    if not k % 2: ans += tmp
    else: ans -= tmp
    ans %= MOD

print(ans)

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

fact = [0 for _ in range(MAX)]
factinv = [0 for _ in range(MAX)]

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

factinv[MAX - 1] = pow(fact[MAX - 1], MOD - 2, MOD)
for k in range(MAX - 1, 0, -1):
    factinv[k - 1] = factinv[k]*k
    factinv[k - 1] %= MOD

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

def nPk(n, k):
    return fact[n]*factinv[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)*nPk(m - k, n - k) % MOD
    if not k % 2: ans += tmp
    else: ans -= tmp
    ans %= MOD

print(ans)

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

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

def init(n):
    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)

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

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

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

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

num = 500010
fact = [1] * (num+1)
ifact = [1] * (num+1)

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

def nPr(n,r):
    if r > n:
        return 0
    return (fact[n] * ifact[n-r]) % MOD
def nCr(n,r):
    if r > n:
        return 0
    return (fact[n] * ifact[r] * ifact[n-r]) % MOD

ans = nPr(M,N) ** 2
sign = -1
for n in range(1,N+1):
    t = (nPr(M,n) * nCr(N,n) * (nPr(M-n,N-n) ** 2))
    ans += sign * t
    sign *= -1
    ans %= MOD
print(ans)

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


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


def comb(n, r):
    global MOD, modFacts, invs
    return (modFacts[n] * invs[n - r] * invs[r]) % MOD


def perm(n, r):
    global MOD, modFacts, invs
    return (modFacts[n] * invs[n - r]) % MOD


N, M = map(int, input().split())
modFacts, invs = prepare(M)
ans = 0
flag = 1
for i in range(N + 1):
    a = comb(N, i)
    b = perm(M, i)
    c = pow(perm(M - i, N - i), 2, MOD)
    cnt = (a * b * c) % MOD
    ans = (ans + flag * cnt) % MOD
    flag *= -1
print(ans)

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

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

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

ans = 0
for k in range(N+1):
  if k%2 == 0:
    ans = (ans+factinv[N-k]*factinv[k]%mod*fact[M-k]%mod)%mod
  else:
    ans = (ans-factinv[N-k]*factinv[k]%mod*fact[M-k]%mod)%mod
ans = ans*fact[N]%mod*fact[M]%mod*factinv[M-N]%mod*factinv[M-N]%mod
print(ans)

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, readline().split())

ans = 0

for i in range(N+1):
    sign = -1 if i&1 else 1
    ans = (ans + sign*comb(M, i)*comb(N, i)*frac[i]*comb(M-i, N-i)**2*frac[N-i]**2)%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

def p(n,r,mod):
    return fact[n]*factinv[n-r] % mod
mod = 10 ** 9 + 7
maxi = 5*10 ** 5  
fact = [1, 1]  
factinv = [1, 1]  
inv = [0, 1]  
 
for i in range(2, maxi + 1):
    fact.append((fact[-1] * i) % mod)
    inv.append((-inv[mod % i] * (mod // i)) % mod)
    factinv.append((factinv[-1] * inv[-1]) % mod)

N,M=map(int,input().split())    
    
ans=0
for i in range(N+1):
    ans += cmb(N,i,mod) * (-1)**i * p(M,i,mod) * (p(M-i,N-i,mod))**2
    ans %= mod
print(ans)

Solution (Python):
N,M=map(int,input().split())
list=[0 for i in range(5*10**5+1)]
mod=10**9+7
list[0]=1
for i in range(len(list)-1):
  list[i+1]=list[i]*(i+1)%mod
def c(m,n):
  a=list[m]*pow(list[m-n],mod-2,mod)*pow(list[n],mod-2,mod)
  return a
def p(m,n):
  a=list[m]*pow(list[m-n],mod-2,mod)
  return a
ans=pow(p(M,N),2)
g=pow(list[M-N],mod-2,mod)
g2=pow(g,2,mod)
nn=list[N]
mm=list[M]
nmg=nn*mm*g2%mod
for i in range(1,N+1):
  ans+=nmg*pow(list[N-i],mod-2,mod)*pow(list[i],mod-2,mod)*list[M-i]*(-1)**i
  ans%=mod
print(ans)

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

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


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

ans = 0

for k in range(N+1):
    a = (-1)**k
    a *= R[N]*R[M]*R[M-k]
    a %= mod
    a *= pow(R[k],mod-2,mod)
    a *= pow(R[N-k],mod-2,mod)
    a %= mod
    a *= (pow(R[M-N],mod-2,mod))**2
    a %= mod
    ans += a

print(ans%mod)

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

    def 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



N, M = map(int, input().split())
MOD=10 ** 9 + 7
C = Combination(M, MOD)
sgn = 1
res = 0
for k in range(0, N+1):
    res += sgn*C.C(N,k)*C.P(M-k,N-k)
    res %= MOD
    sgn *= -1
print(res*C.P(M,N)%MOD)

Solution (Python):
from functools import reduce

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

mod = 10 ** 9 + 7
ans = 0
c = 1
p1 = 1
p2 = reduce(lambda ac, v: ac * v % mod, range(m, m - n, -1), 1)
for i in range(n + 1):
    if i > 0:
        c = (c * (n - i + 1) * pow(i, mod - 2, mod)) % mod
        p1 = (p1 * (m - i + 1)) % mod
        p2 = (p2 * pow(m - i + 1, mod - 2, mod)) % mod
    x = c * p1 * (p2 ** 2)
    if i % 2 == 0:
        ans += x
    else:
        ans -= x
    ans %= mod
print(ans)

Solution (Python):
MOD=10**9+7
def facinv(N):
    fac,finv,inv=[0]*(N+1),[0]*(N+1),[0]*(N+1)
    fac[0]=1;fac[1]=1;finv[0]=1;finv[1]=1;inv[1]=1
    for i in range(2,N+1):
        fac[i]=fac[i-1]*i%MOD
        inv[i]=MOD-inv[MOD%i]*(MOD//i)%MOD
        finv[i]=finv[i-1]*inv[i]%MOD
    return fac,finv,inv

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

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

def p(r):
    return 1 if r%2 else -1

N,M=map(int,input().split())
fac,finv,inv=facinv(M)

ans=PER(M,N)**2%MOD
res=0
for r in range(1,N+1):
    f=COM(N,r)*PER(M,r)%MOD*(PER(M-r,N-r)**2%MOD)%MOD
    res=(res+(p(r)*f)+MOD)%MOD
print((ans-res+MOD)%MOD)

Solution (Python):
MAX = 5*10**5+100
MOD = 10**9+7
fact = [0]*MAX 
inv = [0]*MAX 
finv = [0]*MAX 
fact[0] = 1
fact[1] = 1
finv[0] = 1
finv[1] = 1
inv[1] = 1
    
for i in range(2, MAX):
    fact[i] = fact[i-1]*i%MOD
    inv[i] = MOD-inv[MOD%i]*(MOD//i)%MOD
    finv[i] = finv[i-1]*inv[i]%MOD

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

N, M = map(int, input().split())
ans = C(M, N)*fact[N]

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

Solution (Python):
MOD=10**9+7
def facinv(N):
    fac,finv,inv=[0]*(N+1),[0]*(N+1),[0]*(N+1)
    fac[0]=1;fac[1]=1;finv[0]=1;finv[1]=1;inv[1]=1
    for i in range(2,N+1):
        fac[i]=fac[i-1]*i%MOD
        inv[i]=MOD-inv[MOD%i]*(MOD//i)%MOD
        finv[i]=finv[i-1]*inv[i]%MOD
    return fac,finv,inv

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

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

N,M=map(int,input().split())
fac,finv,inv=facinv(M)

ans=PER(M,N)**2%MOD
res=0
for r in range(1,N+1):
    f=COM(N,r)*PER(M,r)%MOD*(PER(M-r,N-r)**2%MOD)%MOD
    res=(res+(pow(-1,r-1)*f)+MOD)%MOD
print((ans-res+MOD)%MOD)

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

def cmb(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):
    return g1[n] * g2[n-r] % MOD

maxN = 5*10**5
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 solve():
    N, M = map(int, input().split())
    ans = 0
    for k in range(N+1):
        x = 0
        tmp = p(M-k, N-k)
        x = cmb(N, k) * p(M, k) * tmp * tmp
        if k%2==1:
            x *= -1
        ans = (ans + x) % MOD
    print(ans)

solve()

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

def calc(n):
    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

factorials, invs = calc(M)
mPn = (factorials[M] * invs[M - N]) % mod
ans = (mPn * mPn) % mod
for k in range(1, N + 1):
    v = pow(-1, k - 1) * factorials[N] * factorials[M - k] * invs[N - k] * invs[k] * invs[M - N] % mod
    ans -= v * mPn
print(ans % mod)

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(max(N, M))

ans = 0
for i in range(N + 1):
    Ti = (modFacts[N] * invs[i] * invs[N - i]) % MOD
    Ti *= (modFacts[M] * invs[M - i]) % MOD
    Ti %= MOD
    Ti *= pow(modFacts[M - i] * invs[(M - i) - (N - i)], 2, MOD)
    Ti %= MOD
    ans += pow(-1, i) * Ti
    ans %= MOD

print(ans)

Solution (Python):



def p(n, k, mod):
  x = 1
  for ik in range(k):
    x *= (n - ik)
    x %= mod
  return x

MAX = 10 ** 9 + 7
def main():
  n, m = map(int, input().strip().split())
  s = 0
  nck = 1
  mpk = 1
  mkpnk = p(m, n, MAX)
  for k in range(0, n+1):
    if 0 < k <= n:
      nck *= (n - (k-1)) * pow(k, MAX-2, MAX) 
      nck %= MAX
    if 0 < k <= m:
      mpk *= (m - (k-1))
      mpk %= MAX
    if 0 < k <= n:
      mkpnk *= pow(m - (k - 1), MAX-2, MAX) 
      mkpnk %= MAX
    s += nck * ( (-1) ** k) * mpk * (mkpnk) ** 2
    s %= MAX
  print(s)
  

if __name__=='__main__':
  main()

Solution (Python):
import sys
import numpy as np

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

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 = map(int, readline().split())
    mod = 10**9 + 7
    fact, fact_inv, inv = fact_table(n+1, mod)
    comb = fact[n] * fact_inv % mod * fact_inv[::-1] % mod
    comb[1::2] *= (-1)
    comb = comb[::-1]

    p = 1
    ans = 0
    for k in range(n+1):
        ans += comb[k] * p % mod
        if k != n:
            p = p * (m-n+1+k)%mod
    ans %= mod
    print( p * ans % mod )
    return

main()

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

_fact = [1,1] + [0] * MAX
_inv = [0,1] + [0] * MAX
_fact_inv = [1,1] + [0] * MAX
for k in range(2, MAX + 1):
    _fact[k] = _fact[k-1] * k % MOD
    _inv[k] = (-_inv[MOD % k] * (MOD // k)) % MOD
    _fact_inv[k] = _fact_inv[k-1] * _inv[k] % MOD

def P(m, n):
    return _fact[m] * _fact_inv[m-n] % MOD

def C(m, n):
    return _fact[m] * _fact_inv[m-n] * _fact_inv[n] % MOD

def solve(n, m):
    res = 0
    for k in range(n+1):
        sign = 1 - 2 * (k%2)
        cur = sign * C(n, k) * P(m, k) * (P(m-k, n-k) ** 2)
        res = (res + cur) % MOD
    return res

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

Solution (Python):
mod = 10 ** 9 + 7
def comb(n, r):
    if n < r:return 0
    if n < 0 or r < 0:return 0
    return fa[n] * fi[r] % mod * fi[n - r] % mod
def per(n, r):
    if n < r:return 0
    if n < 0 or r < 0:return 0
    return fa[n] * fi[n - r] % mod
n, m = map(int, input().split())
fa = [1] * (m + 1)
fi = [1] * (m + 1)
for i in range(1, m + 1):
    fa[i] = fa[i - 1] * i % mod
    fi[i] = pow(fa[i], mod - 2, mod)
ans = 0
for i in range(n + 1):
    if i % 2 == 0:
        ans += comb(n, i) * per(m, i) % mod * (per(m - i, n - i) ** 2) % mod
        ans %= mod
    else:
        ans += comb(n, i) * -1 * per(m, i) % mod * (per(m - i, n - i) ** 2) % mod
        ans % mod
print(ans % mod)

Solution (Python):
N, M=list(map(int,input().split()))
p=10**9+7
fact = [1, 1]  
factinv = [1, 1]  
inv = [0, 1]  
for i in range(2, M+1):
    fact.append((fact[-1] * i) % p) 
    inv.append((-inv[p % i] * (p // i)) % p)
    factinv.append((factinv[-1] * inv[-1]) % p)
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
C=0
for s in range(N+1):
    C=(C+(-1)**s*cmb(N,s,p)*pmt(M,N,p)*pmt(M-s,N-s,p)) %p
print(C)

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

mod = 10**9+7
N = 10**6  
fact  = [1, 1]  
fact_inv = [1, 1]  
inv = [0, 1]  


def comb(n, r, mod):
    if (r<0 or r>n):
        return 0
    r = min(r, n-r)
    return fact[n] * fact_inv[r] * fact_inv[n-r] % mod

def perm(n, r, mod):
    if (r<0 or r>n):
        return 0
    return fact[n] *  fact_inv[n-r] % mod


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

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

print(ans)

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



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

mod_fact_inv = [pow(x, MOD-2, MOD) for x in mod_factorial]

def modperm(n, k):
    return mod_factorial[n] * mod_fact_inv[n-k] % MOD

def modcomb(n, k):
    return modperm(n, k) * mod_fact_inv[k] % MOD

sigma = 0
for i in range(N+1):
    x = modcomb(N, i)
    x %= MOD
    x *= modperm(M, i)
    x %= MOD
    x *= modperm(M-i, N-i)
    x %= MOD
    x *= modperm(M-i, N-i)
    x %= MOD
    x *= -1 if i&1 else 1
    sigma += x

sigma %= MOD
print(sigma)

Solution (Python):
mod=10**9+7
n=10**6
inv_t=[0]+[1]
for i in range(2,n):
  inv_t+=[inv_t[mod%i]*(mod-int(mod/i))%mod]

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

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

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

N,M=map(int, input().split())
def f(x):
  e=cmb(N,x)*prm(M,x)*prm(M-x,N-x)*prm(M-x,N-x)
  e%=mod
  return e

ans=0
for i in range(N+1):
  g=f(i)*((-1)**i)
  ans+=g
print(ans%mod)

Solution (Python):
mod=10**9+7
n=10**6
inv_t=[0]+[1]
for i in range(2,n):
  inv_t+=[inv_t[mod%i]*(mod-int(mod/i))%mod]

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

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

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

N,M=map(int, input().split())
def f(x):
  e=cmb(N,x)*prm(M,x)*prm(M-x,N-x)*prm(M-x,N-x)
  e%=mod
  return e

ans=0
for i in range(N+1):
  g=f(i)*((-1)**i)
  ans+=g
print(ans%mod)

Solution (Python):
MOD = 1000000007
class mint:
    def __init__(self,x=0):
        self.x = x % MOD
    def __str__(self):
        return str(self.x)

    __repr__ = __str__

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

    def __sub__(self, other):
        return(
            mint(self.x - other.x) if isinstance(other, mint) else mint(self.x - other)
        )
    
    def __mul__(self, other):
        return(
            mint(self.x * other.x) if isinstance(other, mint) else mint(self.x * other)
        )

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

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


    __radd__ = __add__

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

    __rmul__ = __mul__

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

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

MAX = 1000005
fac = [1,1]
finv = [1,1]
inv = [0,1]
def COMinit():
    for i in range(2,MAX):
        fac.append((fac[-1] * i) % MOD)
        inv.append((-inv[MOD%i] * (MOD//i)) % MOD)
        finv.append((finv[-1] * inv[-1]) % MOD)

def 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 

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

N,M = map(int, input().split())
COMinit()
res = mint(0)
now = mint()
for i in range(N+1):
    now = COM(N, i)
    now *= P(M-i, N-i)
    if(i%2 == 1):
        now = -1 * now
    res += now
res *= P(M, N)
print(res)

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


F = 10**6+5
mod = 10**9+7
fact = [1]*F
inv = [1]*F
for i in range(2,F):
    fact[i]=(fact[i-1]*i)%mod
inv[-1]=pow(fact[-1],mod-2,mod)
for i in range(F-2,1,-1):
    inv[i] = (inv[i+1]*(i+1))%mod
def comb(n,k):
    ans=(fact[n]*inv[n-k])%mod
    return (ans*inv[k])%mod

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

print(ans)