module big_integers             ! Big (but not very big) integers.

  !
  ! A very primitive multiple precision module, using the classical
  ! algorithms (long multiplication, long division) and a mere 8-bit
  ! "digit" size. Some procedures might assume integers are in a
  ! two's-complement representation. This module is good enough for us
  ! to fulfill the task.
  !

  ! NOTE: I assume that iachar and achar do not alter the most
  !       significant bit.

  use, intrinsic :: iso_fortran_env, only: int16
  implicit none
  private

  public :: big_integer
  public :: integer2big
  public :: string2big
  public :: big2string
  public :: big_sgn
  public :: big_cmp, big_cmpabs
  public :: big_neg, big_abs
  public :: big_addabs, big_add
  public :: big_subabs, big_sub
  public :: big_mul      ! One might also include a big_muladd.
  public :: big_divrem   ! One could also include big_div and big_rem.
  public :: big_pow
  public :: operator(+)
  public :: operator(-)
  public :: operator(*)
  public :: operator(**)

  type :: big_integer
     ! The representation is sign-magnitude. The radix is 256, which
     ! is not speed-efficient, but which seemed relatively easy to
     ! work with if one were writing in standard Fortran (and assuming
     ! iachar and achar were "8-bit clean").
     logical :: sign = .false.  ! .false. => +sign, .true. => -sign.
     character, allocatable :: bytes(:)
  end type big_integer

  character, parameter :: zero = achar (0)
  character, parameter :: one = achar (1)

  ! An integer type capable of holding an unsigned 8-bit value.
  integer, parameter :: bytekind = int16

  interface operator(+)
     module procedure big_add
  end interface

  interface operator(-)
     module procedure big_neg
     module procedure big_sub
  end interface

  interface operator(*)
     module procedure big_mul
  end interface

  interface operator(**)
     module procedure big_pow
  end interface

