module module_modint
  use, intrinsic :: iso_fortran_env
  implicit none
  integer(int64) :: modu = 1000000007
  type modint
    integer(int64), private :: val = 0
  contains
    procedure, private :: asgn => asgn
    procedure, private :: add => add
    procedure, private :: sub => sub
    procedure, private :: mul => mul
    procedure, private :: div => div
    procedure, private :: pow => pow
    procedure, public :: inv => inv
    procedure, public :: get => get
    generic, public :: assignment(=) => asgn
    generic, public :: operator(+) => add
    generic, public :: operator(-) => sub
    generic, public :: operator(*) => mul
    generic, public :: operator(/) => div
    generic, public :: operator(**) => pow
  end type
  interface modint
    module procedure :: new_modint, new_modint_with_num
  end interface modint
  interface operator(+)
    module procedure :: add_r
  end interface
  interface operator(-)
    module procedure :: sub_r
  end interface
  interface operator(*)
    module procedure :: mul_r
  end interface
  interface operator(/)
    module procedure :: div_r
  end interface
  interface dot_product
    module procedure :: dot_prod
  end interface dot_product
  interface matmul
    module procedure :: matmul1, matmul2, matmul3
  end interface matmul
contains
  pure elemental integer(int64) function to_int64(num) result(res)
    class(*), intent(in) :: num
    select type (num)
    type is (integer(int64))
      res = num
    type is (integer(int32))
      res = int(num, int64)
    type is (integer(int16))
      res = int(num, int64)
    type is (integer(int8))
      res = int(num, int64)
    class default
      res = 0_8
      res = 1_8 / res
    end select
  end
  pure elemental type(modint) function new_modint()
  end
  pure elemental type(modint) function new_modint_with_num(num) result(res)
    class(*), intent(in) :: num
    res%val = to_int64(num)
  end
  pure elemental subroutine asgn(this, other)
    class(modint), intent(inout) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      this%val = other%val
    class default
      val = to_int64(other)
      if (0_8 <= val .and. val < modu) then
        this%val = val
      else
        this%val = modulo(val, modu)
      end if
    end select
  end
  pure elemental type(modint) function add(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      if (0_8 <= this%val + other%val .and. this%val + other%val < modu) then
        res%val = this%val + other%val
      else
        res%val = this%val + other%val - modu
      end if
    class default
      val = to_int64(other)
      if (0_8 <= this%val + val .and. this%val + val < modu) then
        res%val = this%val + val
      else
        res%val = modulo(this%val + val, modu)
      end if
    end select
  end
  pure elemental type(modint) function sub(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      if (0_8 <= this%val - other%val .and. this%val - other%val < modu) then
        res%val = this%val - other%val
      else
        res%val = this%val - other%val + modu
      end if
    class default
      val = to_int64(other)
      if (0_8 <= this%val - val .and. this%val - val < modu) then
        res%val = this%val - val
      else
        res%val = modulo(this%val - val, modu)
      end if
    end select
  end
  pure elemental type(modint) function mul(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: val
    select type (other)
    type is (modint)
      val = other%val
    class default
      val = to_int64(other)
    end select
    if (0_8 <= this%val * val .and. this%val * val < modu) then
      res%val = this%val * val
    else
      res%val = modulo(this%val * val, modu)
    end if
  end
  pure elemental type(modint) function div(this,other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    type(modint) :: tmp
    call asgn(tmp, other)
    call asgn(res, mul(this, inv(tmp)))
  end
  pure elemental integer(int64) function get(this) result(res)
    class(modint), intent(in) :: this
    res = this%val
  end
  pure elemental type(modint) function inv(this) result(res)
    class(modint), intent(in) :: this
    integer(int64) :: a, b, c, n
    integer(int64) :: x, y, z, m
    a = this%val
    b = modu
    c = 0_8
    n = 1_8
    do while (b /= 0_8)
      x = b
      y = mod(a,b)
      z = n-a/b*c
      m = c
      a = x
      b = y
      c = z
      n = m
    end do
    call asgn(res,n)
  end
  pure elemental type(modint) function pow(this, other) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    integer(int64) :: n
    type(modint) :: i
    res%val = 1
    i%val = this%val
    n = to_int64(other)
    do while (n > 0_8)
      if (btest(n, 0)) call asgn(res, mul(res, i))
      call asgn(i, mul(i, i))
      n = rshift(n, 1)
    end do
  end
  pure elemental type(modint) function add_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, add(this, other))
  end
  pure elemental type(modint) function sub_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, sub(this, other))
  end
  pure elemental type(modint) function mul_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, mul(this, other))
  end
  pure elemental type(modint) function div_r(other, this) result(res)
    class(modint), intent(in) :: this
    class(*), intent(in) :: other
    call asgn(res, div(this, other))
  end
  pure type(modint) function dot_prod(x, y) result(res)
    type(modint), intent(in) :: x(:), y(:)
    integer :: i
    if (size(x, 1) /= size(y, 1)) i = to_int64('')
    do i = 1, size(x,1)
      call asgn(res, add(res, mul(x(i), y(i))))
    end do
  end
  pure function matmul1(x, y) result(res)
    type(modint), intent(in) :: x(:, :), y(:)
    type(modint) :: res(size(x, 1))
    integer :: i
    do i = 1, size(x, 1)
      call asgn(res(i), dot_prod(x(i, :), y))
    end do
  end
  pure function matmul2(x, y) result(res)
    type(modint), intent(in) :: x(:), y(:, :)
    type(modint) :: res(size(y, 2))
    integer :: i
    do i = 1, size(y, 2)
      call asgn(res(i), dot_prod(x, y(:, i)))
    end do
  end
  pure function matmul3(x, y) result(res)
    type(modint), intent(in) :: x(:, :), y(:, :)
    type(modint) :: res(size(x, 1), size(y, 2))
    integer :: i
    do i = 1, size(x, 1)
      call asgn(res(i, :), matmul2(x(i, :), y))
    end do
  end
end module module_modint
module modulo_utils
  use module_modint
  implicit none
  type(modint), allocatable :: f(:), invf(:)
contains
  subroutine init(n)
    integer, intent(in) :: n
    integer :: i
    if (allocated(f)) deallocate(f)
    if (allocated(invf)) deallocate(invf)
    allocate(f(0:n), invf(0:n))
    f(0) = 1
    do i = 1, n
      f(i) = f(i - 1) * i
    end do
    invf(n) = f(n)%inv()
    do i = n, 1, -1
      invf(i - 1) = invf(i) * i
    end do
  end
  pure elemental type(modint) function perm(n, k) result(res)
    integer, intent(in) :: n, k
    if (k > n .or. n < 0 .or. k < 0) return
    res = f(n) * invf(n - k)
  end
  pure elemental type(modint) function comb(n, k) result(res)
    integer, intent(in) :: n, k
    res = perm(n, k) * invf(k)
  end
  pure elemental type(modint) function homo(n, k) result(res)
    integer, intent(in) :: n, k
    res = 1
    if (n == 0 .and. k == 0) return
    res = comb(n + k - 1, k)
  end
end module modulo_utils
program neq
  use module_modint
  use modulo_utils
  implicit none
  integer :: n, m, i
  type(modint) :: x, p
  read(*,*) n, m
  call init(m)
  do i = 0, n
    p = perm(m - i, n - i)
    x = x + comb(n, i) * merge(-1, 1, btest(i, 0)) * &
    perm(m, i) * p * p
  end do
  write(*,'(i0)') x%get()
end program neq

Measured memory_bytes: 10548