module n3long

   ! 64-BIT UNSIGNED INTEGER

   implicit none

   type long
     integer(4):: i(2)
   end type
   integer, parameter:: HI = 2, LO = 1

   integer, parameter:: INT4BITS = bit_size(0_4)

   interface assignment(=)
     module procedure long_from_int
     module procedure long_to_int
     module procedure long_from_double
     module procedure long_to_double
   end interface

   interface operator(+)
     module procedure long_add
     module procedure long_add_int
     module procedure long_add_double
   end interface

   interface operator(-)
     module procedure long_sub
   end interface

   interface operator(*)
     module procedure long_mul
   end interface

   interface operator(/)
     module procedure long_div
   end interface

   interface operator(.mod.)
     module procedure long_mod
   end interface

   interface operator(==)
     module procedure long_eq
     module procedure long_eq_int
   end interface

   interface operator(/=)
     module procedure long_ne
     module procedure long_ne_int
   end interface

   interface operator(<)
     module procedure long_lt
   end interface

   interface operator(<=)
     module procedure long_le
   end interface

   interface operator(>)
     module procedure long_gt
   end interface

   interface operator(>=)
     module procedure long_ge
   end interface

   interface int
     module procedure long_to_int_func
   end interface

   interface dble
     module procedure long_to_double_func
   end interface

   interface hex
     module procedure long_hex
   end interface