contains

  elemental function logical2byte (bool) result (byte)
    logical, intent(in) :: bool
    character :: byte
    if (bool) then
       byte = one
    else
       byte = zero
    end if
  end function logical2byte

  elemental function logical2i (bool) result (i)
    logical, intent(in) :: bool
    integer :: i
    if (bool) then
       i = 1
    else
       i = 0
    end if
  end function logical2i

  elemental function byte2i (c) result (i)
    character, intent(in) :: c
    integer :: i
    i = iachar (c)
  end function byte2i

  elemental function i2byte (i) result (c)
    integer, intent(in) :: i
    character :: c
    c = achar (i)
  end function i2byte

  elemental function byte2bk (c) result (i)
    character, intent(in) :: c
    integer(bytekind) :: i
    i = iachar (c, kind = bytekind)
  end function byte2bk

  elemental function bk2byte (i) result (c)
    integer(bytekind), intent(in) :: i
    character :: c
    c = achar (i)
  end function bk2byte

  elemental function bk2i (i) result (j)
    integer(bytekind), intent(in) :: i
    integer :: j
    j = int (i)
  end function bk2i

  elemental function i2bk (i) result (j)
    integer, intent(in) :: i
    integer(bytekind) :: j
    j = int (iand (i, 255), kind = bytekind)
  end function i2bk

  ! Left shift of the least significant 8 bits of a bytekind integer.
  elemental function lshftbk (a, i) result (c)
    integer(bytekind), intent(in) :: a
    integer, intent(in) :: i
    integer(bytekind) :: c
    c = ishft (ibits (a, 0, 8 - i), i)
  end function lshftbk

  ! Right shift of the least significant 8 bits of a bytekind integer.
  elemental function rshftbk (a, i) result (c)
    integer(bytekind), intent(in) :: a
    integer, intent(in) :: i
    integer(bytekind) :: c
    c = ibits (a, i, 8 - i)
  end function rshftbk

  ! Left shift an integer.
  elemental function lshfti (a, i) result (c)
    integer, intent(in) :: a
    integer, intent(in) :: i
    integer :: c
    c = ishft (a, i)
  end function lshfti

  ! Right shift an integer.
  elemental function rshfti (a, i) result (c)
    integer, intent(in) :: a
    integer, intent(in) :: i
    integer :: c
    c = ishft (a, -i)
  end function rshfti

  function integer2big (i) result (a)
    integer, intent(in) :: i
    type(big_integer), allocatable :: a

    !
    ! To write a more efficient implementation of this procedure is
    ! left as an exercise for the reader.
    !

    character(len = 100) :: buffer

    write (buffer, '(I0)') i
    a = string2big (trim (buffer))
  end function integer2big

  function string2big (s) result (a)
    character(len = *), intent(in) :: s
    type(big_integer), allocatable :: a

    integer :: n, i, istart, iend
    integer :: digit

    if ((s(1:1) == '-') .or. s(1:1) == '+') then
       istart = 2
    else
       istart = 1
    end if

    iend = len (s)

    n = (iend - istart + 2) / 2

    allocate (a)
    allocate (a%bytes(n))

    a%bytes = zero
    do i = istart, iend
       digit = ichar (s(i:i)) - ichar ('0')
       if (digit < 0 .or. 9 < digit) error stop
       a = short_multiplication (a, 10)
       a = short_addition (a, digit)
    end do
    a%sign = (s(1:1) == '-')
    call normalize (a)
  end function string2big

  function big2string (a) result (s)
    type(big_integer), intent(in) :: a
    character(len = :), allocatable :: s

    type(big_integer), allocatable :: q
    integer :: r
    integer :: sgn

    sgn = big_sgn (a)
    if (sgn == 0) then
       s = '0'
    else
       q = a
       s = ''
       do while (big_sgn (q) /= 0)
          call short_division (q, 10, q, r)
          s = achar (r + ichar ('0')) // s
       end do
       if (sgn < 0) s = '-' // s
    end if
  end function big2string

  function big_sgn (a) result (sgn)
    type(big_integer), intent(in) :: a
    integer :: sgn

    integer :: n, i

    n = size (a%bytes)
    i = 1
    sgn = 1234
    do while (sgn == 1234)
       if (i == n + 1) then
          sgn = 0
       else if (a%bytes(i) /= zero) then
          if (a%sign) then
             sgn = -1
          else
             sgn = 1
          end if
       else
          i = i + 1
       end if
    end do
  end function big_sgn

  function big_cmp (a, b) result (cmp)
    type(big_integer(*)), intent(in) :: a, b
    integer :: cmp

    if (a%sign) then
       if (b%sign) then
          cmp = -big_cmpabs (a, b)
       else
          cmp = -1
       end if
    else
       if (b%sign) then
          cmp = 1
       else
          cmp = big_cmpabs (a, b)
       end if
    end if
  end function big_cmp

  function big_cmpabs (a, b) result (cmp)
    type(big_integer(*)), intent(in) :: a, b
    integer :: cmp

    integer :: n, i
    integer :: ia, ib

    cmp = 1234
    n = max (size (a%bytes), size (b%bytes))
    i = n
    do while (cmp == 1234)
       if (i == 0) then
          cmp = 0
       else
          ia = byteval (a, i)
          ib = byteval (b, i)
          if (ia < ib) then
             cmp = -1
          else if (ia > ib) then
             cmp = 1
          else
             i = i - 1
          end if
       end if
    end do
  end function big_cmpabs

  function big_neg (a) result (c)
    type(big_integer), intent(in) :: a
    type(big_integer), allocatable :: c
    c = a
    c%sign = .not. c%sign
  end function big_neg

  function big_abs (a) result (c)
    type(big_integer), intent(in) :: a
    type(big_integer), allocatable :: c
    c = a
    c%sign = .false.
  end function big_abs

  function big_add (a, b) result (c)
    type(big_integer), intent(in) :: a
    type(big_integer), intent(in) :: b
    type(big_integer), allocatable :: c

    logical :: sign

    if (a%sign) then
       if (b%sign) then      ! a <= 0, b <= 0
          c = big_addabs (a, b)
          sign = .true.
       else                  ! a <= 0, b >= 0
          c = big_subabs (a, b)
          sign = .not. c%sign
       end if
    else
       if (b%sign) then      ! a >= 0, b <= 0
          c = big_subabs (a, b)
          sign = c%sign
       else                  ! a >= 0, b >= 0
          c = big_addabs (a, b)
          sign = .false.
       end if
    end if
    c%sign = sign
  end function big_add

  function big_sub (a, b) result (c)
    type(big_integer), intent(in) :: a
    type(big_integer), intent(in) :: b
    type(big_integer), allocatable :: c

    logical :: sign

    if (a%sign) then
       if (b%sign) then      ! a <= 0, b <= 0
          c = big_subabs (a, b)
          sign = .not. c%sign
       else                  ! a <= 0, b >= 0
          c = big_addabs (a, b)
          sign = .true.
       end if
    else
       if (b%sign) then      ! a >= 0, b <= 0
          c = big_addabs (a, b)
          sign = .false.
       else                  ! a >= 0, b >= 0
          c = big_subabs (a, b)
          sign = c%sign
       end if
    end if
    c%sign = sign
  end function big_sub

  function big_addabs (a, b) result (c)
    type(big_integer), intent(in) :: a, b
    type(big_integer), allocatable :: c

    ! Compute abs(a) + abs(b).

    integer :: n, nc, i
    logical :: carry
    type(big_integer), allocatable :: tmp

    n = max (size (a%bytes), size (b%bytes))
    nc = n + 1

    allocate(tmp)
    allocate(tmp%bytes(nc))

    call add_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
    do i = 2, n
       call add_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
    end do
    tmp%bytes(nc) = logical2byte (carry)
    call normalize (tmp)
    c = tmp
  end function big_addabs

  function big_subabs (a, b) result (c)
    type(big_integer), intent(in) :: a, b
    type(big_integer), allocatable :: c

    ! Compute abs(a) - abs(b). The result is signed.

    integer :: n, i
    logical :: carry
    type(big_integer), allocatable :: tmp

    n = max (size (a%bytes), size (b%bytes))
    allocate(tmp)
    allocate(tmp%bytes(n))

    if (big_cmpabs (a, b) >= 0) then
       tmp%sign = .false.
       call sub_bytes (get_byte (a, 1), get_byte (b, 1), .false., tmp%bytes(1), carry)
       do i = 2, n
          call sub_bytes (get_byte (a, i), get_byte (b, i), carry, tmp%bytes(i), carry)
       end do
    else
       tmp%sign = .true.
       call sub_bytes (get_byte (b, 1), get_byte (a, 1), .false., tmp%bytes(1), carry)
       do i = 2, n
          call sub_bytes (get_byte (b, i), get_byte (a, i), carry, tmp%bytes(i), carry)
       end do
    end if
    call normalize (tmp)
    c = tmp
  end function big_subabs

  function big_mul (a, b) result (c)
    type(big_integer), intent(in) :: a, b
    type(big_integer), allocatable :: c

    !
    ! This is Knuth, Volume 2, Algorithm 4.3.1M.
    !

    integer :: na, nb, nc
    integer :: i, j
    integer :: ia, ib, ic
    integer :: carry
    type(big_integer), allocatable :: tmp

    na = size (a%bytes)
    nb = size (b%bytes)
    nc = na + nb + 1

    allocate (tmp)
    allocate (tmp%bytes(nc))

    tmp%bytes = zero
    j = 1
    do j = 1, nb
       ib = byte2i (b%bytes(j))
       if (ib /= 0) then
          carry = 0
          do i = 1, na
             ia = byte2i (a%bytes(i))
             ic = byte2i (tmp%bytes(i + j - 1))
             ic = (ia * ib) + ic + carry
             tmp%bytes(i + j - 1) = i2byte (iand (ic, 255))
             carry = ishft (ic, -8)
          end do
          tmp%bytes(na + j) = i2byte (carry)
       end if
    end do
    tmp%sign = (a%sign .neqv. b%sign)
    call normalize (tmp)
    c = tmp
  end function big_mul

  subroutine big_divrem (a, b, q, r)
    type(big_integer), intent(in) :: a, b
    type(big_integer), allocatable, intent(inout) :: q, r

    !
    ! Division with a remainder that is never negative. Equivalently,
    ! this is floor division if the divisor is positive, and ceiling
    ! division if the divisor is negative.
    !
    ! See Raymond T. Boute, "The Euclidean definition of the functions
    ! div and mod", ACM Transactions on Programming Languages and
    ! Systems, Volume 14, Issue 2, pp. 127-144.
    ! https://doi.org/10.1145/128861.128862
    !

    call nonnegative_division (a, b, .true., .true., q, r)
    if (a%sign) then
       if (big_sgn (r) /= 0) then
          q = short_addition (q, 1)
          r = big_sub (big_abs (b), r)
       end if
       q%sign = .not. b%sign
    else
       q%sign = b%sign
    end if
  end subroutine big_divrem

  function short_addition (a, b) result (c)
    type(big_integer), intent(in) :: a
    integer, intent(in) :: b
    type(big_integer), allocatable :: c

    ! Compute abs(a) + b.

    integer :: na, nc, i
    logical :: carry
    type(big_integer), allocatable :: tmp

    na = size (a%bytes)
    nc = na + 1

    allocate(tmp)
    allocate(tmp%bytes(nc))

    call add_bytes (a%bytes(1), i2byte (b), .false., tmp%bytes(1), carry)
    do i = 2, na
       call add_bytes (a%bytes(i), zero, carry, tmp%bytes(i), carry)
    end do
    tmp%bytes(nc) = logical2byte (carry)
    call normalize (tmp)
    c = tmp
  end function short_addition

  function short_multiplication (a, b) result (c)
    type(big_integer), intent(in) :: a
    integer, intent(in) :: b
    type(big_integer), allocatable :: c

    integer :: i, na, nc
    integer :: ia, ic
    integer :: carry
    type(big_integer), allocatable :: tmp

    na = size (a%bytes)
    nc = na + 1

    allocate (tmp)
    allocate (tmp%bytes(nc))

    tmp%sign = a%sign
    carry = 0
    do i = 1, na
       ia = byte2i (a%bytes(i))
       ic = (ia * b) + carry
       tmp%bytes(i) = i2byte (iand (ic, 255))
       carry = ishft (ic, -8)
    end do
    tmp%bytes(nc) = i2byte (carry)
    call normalize (tmp)
    c = tmp
  end function short_multiplication

  ! Division without regard to signs.
  subroutine nonnegative_division (a, b, want_q, want_r, q, r)
    type(big_integer), intent(in) :: a, b
    logical, intent(in) :: want_q, want_r
    type(big_integer), intent(inout), allocatable :: q, r

    integer :: na, nb
    integer :: remainder

    na = size (a%bytes)
    nb = size (b%bytes)

    ! It is an error if b has "significant" zero-bytes or is equal to
    ! zero.
    if (b%bytes(nb) == zero) error stop

    if (nb == 1) then
       if (want_q) then
          call short_division (a, byte2i (b%bytes(1)), q, remainder)
       else
          block
            type(big_integer), allocatable :: bit_bucket
            call short_division (a, byte2i (b%bytes(1)), bit_bucket, remainder)
          end block
       end if
       if (want_r) then
          if (allocated (r)) deallocate (r)
          allocate (r)
          allocate (r%bytes(1))
          r%bytes(1) = i2byte (remainder)
       end if
    else
       if (na >= nb) then
          call long_division (a, b, want_q, want_r, q, r)
       else
          if (want_q) q = string2big ("0")
          if (want_r) r = a
       end if
    end if
  end subroutine nonnegative_division

  subroutine short_division (a, b, q, r)
    type(big_integer), intent(in) :: a
    integer, intent(in) :: b
    type(big_integer), intent(inout), allocatable :: q
    integer, intent(inout) :: r

    !
    ! This is Knuth, Volume 2, Exercise 4.3.1.16.
    !
    ! The divisor is assumed to be positive.
    !

    integer :: n, i
    integer :: ia, ib, iq
    type(big_integer), allocatable :: tmp

    ib = b
    n = size (a%bytes)

    allocate (tmp)
    allocate (tmp%bytes(n))

    r = 0
    do i = n, 1, -1
       ia = (256 * r) + byte2i (a%bytes(i))
       iq = ia / ib
       r = mod (ia, ib)
       tmp%bytes(i) = i2byte (iq)
    end do
    tmp%sign = a%sign
    call normalize (tmp)
    q = tmp
  end subroutine short_division

  subroutine long_division (a, b, want_quotient, want_remainder, quotient, remainder)
    type(big_integer), intent(in) :: a, b
    logical, intent(in) :: want_quotient, want_remainder
    type(big_integer), intent(inout), allocatable :: quotient
    type(big_integer), intent(inout), allocatable :: remainder

    !
    ! This is Knuth, Volume 2, Algorithm 4.3.1D.
    !
    ! We do not deal here with the signs of the inputs and outputs.
    !
    ! It is assumed size(a%bytes) >= size(b%bytes), and that b has no
    ! leading zero-bytes and is at least two bytes long. If b is one
    ! byte long and nonzero, use short division.
    !

    integer :: na, nb, m, n
    integer :: num_lz, num_nonlz
    integer :: j
    integer :: qhat
    logical :: carry

    !
    ! We will NOT be working with VERY large numbers, and so it will
    ! be safe to put temporary storage on the stack. (Note: your
    ! Fortran might put this storage in a heap instead of the stack.)
    !
    !    v = b, normalized to put its most significant 1-bit all the
    !           way left.
    !
    !    u = a, shifted left by the same amount as b.
    !
    !    q = the quotient.
    !
    ! The remainder, although shifted left, will end up in u.
    !
    integer(bytekind) :: u(0:size (a%bytes) + size (b%bytes))
    integer(bytekind) :: v(0:size (b%bytes) - 1)
    integer(bytekind) :: q(0:size (a%bytes) - size (b%bytes))

    na = size (a%bytes)
    nb = size (b%bytes)

    n = nb
    m = na - nb

    ! In the most significant byte of the divisor, find the number of
    ! leading zero bits, and the number of bits after that.
    block
      integer(bytekind) :: tmp
      tmp = byte2bk (b%bytes(n))
      num_nonlz = bit_size (tmp) - leadz (tmp)
      num_lz = 8 - num_nonlz
    end block

    call normalize_v (b%bytes) ! Make the most significant bit of v be one.
    call normalize_u (a%bytes) ! Shifted by the same amount as v.

    ! Assure ourselves that the most significant bit of v is a one.
    if (.not. btest (v(n - 1), 7)) error stop

    do j = m, 0, -1
       call calculate_qhat (qhat)
       call multiply_and_subtract (carry)
       q(j) = i2bk (qhat)
       if (carry) call add_back
    end do

    if (want_quotient) then
       if (allocated (quotient)) deallocate (quotient)
       allocate (quotient)
       allocate (quotient%bytes(m + 1))
       quotient%bytes = bk2byte (q)
       call normalize (quotient)
    end if

    if (want_remainder) then
       if (allocated (remainder)) deallocate (remainder)
       allocate (remainder)
       allocate (remainder%bytes(n))
       call unnormalize_u (remainder%bytes)
       call normalize (remainder)
    end if

  contains

    subroutine normalize_v (b_bytes)
      character, intent(in) :: b_bytes(n)

      !
      ! Normalize v so its most significant bit is a one. Any
      ! normalization factor that achieves this goal will suffice; we
      ! choose 2**num_lz. (Knuth uses (2**32) div (y[n-1] + 1).)
      !
      ! Strictly for readability, we use linear stack space for an
      ! intermediate result.
      !

      integer :: i
      integer(bytekind) :: btmp(0:n - 1)

      btmp = byte2bk (b_bytes)

      v(0) = lshftbk (btmp(0), num_lz)
      do i = 1, n - 1
         v(i) = ior (lshftbk (btmp(i), num_lz), &
              &      rshftbk (btmp(i - 1), num_nonlz))
      end do
    end subroutine normalize_v

    subroutine normalize_u (a_bytes)
      character, intent(in) :: a_bytes(m + n)

      !
      ! Shift a leftwards to get u. Shift by as much as b was shifted
      ! to get v.
      !
      ! Strictly for readability, we use linear stack space for an
      ! intermediate result.
      !

      integer :: i
      integer(bytekind) :: atmp(0:m + n - 1)

      atmp = byte2bk (a_bytes)

      u(0) = lshftbk (atmp(0), num_lz)
      do i = 1, m + n - 1
         u(i) = ior (lshftbk (atmp(i), num_lz), &
              &      rshftbk (atmp(i - 1), num_nonlz))
      end do
      u(m + n) = rshftbk (atmp(m + n - 1), num_nonlz)
    end subroutine normalize_u

    subroutine unnormalize_u (r_bytes)
      character, intent(out) :: r_bytes(n)

      !
      ! Strictly for readability, we use linear stack space for an
      ! intermediate result.
      !

      integer :: i
      integer(bytekind) :: rtmp(0:n - 1)

      do i = 0, n - 1
         rtmp(i) = ior (rshftbk (u(i), num_lz), &
              &         lshftbk (u(i + 1), num_nonlz))
      end do
      rtmp(n - 1) = rshftbk (u(n - 1), num_lz)

      r_bytes = bk2byte (rtmp)
    end subroutine unnormalize_u

    subroutine calculate_qhat (qhat)
      integer, intent(out) :: qhat

      integer :: itmp, rhat
      logical :: adjust

      itmp = ior (lshfti (bk2i (u(j + n)), 8), &
           &      bk2i (u(j + n - 1)))
      qhat = itmp / bk2i (v(n - 1))
      rhat = mod (itmp, bk2i (v(n - 1)))
      adjust = .true.
      do while (adjust)
         if (rshfti (qhat, 8) /= 0) then
            continue
         else if (qhat * bk2i (v(n - 2)) &
              &     > ior (lshfti (rhat, 8), &
              &            bk2i (u(j + n - 2)))) then
            continue
         else
            adjust = .false.
         end if
         if (adjust) then
            qhat = qhat - 1
            rhat = rhat + bk2i (v(n - 1))
            if (rshfti (rhat, 8) == 0) then
               adjust = .false.
            end if
         end if
      end do
    end subroutine calculate_qhat

    subroutine multiply_and_subtract (carry)
      logical, intent(out) :: carry

      integer :: i
      integer :: qhat_v
      integer :: mul_carry, sub_carry
      integer :: diff

      mul_carry = 0
      sub_carry = 0
      do i = 0, n
         ! Multiplication.
         qhat_v = mul_carry
         if (i /= n) qhat_v = qhat_v + (qhat * bk2i (v(i)))
         mul_carry = rshfti (qhat_v, 8)
         qhat_v = iand (qhat_v, 255)

         ! Subtraction.
         diff = bk2i (u(j + i)) - qhat_v + sub_carry
         sub_carry = -(logical2i (diff < 0)) ! Carry 0 or -1.
         u(j + i) = i2bk (diff)
      end do
      carry = (sub_carry /= 0)
    end subroutine multiply_and_subtract

    subroutine add_back
      integer :: i, carry, sum

      q(j) = q(j) - 1_bytekind
      carry = 0
      do i = 0, n - 1
         sum = bk2i (u(j + i)) + bk2i (v(i)) + carry
         carry = ishft (sum, -8)
         u(j + i) = i2bk (sum)
      end do
    end subroutine add_back

  end subroutine long_division

  function big_pow (a, i) result (c)
    type(big_integer), intent(in) :: a
    integer, intent(in) :: i
    type(big_integer), allocatable :: c

    type(big_integer), allocatable :: base
    integer :: exponent, exponent_halved
    integer :: j, last_set

    if (i < 0) error stop

    if (i == 0) then
       c = integer2big (1)
    else
       last_set = bit_size (i) - leadz (i)
       base = a
       exponent = i
       c = integer2big (1)
       do j = 0, last_set - 1
          exponent_halved = exponent / 2
          if (2 * exponent_halved /= exponent) then
             c = c * base
          end if
          exponent = exponent_halved
          base = base * base
       end do
    end if
  end function big_pow

  subroutine add_bytes (a, b, carry_in, c, carry_out)
    character, intent(in) :: a, b
    logical, value :: carry_in
    character, intent(inout) :: c
    logical, intent(inout) :: carry_out

    integer :: ia, ib, ic

    ia = byte2i (a)
    if (carry_in) ia = ia + 1
    ib = byte2i (b)
    ic = ia + ib
    c = i2byte (iand (ic, 255))
    carry_out = (ic >= 256)
  end subroutine add_bytes

  subroutine sub_bytes (a, b, carry_in, c, carry_out)
    character, intent(in) :: a, b
    logical, value :: carry_in
    character, intent(inout) :: c
    logical, intent(inout) :: carry_out

    integer :: ia, ib, ic

    ia = byte2i (a)
    ib = byte2i (b)
    if (carry_in) ib = ib + 1
    ic = ia - ib
    carry_out = (ic < 0)
    if (carry_out) ic = ic + 256
    c = i2byte (iand (ic, 255))
  end subroutine sub_bytes

  function get_byte (a, i) result (byte)
    type(big_integer), intent(in) :: a
    integer, intent(in) :: i
    character :: byte

    if (size (a%bytes) < i) then
       byte = zero
    else
       byte = a%bytes(i)
    end if
  end function get_byte

  function byteval (a, i) result (v)
    type(big_integer), intent(in) :: a
    integer, intent(in) :: i
    integer :: v

    if (size (a%bytes) < i) then
       v = 0
    else
       v = byte2i (a%bytes(i))
    end if
  end function byteval

  subroutine normalize (a)
    type(big_integer), intent(inout) :: a

    logical :: done
    integer :: i
    character, allocatable :: fewer_bytes(:)

    ! Shorten to the minimum number of bytes.
    i = size (a%bytes)
    done = .false.
    do while (.not. done)
       if (i == 1) then
          done = .true.
       else if (a%bytes(i) /= zero) then
          done = .true.
       else
          i = i - 1
       end if
    end do
    if (i /= size (a%bytes)) then
       allocate (fewer_bytes (i))
       fewer_bytes = a%bytes(1:i)
       call move_alloc (fewer_bytes, a%bytes)
    end if

    ! If the magnitude is zero, then clear the sign bit.
    if (size (a%bytes) == 1) then
       if (a%bytes(1) == zero) then
          a%sign = .false.
       end if
    end if
  end subroutine normalize

