Fortran 95 source code to compute modular exponetiation iteratively.

modular exponentiation_1.f95

!****************************************************
! べき剰余計算(iterative)
!
!   date          name            version
!   2018.10.19    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
program modulare_exponentiation
  implicit none
  ! SP: 単精度(4), DP: 倍精度(8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
  integer(SP), parameter :: b = 12345
  integer(SP), parameter :: e = 6789
  integer(SP), parameter :: m = 4567

  print '(I0, "^", I0, " mod ", I0 " = ", I0)', &
      & b, e, m, calc_mod_exp(b, e, m)

  stop
contains
  ! べき剰余計算
  ! * 奇数判定にはビットごとの論理積を使用
  ! * 2での除算には1ビット右シフトを使用
  !
  ! :param(in) integer(4)  b: B
  ! :param(in) integer(4)  e: E
  ! :param(in) integer(4)  m: M
  ! :return    integer(4) me: べき剰余
  function calc_mod_exp(b, e, m) result(me)
    implicit none
    integer(SP), intent(in) :: b, e, m
    integer(SP) :: me
    integer(SP) :: b_w, e_w, m_w  ! 計算用

    b_w = b
    e_w = e
    m_w = m
    me = 1
    do while (e_w > 0)
      if (iand(e_w, 1) == 1) me = mod(me * b_w, m_w)
      e_w = rshift(e_w, 1)
      b_w = mod(b_w * b_w, m_w)
    end do
  end function calc_mod_exp
end program modulare_exponentiation