Fortran 95 source code to compute modular exponetiation recursively.

modular exponentiation_2.f95

!****************************************************
! べき剰余計算(recursive)
!
!   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: べき剰余
  recursive function calc_mod_exp(b, e, m) result(me)
    implicit none
    integer(SP), intent(in) :: b, e, m
    integer(SP) :: me

    me = 1
    if (e == 0) return
    me = calc_mod_exp(b, rshift(e, 1), m)
    me = mod(me * me, m)
    if (iand(e, 1) == 1) me = mod(me * b, m)
  end function calc_mod_exp
end program modulare_exponentiation