end module big_integers

program modular_exponentiation_task
  use, non_intrinsic :: big_integers
  implicit none

  type(big_integer), allocatable :: zero, one, two
  type(big_integer), allocatable :: a, b, modulus

  zero = integer2big (0)
  one = integer2big (1)
  two = integer2big (2)

  a = string2big ("2988348162058574136915891421498819466320163312926952423791023078876139")
  b = string2big ("2351399303373464486466122544523690094744975233415544072992656881240319")
  modulus = string2big ("10") ** 40
  write (*,*) big2string (modular_pow (a, b, modulus))

contains

  ! The right-to-left binary method,
  ! https://en.wikipedia.org/w/index.php?title=Modular_exponentiation&oldid=1136216610#Right-to-left_binary_method
  function modular_pow (base, exponent, modulus) result (retval)
    type(big_integer), intent(in) :: base, exponent, modulus
    type(big_integer), allocatable :: retval

    type(big_integer), allocatable :: bas, expnt, remainder, bit_bucket

    if (big_sgn (modulus - one) == 0) then
       retval = zero
    else
       retval = one
       bas = base
       expnt = exponent
       do while (big_sgn (expnt) /= 0)
          call big_divrem (expnt, two, expnt, remainder)
          if (big_sgn (remainder) /= 0) then
             call big_divrem (retval * bas, modulus, bit_bucket, retval)
          end if
          call big_divrem (bas * bas, modulus, bit_bucket, bas)
       end do
    end if
  end function modular_pow

end program modular_exponentiation_task