contains

  subroutine expand16(val, xpand)
    type(long), intent(in):: val
    integer, intent(out):: xpand(4)
    xpand(1) = ibits(val%i(LO),  0, 16)
    xpand(2) = ibits(val%i(LO), 16, 16)
    xpand(3) = ibits(val%i(HI),  0, 16)
    xpand(4) = ibits(val%i(HI), 16, 16)
  end subroutine

  subroutine unexpand16(xpand, out)
    integer, intent(in):: xpand(4)
    type(long), intent(out):: out
    integer:: val(4), sur(4)
    val(1) = ibits(xpand(1), 0, 16)
    sur(1) = ishft(xpand(1), -16)
    val(2) = ibits(xpand(2), 0, 16) + sur(1)
    sur(2) = ishft(xpand(2), -16) + ishft(val(2), -16)
    val(3) = ibits(xpand(3), 0, 16) + sur(2)
    sur(3) = ishft(xpand(3), -16) + ishft(val(3), -16)
    val(4) =       xpand(4)         + sur(3)
    out%i(LO) = ior(ibits(val(1), 0, 16), ishft(ibits(val(2), 0, 16), 16))
    out%i(HI) = ior(ibits(val(3), 0, 16), ishft(val(4), 16))
  end subroutine

  subroutine expand8(val, xpand)
    type(long), intent(in):: val
    integer, intent(out):: xpand(8)
    xpand(1) = ibits(val%i(LO),  0, 8)
    xpand(2) = ibits(val%i(LO),  8, 8)
    xpand(3) = ibits(val%i(LO), 16, 8)
    xpand(4) = ibits(val%i(LO), 24, 8)
    xpand(5) = ibits(val%i(HI),  0, 8)
    xpand(6) = ibits(val%i(HI),  8, 8)
    xpand(7) = ibits(val%i(HI), 16, 8)
    xpand(8) = ibits(val%i(HI), 24, 8)
  end subroutine

  subroutine unexpand8(xpand, out)
    integer, intent(in):: xpand(8)
    type(long), intent(out):: out
    integer:: val(8), sur(8), i
    val(1) = ibits(xpand(1), 0, 8)
    sur(1) = ishft(xpand(1), -8)
    do, i = 2, 7
      val(i) = ibits(xpand(i), 0, 8) + sur(i - 1)
      sur(i) = ishft(xpand(i), -8) + ishft(val(i), -8)
      val(i) = ibits(val(i), 0, 8)
    enddo
    val(8) = xpand(8) + sur(7)
    out%i(LO) = ior( &
      ior(val(1), ishft(val(2), 8)), &
      ior(ishft(val(3), 16), ishft(val(4), 24)))
    out%i(HI) = ior( &
      ior(val(5), ishft(val(6), 8)), &
      ior(ishft(val(7), 16), ishft(val(8), 24)))
  end subroutine

  subroutine long_by_muladd(x, i, j, k)
    type(long), intent(out):: x
    integer, intent(in):: i, j, k
    type(long):: y, z
    x = i
    y = j
    z = k
    x = x * y + z
  end subroutine

  ! NEGATIVE rhs IS TREATED AS MODULO
  !
  subroutine long_from_int(lhs, rhs)
    type(long), intent(out):: lhs
    integer, intent(in):: rhs
    if (bit_size(rhs) <= INT4BITS) then
      lhs%i(LO) = rhs
      lhs%i(HI) = 0
    else
      lhs%i(LO) = ibits(rhs, 0, INT4BITS)
      lhs%i(HI) = ibits(rhs, INT4BITS, bit_size(rhs) - INT4BITS)
    endif
  end subroutine

  ! CAST OPERATION MAY CAUSE CHANGE IN NUMBER
  !
  subroutine long_to_int(lhs, rhs)
    integer, intent(out):: lhs
    type(long), intent(in):: rhs
    lhs = rhs%i(LO)
  end subroutine

  integer function long_to_int_func(i) result(result)
    type(long), intent(in):: i
     call long_to_int(result, i)
  end function

  ! NEGATIVE INPUT: ABS VALUE
  !
  subroutine long_from_double(lhs, rhs)
    type(long), intent(out):: lhs
    double precision, intent(in):: rhs
    double precision:: x
    integer:: i, e(4)
    if (radix(rhs) /= 2) then
      stop 'radix of REAL must be two'
    endif
    x = abs(rhs)
    do, i = 4, 1, -1
      e(i) = ibits(floor(scale(x, -((i - 1) * 16))), 0, 16)
    enddo
    call unexpand16(e, lhs)
  end subroutine

  subroutine long_to_double(lhs, rhs)
    double precision, intent(out):: lhs
    type(long), intent(in):: rhs
    integer:: e(4), i
    call expand16(rhs, e)
    lhs = 0.0d0
    do, i = 1, 4
      lhs = lhs + scale(dble(e(i)), ((i-1)*16))
    enddo
  end subroutine

  double precision function long_to_double_func(i) result(result)
    type(long), intent(in):: i
     call long_to_double(result, i)
  end function

  type(long) function long_add(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    integer:: a(4), b(4), c(4)
    call expand16(lhs, a)
    call expand16(rhs, b)
    c(:) = a(:) + b(:)
    call unexpand16(c, result)
  end function

  type(long) function long_add_int(lhs, rhs) result(result)
    type(long), intent(in):: lhs
    integer, intent(in):: rhs
    type(long):: rlong
    rlong = rhs
    result = lhs + rlong
  end function

  type(long) function long_add_double(lhs, rhs) result(result)
    type(long), intent(in):: lhs
    double precision, intent(in):: rhs
    type(long):: rlong
    rlong = rhs
    result = lhs + rlong
  end function

  type(long) function long_sub(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    integer:: a(4), b(4), c(4), i
    call expand16(lhs, a)
    call expand16(rhs, b)
    c(:) = a(:) - b(:)
    do, i = 1, 3
      if (c(i) < 0) then
        c(i + 1) = c(i + 1) - 1
        c(i) = c(i) + (2 ** 16)
      endif
    enddo
    call unexpand16(c, result)
  end function

  type(long) function long_mul(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    integer:: a(8), b(8), c(8)
    integer:: i, j, k, r
    call expand8(lhs, a)
    call expand8(rhs, b)
    c(:) = 0
    do, i = 1, 8
      do, j = 1, (8 - (i - 1))
        k = i + j - 1
        c(k) = c(k) + a(i) * b(j)
      enddo
    enddo
    call unexpand8(c, result)
  end function

  type(long) function long_div(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    type(long):: q
    if (all(rhs%i == 0) .or. all(lhs%i == 0)) then
      result%i(:) = 0
      return
    endif
    if (lhs%i(HI) == 0 .and. lhs%i(LO) > 0 .and. &
    & rhs%i(HI) == 0 .and. rhs%i(LO) > 0) then
      result%i(HI) = 0
      result%i(LO) = lhs%i(LO) / rhs%i(LO)
      return
    endif
    result = DBLE(lhs) / DBLE(rhs)
    do
      q = lhs - rhs * result
      if (q < rhs) exit
      result = result + 1
    enddo
  end function

  type(long) function long_mod(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    type(long):: q
    if (rhs%i(HI) == 0 .and. rhs%i(LO) == 0) then
      result = lhs
    else if (lhs%i(HI) == 0 .and. rhs%i(HI) == 0 .and. rhs%i(LO) > 0) then
      result%i(HI) = 0
      result%i(LO) = mod(ibits(lhs%i(LO), 0, 31), rhs%i(LO)) + &
      & mod(mod(ibset(0, 31), rhs%i(LO)) * 2, rhs%i(LO))
    else
      result = lhs - rhs * long_div(lhs, rhs)
    endif
  end function

  subroutine long_shift(x, n)
    type(long), intent(inout):: x
    integer, intent(in):: n
    if (n <= -64) then
      x%i(:) = 0
    else if (n <= -32) then
      x%i(LO) = ishft(x%i(HI), n + 32)
      x%i(HI) = 0
    else if (n < 0) then
      x%i(LO) = ior(ishft(x%i(LO), n), ishft(x%i(HI), n + 32))
      x%i(HI) = ishft(x%i(HI), n)
    else if (n >= 64) then
      x%i(:) = 0
    else if (n >= 32) then
      x%i(Hi) = ishft(x%i(LO), n - 32)
      x%i(LO) = 0
    else if (n > 0) then
      x%i(HI) = ior(ishft(x%i(LO), n - 32), ishft(x%i(HI), n))
      x%i(LO) = ishft(x%i(LO), n)
    else
      continue
    endif
  end subroutine

  subroutine long_or(x, n)
    type(long), intent(inout):: x
    integer, intent(in):: n
    if (BIT_SIZE(n) > 32) then
      x%i(LO) = ior(x%i(LO), ibits(n, 0, 32))
      x%i(HI) = ior(x%i(HI), ishft(n, -32))
    else
      x%i(LO) = ior(x%i(LO), n)
    endif
  end subroutine

  logical function long_eq(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    result = all(lhs%i(:) == rhs%i(:))
  end function

  logical function long_eq_int(lhs, rhs) result(result)
    type(long), intent(in):: lhs
    integer, intent(in):: rhs
    type(long):: a
    a = rhs
    result = (lhs == a)
  end function

  logical function long_ne(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    result = any(lhs%i(:) /= rhs%i(:))
  end function

  logical function long_ne_int(lhs, rhs) result(result)
    type(long), intent(in):: lhs
    integer, intent(in):: rhs
    type(long):: a
    a = rhs
    result = (lhs /= a)
  end function

  logical function long_lt(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    integer:: a(4), b(4), i
    call expand16(lhs, a)
    call expand16(rhs, b)
    do, i = 4, 1, -1
      if (a(i) /= b(i)) then
        result = (a(i) .lt. b(i))
        return
      endif
    enddo
    result = .false.
  end function

  logical function long_le(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    integer:: a(4), b(4), i
    call expand16(lhs, a)
    call expand16(rhs, b)
    do, i = 4, 1, -1
      if (a(i) /= b(i)) then
        result = (a(i) .lt. b(i))
        return
      endif
    enddo
    result = .true.
  end function

  logical function long_gt(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    integer:: a(4), b(4), i
    call expand16(lhs, a)
    call expand16(rhs, b)
    do, i = 4, 1, -1
      if (a(i) /= b(i)) then
        result = (a(i) .gt. b(i))
        return
      endif
    enddo
    result = .false.
  end function

  logical function long_ge(lhs, rhs) result(result)
    type(long), intent(in):: lhs, rhs
    integer:: a(4), b(4), i
    call expand16(lhs, a)
    call expand16(rhs, b)
    do, i = 4, 1, -1
      if (a(i) /= b(i)) then
        result = (a(i) .gt. b(i))
        return
      endif
    enddo
    result = .true.
  end function

  character(16) function long_hex(li) result(result)
    use n3fmt, only: hex
    type(long), intent(in):: li
    integer:: ilo, ihi
    result(1:8) = hex(li%i(HI))
    result(9:16) = hex(li%i(LO))
  end function

  subroutine long_dfmt_hex(buf, x, stat)
    use n3fmt, only: hexdigit2val
    character(*), intent(in):: buf
    type(long), intent(out):: x
    integer, intent(out):: stat
    logical:: init
    integer:: digit
    x = 0
    stat = 0
    init = .true.
    do
      stat = stat + 1
      if (stat > len(buf)) exit
      if (init) then
        select case(buf(stat:stat))
        case('0':'9', 'A':'F', 'a':'f')
          init = .false.
        case(' ')
          cycle
        case default
          exit
        end select
      endif
      call hexdigit2val(buf(stat:stat), digit)
      if (digit == -1) exit
      call long_shift(x, 4)
      call long_or(x, digit)
    enddo
  end subroutine

end module