Fortran 95 source code to calculate binomial coefficients.

binom_coeff.f95

!****************************************************
! 二項係数の計算
!
! DATE          AUTHOR          VERSION
! 2019.12.11    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2019 mk-mode.com All Rights Reserved.
!****************************************************
!
module cst
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
end module cst

module cmn
  use cst
  implicit none

  private
  public :: get_arg

contains
  ! コマンドライン引数(n, r)の取得
  !
  ! :param(out) integer(4) n
  ! :param(out) integer(4) r
  subroutine get_arg(n, r)
    implicit none
    integer(SP), intent(out) :: n, r
    character(8) :: a
    integer(SP)  :: ios

    n = 0
    r = 0
    if (iargc() == 2) then
      call getarg(1, a)
      read (a, *, iostat = ios) n
      if (ios /= 0) return
      call getarg(2, a)
      read (a, *, iostat = ios) r
      if (ios /= 0) return
    end if
  end subroutine get_arg

end module cmn

module comp
  use cst
  use FMZM
  implicit none
  private
  public :: binom_1, binom_2, binom_3, binom_4

contains
  ! 二項係数の計算(1)
  !   計算式: (n r) = n! / r!(n -r)!
  !
  ! :param(in) integer(4) n
  ! :param(in) integer(4) r
  ! :return    type(IM)   b
  function binom_1(n, r) result(b)
    use cst
    implicit none
    integer(SP), intent(in) :: n, r
    type(IM) :: b

    b = TO_IM('1')
    if (r == 0 .or. r == n) return
    b =  fact(n) / (fact(r) * fact(n - r))
  end function binom_1

  ! 二項係数の計算(2)
  !   計算式: (n r) = ((n - 1) r) + ((n - 1) (k - 1))
  !           (recursively)
  !
  ! :param(in) integer(4) n
  ! :param(in) integer(4) r
  ! :return    type(IM)   b
  recursive function binom_2(n, r) result(b)
    use cst
    implicit none
    integer(SP), intent(in) :: n, r
    type(IM) :: b

    b = TO_IM('1')
    if (r == 0 .or. r == n) return
    b = binom_2(n - 1, r) + binom_2(n - 1, r - 1)
  end function binom_2

  ! 二項係数の計算(3)
  !   計算式: (n r) = (n / r) * ((n - 1) (k - 1))
  !           (recursively)
  !
  ! :param(in) integer(4) n
  ! :param(in) integer(4) r
  ! :return    type(IM)   b
  recursive function binom_3(n, r) result(b)
    use cst
    implicit none
    integer(SP), intent(in) :: n, r
    type(IM) :: b

    b = TO_IM('1')
    if (r == 0 .or. r == n) return
    b = n * binom_3(n - 1, r - 1) / r
  end function binom_3

  ! 二項係数の計算(4)
  !   計算式: (n r) = Π(n - i + 1) / i  (i = 1, ..., r)
  !
  ! :param(in) integer(4) n
  ! :param(in) integer(4) r
  ! :return    type(IM)   b
  function binom_4(n, r) result(b)
    use cst
    implicit none
    integer(SP), intent(in) :: n, r
    integer(SP) :: i
    type(IM) :: b

    b = TO_IM('1')
    if (r == 0 .or. r == n) return
    do i = 1, r
      b = b * (n - i + 1) / i
    end do
  end function binom_4

  ! 階乗の計算
  !
  ! :param(in) integer(4)  x
  ! :return    type(IM)
  type(IM) function fact(x)
    use cst
    implicit none
    integer(SP), intent(in) :: x
    integer(SP) :: i

    fact = TO_IM('1')
    if (x == 0) return
    do i = 1, x
      fact = fact * i
    end do
  end function fact
end module comp

program binom_coeff
  use cst
  use cmn
  use comp
  use FMZM
  implicit none
  integer(SP) :: n, r
  type(IM)    :: b

  ! コマンドライン引数(n, r)取得
  call get_arg(n, r)
  if (n == 0 .and. r == 0) stop
  if (n < 0 .or. r < 0 .or. n < r) stop
  write (*, '("(", I5, " ", I5, ") = ")', advance='no') n, r

  ! 計算
  if (r * 2 > n) r = n - r
  b = binom_1(n, r)  ! binom_1 or 2 or 3 or 4
  call IMPRINT(b)
end program binom_coeff