Skip to content

Instantly share code, notes, and snippets.

@komasaru
Created January 27, 2020 02:44
Show Gist options
  • Save komasaru/d4a44bc15bfe812fef73682cd14e52d8 to your computer and use it in GitHub Desktop.
Save komasaru/d4a44bc15bfe812fef73682cd14e52d8 to your computer and use it in GitHub Desktop.
Fortran 95 source code to calculate binomial coefficients.
!****************************************************
! 二項係数の計算
!
! 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment