       module unsigned_32
!
! Fortran 90 implementation of unsigned arithmetic - 32 bits
! (Assumes integers at least 32 bits long are available,
!  and are the default integer type)
! _____________________________________________________________________
! All rights to this code waived, so that it may be freely distributed
! as public domain software subject to the condition that these 5 lines
! are verbatim reproduced. Originally written by Michel Olagnon, from
! Ifremer, France, who would be pleased to receive your comments and
! corrections.
! _____________________________________________________________________
! M. Olagnon (Michel.Olagnon@ifremer.fr)
! version 1.3 of 6/30/94
!
! Main differences from version 1.1 : argument naming for intrinsics
!   is now consistent with standard.
!                                     elemental functions are provided
!   for rank 1 arrays.
!                                     missing intrinsics have been added
!
! _____________________________________________________________________
! Note: If you don't like code to start in column 7, remember that,
!       had Diophantes left a 6 characters margin, then mathematicians
!       might have spared much efforts on A**N = B**N + C**N ...
!       My margin is wide to let you put your corrections there.
! _____________________________________________________________________
!
       private
       public :: uint_32 
       public :: assignment (=)
       public :: operator (==)
       public :: operator (/=)
       public :: operator (>=)
       public :: operator (>)
       public :: operator (<=)
       public :: operator (<)
       public :: operator (+)
       public :: operator (-)
       public :: operator (*)
       public :: operator (/)
       public :: operator (**)
       public :: abs
       public :: bit_size
       public :: btest
       public :: char
       public :: cmplx
       public :: dble
       public :: digits
       public :: dim
       public :: dot_product
       public :: huge
       public :: iachar
       public :: iand
       public :: ibclr
       public :: ibits
       public :: ibset
       public :: ichar
       public :: ieor
       public :: int
       public :: ior
       public :: ishft
       public :: ishftc
       public :: matmul
       public :: max
       public :: maxloc
       public :: maxval
       public :: min
       public :: minloc
       public :: minval
       public :: mod
       public :: modulo
       public :: mvbits
       public :: not
       public :: product
       public :: radix
       public :: random_number
       public :: random_seed
       public :: range
       public :: real
       public :: sum
       public :: tiny
       public :: uabs
       public :: uachar
       public :: uchar
       public :: uint, unsigned
!
       type uint_32
          integer (selected_int_kind (9)) :: val
       end type uint_32
!
       character (len=10), parameter :: zdgt = "0123456789"
                                     !  decimal digits
       integer, parameter            :: lw = 32
                                     !  word length
       integer, parameter            :: nd = 10
                                     !  # of dec. digits in lw-bit word
       integer, parameter            :: lwm1 = (lw - 1)
       integer, parameter            :: lh = (lw / 2)
       integer, parameter            :: lrk = kind (1.0d0)
                                     ! largest available real kind
       real (selected_real_kind (nd, nd)), parameter :: one = 1.0
       real (selected_real_kind (nd, nd)), parameter :: two = 2.0
       real (selected_real_kind (nd, nd)), parameter :: x231 = two**lwm1
       real (selected_real_kind (nd, nd)), parameter :: x232 = two**lw
       type (uint_32), parameter     :: uten = uint_32 (10)
       type (uint_32), parameter     :: umul = uint_32 (843314861)
                                     ! random generator multiplier
       type (uint_32), parameter     :: uadd = uint_32 (453816693)
                                     ! random generator add. cnst.
       type (uint_32)                :: usee = uint_32 (19561003)
                                     ! random generator seed
!
       interface assignment (=)
          module procedure u32_itu, u32_i1tu, u32_i2tu,                  &
     &                     u32a_itu, u32a_i1tu, u32a_i2tu,               &
     &                     u32_uti, u32_uti1, u32_uti2,                  &
     &                     u32a_uti, u32a_uti1, u32a_uti2,               &
     &                     u32_rtu, u32_dtu, u32_utr, u32_utd,           &
     &                     u32a_rtu, u32a_dtu, u32a_utr, u32a_utd,       &
     &                     u32_ctu, u32_ktu, u32_utc, u32_utk,           &
     &                     u32a_ctu, u32a_ktu, u32a_utc, u32a_utk,       &
     &                     u32_stu, u32_uts,                             &
     &                     u32a_stu, u32a_uts
       end interface
!
       interface operator (==)
          module procedure u32_teq, u32a_teq
       end interface
!
       interface operator (/=)
          module procedure u32_tne, u32a_tne
       end interface
!
       interface operator (>=)
          module procedure u32_tge, u32a_tge
       end interface
!
       interface operator (>)
          module procedure u32_tgt, u32a_tgt
       end interface
!
       interface operator (<=)
          module procedure u32_tle, u32a_tle
       end interface
!
       interface operator (<)
          module procedure u32_tlt, u32a_tlt
       end interface
!
       interface operator (+)
          module procedure u32_add, u32a_add
       end interface
!
       interface operator (-)
          module procedure u32_sub, u32a_sub
       end interface
!
       interface operator (*)
          module procedure u32_mul, u32a_mul
       end interface
!
       interface operator (/)
          module procedure u32_div, u32a_div
       end interface
!
       interface operator (**)
          module procedure u32_exp, u32a_exp
       end interface
!
       interface abs
          module procedure u32_nop, u32a_nop
       end interface
!
       interface achar
          module procedure u32_ach, u32a_ach
       end interface
!
       interface bit_size
          module procedure u32_bsz
       end interface
!
       interface btest
          module procedure u32_tbt, u32a_tbt
       end interface
!
       interface char
          module procedure u32_chr, u32a_chr
       end interface
!
       interface cmplx
          module procedure u32_cpl, u32a_cpl
       end interface
!
       interface dble
          module procedure u32_dbl, u32a_dbl
       end interface
!
       interface digits
          module procedure u32_dgt
       end interface
!
       interface dim
          module procedure u32_dim, u32a_dim
       end interface
!
       interface dot_product
          module procedure u32_dot
       end interface
!
       interface huge
          module procedure u32_hge
       end interface
!
       interface iachar
          module procedure u32_uac, u32a_uac
       end interface
!
       interface iand
          module procedure u32_and, u32a_and
       end interface
!
       interface ibclr
          module procedure u32_bcl, u32a_bcl
       end interface
!
       interface ibits
          module procedure u32_bts, u32a_bts
       end interface
!
       interface ibset
          module procedure u32_bst, u32a_bst
       end interface
!
       interface ichar
          module procedure u32_uch, u32a_uch
       end interface
!
       interface ieor
          module procedure u32_xor, u32a_xor
       end interface
!
       interface int
          module procedure u32_int, u32a_int
       end interface
!
       interface ishft
          module procedure u32_sft, u32a_sft
       end interface
!
       interface ishftc
          module procedure u32_sfc, u32a_sfc
       end interface
!
       interface matmul
          module procedure u32_mmm, u32_vmm, u32_mmv
       end interface
!
       interface max
          module procedure u32_max, u32a_max
       end interface
!
       interface maxloc
          module procedure u32_mal
       end interface
!
       interface maxval
          module procedure u32_mav
       end interface
!
       interface min
          module procedure u32_min, u32a_min
       end interface
!
       interface minloc
          module procedure u32_mil
       end interface
!
       interface minval
          module procedure u32_miv
       end interface
!
       interface mod
          module procedure u32_mod, u32a_mod
       end interface
!
       interface modulo
          module procedure u32_mod, u32a_mod
       end interface
!
       interface mvbits
          module procedure u32_mvb, u32a_mvb
       end interface
!
       interface not
          module procedure u32_not, u32a_not
       end interface
!
       interface product
          module procedure u32_prd
       end interface
!
       interface radix
          module procedure u32_rdx
       end interface
!
       interface random_number
          module procedure u32_rnd, u32a_rnd
       end interface
!
       interface random_seed
          module procedure u32_see
       end interface
!
       interface range
          module procedure u32_rng
       end interface
!
       interface real
          module procedure u32_rea, u32a_rea
       end interface
!
       interface sum
          module procedure u32_sum
       end interface
!
       interface tiny
          module procedure u32_tny
       end interface
!
       interface uabs
          module procedure u32_nop, u32a_nop
       end interface
!
       interface uachar
          module procedure u32_uac, u32a_uac
       end interface
!
       interface uchar
          module procedure u32_uch, u32a_uch
       end interface
!
       interface uint
          module procedure u32_unsi1, u32_unsi2, u32_unsi,               &
     &                     u32a_unsi1, u32a_unsi2, u32a_unsi,            &
     &                     u32_unsr, u32_unsd,                           &
     &                     u32a_unsr, u32a_unsd,                         &
     &                     u32_unsc, u32_unsk,                           &
     &                     u32a_unsc, u32a_unsk,                         &
     &                     u32_unss,                                     &
     &                     u32a_unss
       end interface
!
       interface unsigned
          module procedure u32_unsi1, u32_unsi2, u32_unsi,               &
     &                     u32a_unsi1, u32a_unsi2, u32a_unsi,            &
     &                     u32_unsr, u32_unsd,                           &
     &                     u32a_unsr, u32a_unsd,                         &
     &                     u32_unsc, u32_unsk,                           &
     &                     u32a_unsc, u32a_unsk,                         &
     &                     u32_unss,                                     &
     &                     u32a_unss
       end interface
!
! __________________________________________________________________
!
       contains
!
! Section 1. Scalar operations
!
! _______nothing to do______________________________________________
! _____________________________________________________________________
!
         function u32_nop (a)
            type (uint_32), intent (in)    :: a
            type (uint_32)                 :: u32_nop
            u32_nop = a
         end function u32_nop
!
! _______integer to unsigned assignment_____________________________
! _____________________________________________________________________
!
         subroutine u32_itu (u, i)
            type (uint_32), intent (out)                 :: u
            integer, intent (in)                         :: i
            intrinsic ibits, bit_size, int, min
!
            u%val = ibits (i, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_itu
!
         subroutine u32_i1tu (u, i)
            type (uint_32), intent (out)                 :: u
            integer (selected_int_kind (2)), intent (in) :: i
            intrinsic ibits, bit_size, int, min
!
            u%val = ibits (i, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_i1tu
!
         subroutine u32_i2tu (u, i)
            type (uint_32), intent (out)                 :: u
            integer (selected_int_kind (4)), intent (in) :: i
            intrinsic ibits, bit_size, int, min
!
            u%val = ibits (i, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_i2tu
!
         subroutine u32_i4tu (u, i)
            type (uint_32), intent (out)                 :: u
            integer (selected_int_kind (9)), intent (in) :: i
            intrinsic ibits, bit_size, int, min
!
            u%val = ibits (i, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_i4tu
!
! _______explicit integer to unsigned conversions______________________
!
         function u32_unsi1 (i)
            integer (selected_int_kind (2)), intent (in) :: i
            type (uint_32)                               :: u32_unsi1
            intrinsic ibits, bit_size, int, min
!
            u32_unsi1%val = ibits (i, 0, min (lw, int (bit_size (i),     &
     &                                         kind = kind (lw))))
         end function u32_unsi1
!
         function u32_unsi2 (i)
            integer (selected_int_kind (4)), intent (in) :: i
            type (uint_32)                               :: u32_unsi2
            intrinsic ibits, bit_size, int, min
!
            u32_unsi2%val = ibits (i, 0, min (lw, int (bit_size (i),     &
     &                                         kind = kind (lw))))
         end function u32_unsi2
!
         function u32_unsi4 (i)
            integer (selected_int_kind (9)), intent (in) :: i
            type (uint_32)                               :: u32_unsi4
            intrinsic ibits, bit_size, int, min
!
            u32_unsi4%val = ibits (i, 0, min (lw, int (bit_size (i),     &
     &                                         kind = kind (lw))))
         end function u32_unsi4
!
         function u32_unsi (i)
            integer, intent (in) :: i
            type (uint_32)       :: u32_unsi
            intrinsic ibits, bit_size, int, min
!
            u32_unsi%val = ibits (i, 0, min (lw, int (bit_size (i),      &
     &                                         kind = kind (lw))))
         end function u32_unsi
!
! _______unsigned to integer assignment________________________________
!
         subroutine u32_uti (i, u)
            integer, intent (out)                         :: i
            type (uint_32), intent (in)                   :: u
            intrinsic ibits, bit_size, int, min
!
            i = ibits (u%val, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_uti
!
         subroutine u32_uti1 (i, u)
            integer (selected_int_kind (2)), intent (out) :: i
            type (uint_32), intent (in)                   :: u
            intrinsic ibits, bit_size, int, min
!
            i = ibits (u%val, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_uti1
!
         subroutine u32_uti2 (i, u)
            integer (selected_int_kind (4)), intent (out) :: i
            type (uint_32), intent (in)                   :: u
            intrinsic ibits, bit_size, int, min
!
            i = ibits (u%val, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_uti2
!
         subroutine u32_uti4 (i, u)
            integer (selected_int_kind (9)), intent (out) :: i
            type (uint_32), intent (in)                   :: u
            intrinsic ibits, bit_size, int, min
!
            i = ibits (u%val, 0, min (lw, int (bit_size (i),             &
     &                                         kind = kind (lw))))
         end subroutine u32_uti4
!
! _______explicit unsigned to integer conversion________________________
!
         integer function u32_int (u)
            type (uint_32), intent (in)                   :: u
            intrinsic ibits, bit_size, int, min
!
            u32_int = ibits (u%val, 0, min (lw, int (bit_size (i),       &
     &                                         kind = kind (lw))))
         end function u32_int
!
! _______real to unsigned assignment________________________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32_rtu (u, r)
            type (uint_32), intent (out)                 :: u
            real, intent (in)                            :: r
            integer            :: jwrk
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, ibset, floor, modulo, ibits, min, bit_size
!
            if (abs (r) < one) then
               u%val = 0
            else
               xwrk = r
               xwrk = modulo (xwrk, x232)
               if (xwrk == x231) then
                  jwrk = ibset (0, lwm1)
               elseif (xwrk > x231) then
                  xwrk = xwrk - x232
                  jwrk = floor (xwrk)
               else
                  jwrk = floor (xwrk)
               endif
               u%val = ibits (jwrk, 0, min (lw, bit_size (jwrk)))
            endif
         end subroutine u32_rtu
!
! _______explicit real to unsigned conversion__________________________
!
         function u32_unsr (r)
            real, intent (in) :: r
            type (uint_32)    :: u32_unsr
            type (uint_32)    :: uwrk
!
            call u32_rtu (uwrk, r)
            u32_unsr = uwrk
         end function u32_unsr
!
! _______unsigned to real assignment___________________________________
!
         subroutine u32_utr (r, u)
            real, intent (out)                          :: r
            type (uint_32), intent (in)                 :: u
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic real, ibits
!
            xwrk = real (ibits (u%val, 0, lwm1), kind = kind (xwrk))
            if (u32_tbt (u, lwm1)) then
               xwrk = x231 + xwrk
            endif
            r = xwrk
         end subroutine u32_utr
!
! _______explicit unsigned to real conversion_______________________
!
         function u32_rea (a)
            real                                        :: u32_rea
            type (uint_32), intent (in)                 :: a
            real                                        :: r
!
            call u32_utr (r, a)
            u32_rea = r
         end function u32_rea
!
! _______double to unsigned assignment______________________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32_dtu (u, r)
            type (uint_32), intent (out)                 :: u
            real (kind (1.0d0)), intent (in)             :: r
            integer            :: jwrk
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, floor, ibits, ibset, modulo, min, bit_size
!
            if (abs (r) < one) then
               u%val = 0
            else
               xwrk = r
               xwrk = modulo (xwrk, x232)
               if (xwrk == x231) then
                  jwrk = ibset (0, lwm1)
               elseif (xwrk > x231) then
                  xwrk = xwrk - x232
                  jwrk = floor (xwrk)
               else
                  jwrk = floor (xwrk)
               endif
               u%val = ibits (jwrk, 0, min (lw, bit_size (jwrk)))
            endif
         end subroutine u32_dtu
!
! _______explicit double to unsigned conversion________________________
!
         function u32_unsd (r)
            real (kind = kind (1.0d0)), intent (in) :: r
            type (uint_32)                          :: u32_unsd
            type (uint_32)                          :: uwrk
!
            call u32_dtu (uwrk, r)
            u32_unsd = uwrk
         end function u32_unsd
!
! _______unsigned to double assignment______________________________
!
         subroutine u32_utd (r, u)
            real (kind (1.0d0)), intent (out)           :: r
            type (uint_32), intent (in)                 :: u
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic real, ibits
!
            xwrk = real (ibits (u%val, 0, lwm1), kind = kind (xwrk))
            if (u32_tbt (u, lwm1)) then
               xwrk = x231 + xwrk
            endif
            r = xwrk
         end subroutine u32_utd
!
! _______explicit unsigned to double conversion_____________________
!
         function u32_dbl (a)
            real (kind (1.0d0))                         :: u32_dbl
            type (uint_32), intent (in)                 :: a
            real (kind (1.0d0))                         :: r
!
            call u32_utd (r, a)
            u32_dbl = r
         end function u32_dbl
!
! _______complex to unsigned assignment_____________________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32_ctu (u, c)
            type (uint_32), intent (out)                 :: u
            complex, intent (in)                         :: c
            real                                 :: r
            integer            :: jwrk
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, real, modulo, ibset, ibits, floor, min,       &
     &                bit_size
!
            r = real (c, kind = kind (r))
            if (abs (r) < one) then
               u%val = 0
            else
               xwrk = r
               xwrk = modulo (xwrk, x232)
               if (xwrk == x231) then
                  jwrk = ibset (0, lwm1)
               elseif (xwrk > x231) then
                  xwrk = xwrk - x232
                  jwrk = floor (xwrk)
               else
                  jwrk = floor (xwrk)
               endif
               u%val = ibits (jwrk, 0, min (lw, bit_size (jwrk)))
            endif
         end subroutine u32_ctu
!
! _______explicit complex to unsigned conversion_______________________
!
         function u32_unsc (c)
            complex, intent (in) :: c
            type (uint_32)       :: u32_unsc
            type (uint_32)       :: uwrk
!
            call u32_ctu (uwrk, c)
            u32_unsc = uwrk
         end function u32_unsc
!
! _______unsigned to complex assignment______________________________
!
         subroutine u32_utc (c, u)
            complex, intent (out)                       :: c
            type (uint_32), intent (in)                 :: u
            intrinsic cmplx, ibits, real
            real (selected_real_kind (nd, nd))   :: xwrk
!
            xwrk = real (ibits (u%val, 0, lwm1), kind = kind (xwrk))
            if (u32_tbt (u, lwm1)) then
               xwrk = x231 + xwrk
            endif
            c = cmplx (xwrk, 0.0, kind = kind (c))
         end subroutine u32_utc
!
! _______double complex to unsigned assignment______________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32_ktu (u, c)
            type (uint_32), intent (out)                 :: u
            complex (kind (1.0d0)), intent (in)          :: c
            real (kind (1.0d0))                  :: r
            integer            :: jwrk
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, floor, ibits, ibset, real, modulo, min,       &
     &                bit_size
!
            r = real (c, kind = kind (r))
            if (abs (r) < one) then
               u%val = 0
            else
               xwrk = r
               xwrk = modulo (xwrk, x232)
               if (xwrk == x231) then
                  jwrk = ibset (0, lwm1)
               elseif (xwrk > x231) then
                  xwrk = xwrk - x232
                  jwrk = floor (xwrk)
               else
                  jwrk = floor (xwrk)
               endif
               u%val = ibits (jwrk, 0, min (lw, bit_size (jwrk)))
            endif
         end subroutine u32_ktu
!
! _______explicit double complex to unsigned conversion________________
!
         function u32_unsk (c)
            complex (kind = kind (1.0d0)), intent (in) :: c
            type (uint_32)                             :: u32_unsk
            type (uint_32)                             :: uwrk
!
            call u32_ktu (uwrk, c)
            u32_unsk = uwrk
         end function u32_unsk
!
! _______unsigned to double complex assignment_________________________
!
         subroutine u32_utk (c, u)
            complex (kind (1.0d0)), intent (out)        :: c
            type (uint_32), intent (in)                 :: u
            intrinsic cmplx, ibits, real
            real (selected_real_kind (nd, nd))   :: xwrk
!
            xwrk = real (ibits (u%val, 0, lwm1), kind = kind (xwrk))
            if (u32_tbt (u, lwm1)) then
               xwrk = x231 + xwrk
            endif
            c = cmplx (xwrk, 0.0, kind = kind (c))
         end subroutine u32_utk
!
! _______explicit unsigned to complex conversions______________________
!
         function u32_cpl (x, y)
            complex                                :: u32_cpl
            type (uint_32), intent (in)            :: x
            type (uint_32), intent (in), optional  :: y
            intrinsic cmplx
            real (kind = lrk)                      :: rx, ry
!
            rx = u32_dbl (x)
            if (present (y)) then
               ry = u32_dbl (y)
            else
               ry = 0.0
            endif
            u32_cpl = cmplx (rx, ry)
         end function u32_cpl
!
! _______convert string to unsigned_____________________________________
!
         subroutine u32_stu (u, s)
            type (uint_32), intent (out)                 :: u
            character (len=*), intent (in)               :: s
            type (uint_32)     :: uwrk1, uwrk2
            integer            :: i, i0, i1, j
            intrinsic len_trim, index
!
            i0 = 1
            i1 = len_trim (s)
            uwrk1 = uint_32 (0)
            do i = i0, i1
               if (s (i:i) /= ' ' .and. s (i:i) /= '+') then
                  i0 = i
                  exit
               endif
            enddo
            do i = i0, i1
               j = index (zdgt, s (i:i)) - 1
               if (j < 0) then
                  exit
               else
                  uwrk1 = u32_mul (uwrk1, uten)
                  if (j > 0) then
                     uwrk2%val = j
                     uwrk1 = u32_add (uwrk1, uwrk2)
                  endif
               endif
            enddo
            u = uwrk1   
         end subroutine u32_stu
!
! _______explicit string to unsigned conversion________________________
!
         function u32_unss (s)
            character (len = *), intent (in) :: s
            type (uint_32)                   :: u32_unss
            type (uint_32)                   :: uwrk
!
            call u32_stu (uwrk, s)
            u32_unss = uwrk
         end function u32_unss
!
! _______convert unsigned to string____________________________________
!
         subroutine u32_uts (s, u)
            type (uint_32), intent (in)                  :: u
            character (len=*), intent (out)              :: s
            type (uint_32)     :: uwrk, uwrk1, uwrk2
            integer            :: i, i0, i1, j
            intrinsic repeat, len
!
            i0 = 1
            i1 = nd
            uwrk2 = uten
            do i = i0, i1
               if (u32_tlt (u, uwrk2)) then
                  i1 = i
                  exit
               endif
               uwrk2 = u32_mul (uten, uwrk2)
            enddo
            j = len (s)
            if (i1 > j) then
               s = repeat ('*', j)
            else
               s = ' '
               uwrk1 = u
               uwrk2 = u32_exp (uten, i1 - 1)
               i0 = j - i1 + 1
               i1 = j - 1
               do i = i0, i1
                  uwrk = u32_div (uwrk1, uwrk2)
                  s (i:i) = zdgt (uwrk%val+1:uwrk%val+1)
                  uwrk = u32_mul (uwrk, uwrk2)
                  uwrk1 = u32_sub (uwrk1, uwrk)
                  uwrk2 = u32_div (uwrk2, uten)
               enddo
               s (j:j) = zdgt (uwrk1%val+1:uwrk1%val+1)
            endif
         end subroutine u32_uts
!
! _______convert unsigned to char______________________________________
!
         function u32_chr (i)
            type (uint_32), intent (in)                  :: i
            character (len=1)                            :: u32_chr
            intrinsic  char
!
            u32_chr = char (i%val)
         end function u32_chr
!
! _______convert char to unsigned______________________________________
!
         function u32_uch (c)
            character (len=1), intent (in)               :: c
            type (uint_32)                               :: u32_uch
            intrinsic ichar
!
            u32_uch%val = ichar (c)
         end function u32_uch
!
! _______convert unsigned to ascii_____________________________________
!
         function u32_ach (i)
            type (uint_32), intent (in)                  :: i
            character (len=1)                            :: u32_ach
            intrinsic achar
!
            u32_ach = achar (i%val)
         end function u32_ach
!
! _______convert ascii to unsigned_____________________________________
!
         function u32_uac (c)
            character (len=1), intent (in)               :: c
            type (uint_32)                               :: u32_uac
            intrinsic iachar
!
            u32_uac%val = iachar (c)
         end function u32_uac
!
! _______test for equality______________________________________________
!
         logical function u32_teq (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            intrinsic ibits
!
         u32_teq = (ibits (u32a%val, 0, lh) == ibits (u32b%val,  0, lh))&
     &      .and. (ibits (u32a%val, lh, lh) == ibits (u32b%val, lh, lh))
         end function u32_teq
!
! _______test for inequality____________________________________________
!
         logical function u32_tne (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            intrinsic ibits
!
         u32_tne = (ibits (u32a%val, 0, lh) /= ibits (u32b%val,  0, lh))&
     &       .or. (ibits (u32a%val, lh, lh) /= ibits (u32b%val, lh, lh))
         end function u32_tne
!
! _______test for greater than__________________________________________
!
         logical function u32_tgt (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            integer            :: jwrk1, jwrk2
            intrinsic ibits
!
            jwrk1 = ibits (u32a%val, lh, lh)
            jwrk2 = ibits (u32b%val, lh, lh)
            if (jwrk1 == jwrk2) then
               u32_tgt = ibits (u32a%val, 0, lh) >                      &
     &                                          ibits (u32b%val,  0, lh) 
            else
               u32_tgt = jwrk1 > jwrk2
            endif
         end function u32_tgt
!
! _______test for lower than____________________________________________
!
         logical function u32_tlt (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            integer            :: jwrk1, jwrk2
            intrinsic ibits
!
            jwrk1 = ibits (u32a%val, lh, lh)
            jwrk2 = ibits (u32b%val, lh, lh)
            if (jwrk1 == jwrk2) then
               u32_tlt = ibits (u32a%val, 0, lh) <                      &
     &                                          ibits (u32b%val,  0, lh) 
            else
               u32_tlt = jwrk1 < jwrk2
            endif
         end function u32_tlt
!
! _______test for greater or equal_____________________________________
!
         logical function u32_tge (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            integer            :: jwrk1, jwrk2
            intrinsic ibits
!
            jwrk1 = ibits (u32a%val, lh, lh)
            jwrk2 = ibits (u32b%val, lh, lh)
            if (jwrk1 == jwrk2) then
               u32_tge = ibits (u32a%val, 0, lh) >=                     &
     &                                          ibits (u32b%val,  0, lh) 
            else
               u32_tge = jwrk1 > jwrk2
            endif
         end function u32_tge
!
! _______test for lower or equal_______________________________________
!
         logical function u32_tle (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            integer            :: jwrk1, jwrk2
            intrinsic ibits
!
            jwrk1 = ibits (u32a%val, lh, lh)
            jwrk2 = ibits (u32b%val, lh, lh)
            if (jwrk1 == jwrk2) then
               u32_tle = ibits (u32a%val, 0, lh) <=                     &
     &                                          ibits (u32b%val,  0, lh) 
            else
               u32_tle = jwrk1 < jwrk2
            endif
         end function u32_tle
!
! _______addition of two unsigned numbers_________________________
!
         function u32_add (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            type (uint_32)                               :: u32_add
            integer            :: jwrk1, jwrk2
            intrinsic ibits, ishft
!
            jwrk1 = ibits (u32a%val,  0, lh) + ibits (u32b%val,  0, lh)
            jwrk2 = ibits (u32a%val, lh, lh) + ibits (u32b%val, lh, lh) &
     &              + ibits (jwrk1, lh, lh)
            jwrk1 = ibits (jwrk1,  0, lh)
            jwrk2 = ishft (jwrk2, lh)
            u32_add%val = ior (jwrk1, jwrk2)
         end function u32_add
!
! _______substraction of two unsigned numbers____________________
!
         function u32_sub (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            type (uint_32)                               :: u32_sub
            integer            :: jwrk1, jwrk2, jwrkb
            intrinsic ibits, ishft, not, ior
!
            jwrkb = not (ibits (u32b%val,  0, lw))
            jwrk1 = ibits (u32a%val,  0, lh) + ibits (jwrkb,  0, lh) + 1
            jwrk2 = ibits (u32a%val, lh, lh) + ibits (jwrkb, lh, lh) +  &
     &              ibits (jwrk1, lh, lh)
            jwrk1 = ibits (jwrk1,  0, lh)
            jwrk2 = ishft (jwrk2, lh)
            u32_sub%val = ior (jwrk1, jwrk2)
         end function u32_sub
!
! _______multiplication of two unsigned numbers_________________________
!
         function u32_mul (u32a, u32b)
            type (uint_32), intent (in)                  :: u32a, u32b
            type (uint_32)                               :: u32_mul
            integer            :: jwrka1, jwrka2, jwrkb1, jwrkb2
            integer            :: jwrk1, jwrk2
            intrinsic ibits, ishft, ior
!
            jwrka1 = ibits (u32a%val,  0, lh)
            jwrkb1 = ibits (u32b%val,  0, lh)
            jwrka2 = ibits (u32a%val, lh, lh)
            jwrkb2 = ibits (u32b%val, lh, lh)
            jwrk1 = jwrka1 * jwrkb1
            jwrk2 = jwrka1 * jwrkb2 + jwrka2 * jwrkb1 +                 &
     &              ibits (jwrk1, lh, lh)
            jwrk1 = ibits (jwrk1,  0, lh)
            jwrk2 = ishft (jwrk2, lh)
            u32_mul%val = ior (jwrk1, jwrk2)
         end function u32_mul
!
! _______division of two unsigned numbers_______________________________
!
         function u32_div (u32a, u32b)
            type (uint_32), intent (in) :: u32a, u32b
            type (uint_32)              :: u32_div
            type (uint_32) :: uwrk1, uwrk2, uwrk
            intrinsic btest, ibits
!
            uwrk = uint_32 (0)
            if (u32_teq (u32b, uwrk)) then
               u32_div = uint_32 (0) ! Should raise exception, but ...
            elseif (u32_tlt (u32a, u32b)) then
               u32_div = uint_32 (0)
            else
               if (btest (u32a%val, lwm1)) then
!
! ... a large, b large
                  if (btest (u32b%val, lwm1)) then
                     u32_div = uint_32 (1)
!
! ... a large, b small, first divide a/2
                  else
                     uwrk1%val = ibits (u32a%val, 1, lwm1) /            &
     &                                       ibits (u32b%val, 0, lwm1)
                     uwrk2 = uint_32 (2)
                     uwrk2 = u32_mul (uwrk1, uwrk2)
                     uwrk1 = u32_mul (u32b, uwrk2)
                     uwrk  = u32_sub (u32a,  uwrk1)
                     if (u32_tlt (uwrk, u32b)) then
                        u32_div = uwrk2
                     else
                        uwrk = u32_sub (uwrk, u32b)
                        uwrk1%val = 1 + ibits (uwrk%val, 0, lwm1) /     &
     &                                       ibits (u32b%val, 0, lwm1)
                        u32_div = u32_add (uwrk1, uwrk2)
                     endif
                  endif
!
! ... a small, b small
               else
                  u32_div%val = ibits (u32a%val, 0, lwm1) /             &
     &                                       ibits (u32b%val, 0, lwm1)
               endif
            endif
         end function u32_div
!
! _______exponentiation to integer power_______________________________
!
         function u32_exp (u32a, jpow)
            type (uint_32), intent (in) :: u32a
            integer, intent (in)        :: jpow
            type (uint_32)              :: u32_exp
            type (uint_32) :: uwrk, uwrk1
            integer        :: jwrk
            intrinsic mod
!
            if (jpow < 0) then
               uwrk = uint_32 (0)
               if (u32_teq (u32a, uwrk)) then
                  u32_exp = uint_32 (0) ! Should raise exception, but...
               else
                  uwrk = uint_32 (1)
                  if (u32_teq (u32a, uwrk)) then
                     u32_exp = uint_32 (1)
                  else
                     u32_exp = uint_32 (0)
                  endif
               endif
            elseif (jpow == 0) then
               u32_exp = uint_32 (1)
            else
               uwrk = u32a
               uwrk1 = u32a
               jwrk = jpow - 1
               if (jwrk /= 0) then
                  do
                     if (mod (jwrk, 2) /= 0) then
                        uwrk = u32_mul (uwrk, uwrk1)
                     endif
                     uwrk1 = u32_mul (uwrk1, uwrk1)
                     jwrk = jwrk / 2
                     if (jwrk == 0) exit
                  enddo
               endif
               u32_exp = uwrk
            endif
         end function u32_exp
!
! _______bit functions_________________________________________________
! _____________________________________________________________________
!
! _______bit_size______________________________________________________
!
         integer function u32_bsz (i)
            type (uint_32), intent (in) :: i
!
            u32_bsz = lw
         end function u32_bsz
!
! _______btest_________________________________________________________
!
         logical function u32_tbt (i, pos)
            type (uint_32), intent (in) :: i
            integer, intent (in)        :: pos
            intrinsic btest
!
            u32_tbt = btest (i%val, pos)
         end function u32_tbt
!
! _______digits_______________________________________________________
!
         integer function u32_dgt (x)
            type (uint_32), intent (in) :: x
!
            u32_dgt = nd
         end function u32_dgt
!
! _______huge_________________________________________________________
!
         function u32_hge (x)
            type (uint_32), intent (in) :: x
            type (uint_32)              :: u32_hge
            intrinsic ibits
!
            u32_hge%val = ibits (-1, 0, lw)
         end function u32_hge
!
! _______iand__________________________________________________________
!
         function u32_and (i, j)
            type (uint_32), intent (in) :: i, j
            type (uint_32)              :: u32_and
            intrinsic iand
!
            u32_and%val = iand (i%val, j%val)
         end function u32_and
!
! _______ibclr_________________________________________________________
!
         function u32_bcl (i, pos)
            type (uint_32), intent (in) :: i
            integer, intent (in)        :: pos
            type (uint_32)              :: u32_bcl
            intrinsic ibclr
!
            u32_bcl%val = ibclr (i%val, pos)
         end function u32_bcl
!
! _______ibits_________________________________________________________
!
         function u32_bts (i, pos, len)
            type (uint_32), intent (in) :: i
            integer, intent (in)        :: pos, len
            type (uint_32)              :: u32_bts
            intrinsic ibits
!
            u32_bts%val = ibits (i%val, pos, len)
         end function u32_bts
!
! _______ibset_________________________________________________________
!
         function u32_bst (i, pos)
            type (uint_32), intent (in) :: i
            integer, intent (in)        :: pos
            type (uint_32)              :: u32_bst
            intrinsic ibset
!
            u32_bst%val = ibset (i%val, pos)
         end function u32_bst
!
! _______ieor__________________________________________________________
!
         function u32_xor (i, j)
            type (uint_32), intent (in) :: i, j
            type (uint_32)              :: u32_xor
            intrinsic ieor
!
            u32_xor%val = ieor (i%val, j%val)
         end function u32_xor
!
! _______ior___________________________________________________________
!
         function u32_ior (i, j)
            type (uint_32), intent (in) :: i, j
            type (uint_32)              :: u32_ior
            intrinsic ior
!
            u32_ior%val = ior (i%val, j%val)
         end function u32_ior
!
! _______ishft_________________________________________________________
!
         function u32_sft (i, shift)
            type (uint_32), intent (in) :: i
            integer, intent (in)        :: shift
            type (uint_32)              :: u32_sft
            intrinsic ishft
!
            u32_sft%val = ishft (i%val, shift)
         end function u32_sft
!
! _______ishftc________________________________________________________
!
         function u32_sfc (i, shift)
            type (uint_32), intent (in) :: i
            integer, intent (in)        :: shift
            type (uint_32)              :: u32_sfc
            intrinsic ishftc
!
            u32_sfc%val = ishftc (i%val, shift, size = lw)
         end function u32_sfc
!
! _______mvbits________________________________________________________
!
         subroutine u32_mvb (from, frompos, len, to, topos)
            type (uint_32), intent (in)  :: from
            integer, intent (in)         :: frompos, len, topos
            type (uint_32), intent (out) :: to
            intrinsic mvbits
!
            call mvbits (from%val, frompos, len, to%val, topos)
         end subroutine u32_mvb
!
! _______not___________________________________________________________
!
         function u32_not (i)
            type (uint_32), intent (in) :: i
            type (uint_32)              :: u32_not
            intrinsic not
!
            u32_not%val = not (i%val)
         end function u32_not
!
! _______radix_________________________________________________________
!
         integer function u32_rdx (x)
            type (uint_32), intent (in) :: x
!
            u32_rdx = 2
         end function u32_rdx
!
! _______radix_________________________________________________________
!
         integer function u32_rng (x)
            type (uint_32), intent (in) :: x
!
            u32_rng = nd - 1
         end function u32_rng
!
! _______tiny__________________________________________________________
!
         function u32_tny (x)
            type (uint_32), intent (in) :: x
            type (uint_32)              :: u32_tny
!
            u32_tny = uint_32 (1)
         end function u32_tny
!
! _______random_numbers________________________________________________
! _____________________________________________________________________
!
         subroutine u32_rnd (harvest)
            type (uint_32), intent (out) :: harvest
            type (uint_32)               :: utmp
!
! ... build the value from twice 16 highest bits of random 1-2^31
!
            usee = u32_add (u32_mul (usee, umul), uadd)
            usee = u32_bts (usee, 0, lwm1)
            utmp = u32_bts (usee, lh-1, lh)
            usee = u32_add (u32_mul (usee, umul), uadd)
            usee = u32_bts (usee, 0, lwm1)
            harvest = u32_ior (u32_sft (utmp, lh),                      &
     &                         u32_bts (usee, lh-1, lh))
         end subroutine u32_rnd
!
         subroutine u32_see (size, put, get)
            integer, optional, intent (out)        :: size
            type (uint_32), optional, intent (in)  :: put (:)
            type (uint_32), optional, intent (out) :: get (:)
            integer   :: jsee
            intrinsic system_clock
!
            if (present (size)) then
               size = 1
            endif
            if (present (get)) then
               get (1) = usee
            endif
            if (present (put)) then
               usee = put (1)
            else
               call system_clock (count = jsee)
               call u32_itu (usee, jsee)
            endif
         end subroutine u32_see
!
! _______max, min, mod, dim____________________________________________
! _____________________________________________________________________
!
         function u32_max (a1, a2)
            type (uint_32), intent (in) :: a1, a2
            type (uint_32)              :: u32_max
!
            if (u32_tge (a1, a2)) then
               u32_max = a1
            else
               u32_max = a2
            endif
         end function u32_max
!
         function u32_min (a1, a2)
            type (uint_32), intent (in) :: a1, a2
            type (uint_32)              :: u32_min
!
            if (u32_tle (a1, a2)) then
               u32_min = a1
            else
               u32_min = a2
            endif
         end function u32_min
!
         function u32_mod (a, p)
            type (uint_32), intent (in) :: a, p
            type (uint_32)              :: u32_mod
!
            u32_mod = u32_sub (a, u32_mul (u32_div (a, p), p))
         end function u32_mod
!
         function u32_dim (x, y)
            type (uint_32), intent (in) :: x, y
            type (uint_32)              :: u32_dim
!
            if (u32_tge (y, x)) then
               u32_dim = uint_32 (0)
            else
               u32_dim = u32_sub (x, y)
            endif
         end function u32_dim
!
! Section 2. Array operations
! _____________________________________________________________________
!
!
! _______nothing to do______________________________________________
! _____________________________________________________________________
!
         function u32a_nop (a)
            type (uint_32), intent (in), dimension (:) :: a
            type (uint_32), dimension (1:size (a))     :: u32a_nop
!
            u32a_nop = a
         end function u32a_nop
!
! _______integer to unsigned assignment_____________________________
! _____________________________________________________________________
!
         subroutine u32a_itu (u, i)
            integer, intent (in), dimension (:)                   :: i
            type (uint_32), intent (out), dimension (1:size (i))  :: u
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end subroutine u32a_itu
!
         subroutine u32a_i1tu (u, i)
            integer(selected_int_kind(2)), intent(in), dimension(:) :: i
            type (uint_32), intent (out), dimension (1:size (i))  :: u
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end subroutine u32a_i1tu
!
         subroutine u32a_i2tu (u, i)
            integer(selected_int_kind(4)), intent(in), dimension(:) :: i
            type (uint_32), intent (out), dimension (1:size (i))  :: u
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end subroutine u32a_i2tu
!
         subroutine u32a_i4tu (u, i)
            integer(selected_int_kind(9)), intent(in), dimension(:) :: i
            type (uint_32), intent (out), dimension (1:size (i))  :: u
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end subroutine u32a_i4tu
!
! _______explicit integer to unsigned conversions______________________
!
         function u32a_unsi1 (i)
            integer(selected_int_kind(2)), intent(in), dimension(:) :: i
            type (uint_32), dimension (1:size (i))       :: u32a_unsi1
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u32a_unsi1(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end function u32a_unsi1
!
         function u32a_unsi2 (i)
            integer(selected_int_kind(4)), intent(in), dimension(:) :: i
            type (uint_32), dimension (1:size (i))       :: u32a_unsi2
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u32a_unsi2(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end function u32a_unsi2
!
         function u32a_unsi4 (i)
            integer(selected_int_kind(9)), intent(in), dimension(:) :: i
            type (uint_32), dimension (1:size (i))       :: u32a_unsi4
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u32a_unsi4(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end function u32a_unsi4
!
         function u32a_unsi (i)
            integer, intent (in), dimension(:)           :: i
            type (uint_32), dimension (1:size (i))       :: u32a_unsi
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (i)
               u32a_unsi(ival)%val = ibits (i (ival), 0, nbts)
            enddo
         end function u32a_unsi
!
! _______unsigned to integer assignment________________________________
!
         subroutine u32a_uti (i, u)
            type (uint_32), intent (in), dimension (:)    :: u
            integer, intent (out), dimension (1:size (u)) :: i
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (u)
               i (ival) = ibits (u(ival)%val, 0, nbts)
            enddo
         end subroutine u32a_uti
!
         subroutine u32a_uti1 (i, u)
            type (uint_32), intent (in), dimension (:)    :: u
            integer (selected_int_kind (2)), intent (out),              &
     &                             dimension (1:size (u)) :: i
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (u)
               i (ival) = ibits (u(ival)%val, 0, nbts)
            enddo
         end subroutine u32a_uti1
!
         subroutine u32a_uti2 (i, u)
            type (uint_32), intent (in), dimension (:)    :: u
            integer (selected_int_kind (4)), intent (out),              &
     &                             dimension (1:size (u)) :: i
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (u)
               i (ival) = ibits (u(ival)%val, 0, nbts)
            enddo
         end subroutine u32a_uti2
!
         subroutine u32a_uti4 (i, u)
            type (uint_32), intent (in), dimension (:)    :: u
            integer (selected_int_kind (9)), intent (out),              &
     &                             dimension (1:size (u)) :: i
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (u)
               i (ival) = ibits (u(ival)%val, 0, nbts)
            enddo
         end subroutine u32a_uti4
!
! _______explicit unsigned to integer conversion________________________
!
         function u32a_int (u)
            type (uint_32), intent (in), dimension (:)    :: u
            integer, dimension (1:size (u))               :: u32a_int
            integer      :: ival, nbts
            intrinsic min, int, bit_size, ibits
!
            nbts = min (lw, int (bit_size (u32a_int), kind = kind (lw)))
            do ival = 1, size (u)
               u32a_int (ival) = ibits (u(ival)%val, 0, nbts)
            enddo
!
         end function u32a_int
!
! _______real to unsigned assignment________________________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32a_rtu (u, r)
            real, intent (in), dimension (:)                     :: r
            type (uint_32), intent (out), dimension (1:size (r)) :: u
            integer            :: jwrk, ival, nbts
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, min, floor, ibset, bit_size, ibits
!
            nbts = min (lw, bit_size (jwrk))
            do ival = 1, size (r)
               if (abs (r (ival)) < one) then
                  u%val = 0
               else
                  xwrk = r (ival)
                  xwrk = modulo (xwrk, x232)
                  if (xwrk == x231) then
                     jwrk = ibset (0, lwm1)
                  elseif (xwrk > x231) then
                     xwrk = xwrk - x232
                     jwrk = floor (xwrk)
                  else
                     jwrk = floor (xwrk)
                  endif
                  u(ival)%val = ibits (jwrk, 0, nbts)
               endif
            enddo
         end subroutine u32a_rtu
!
! _______explicit real to unsigned conversion__________________________
!
         function u32a_unsr (r)
            real, intent (in), dimension (:)          :: r
            type (uint_32), dimension (1:size (r))    :: u32a_unsr
            type (uint_32), dimension (1:size (r))    :: uwrk
!
            call u32a_rtu (uwrk (:), r (:))
            u32a_unsr (:) = uwrk (:)
         end function u32a_unsr
!
! _______unsigned to real assignment___________________________________
!
         subroutine u32a_utr (r, u)
            type (uint_32), intent (in), dimension (:)  :: u
            real, intent (out), dimension (1:size (u))  :: r
            real (selected_real_kind (nd, nd))   :: xwrk
            integer            :: ival
            intrinsic real, ibits
!
            do ival = 1, size (u)
               xwrk = real (ibits (u(ival)%val, 0, lwm1),               &
     &                      kind = kind (xwrk))
               if (u32_tbt (u (ival), lwm1)) then
                  xwrk = x231 + xwrk
               endif
               r (ival) = xwrk
            enddo
         end subroutine u32a_utr
!
! _______explicit unsigned to real conversion_______________________
!
         function u32a_rea (a)
            type (uint_32), intent (in), dimension (:)  :: a
            real, dimension (1:size (a))                :: u32a_rea
            real, dimension (1:size (a))                :: rwrk
!
            call u32a_utr (rwrk (:), a (:))
            u32a_rea (:) = rwrk (:)
         end function u32a_rea
!
! _______double to unsigned assignment______________________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32a_dtu (u, r)
            real (kind (1.0d0)), intent (in), dimension (:)     :: r
            type (uint_32), intent (out), dimension (1:size(r)) :: u
            integer            :: jwrk, ival, nbts
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, bit_size, min, modulo, ibset, floor, ibits
!
            nbts = min (lw, bit_size (jwrk))
            do ival = 1, size (r)
               if (abs (r (ival)) < one) then
                  u (ival) = uint_32 (0)
               else
                  xwrk = r (ival)
                  xwrk = modulo (xwrk, x232)
                  if (xwrk == x231) then
                     jwrk = ibset (0, lwm1)
                  elseif (xwrk > x231) then
                     xwrk = xwrk - x232
                     jwrk = floor (xwrk)
                  else
                     jwrk = floor (xwrk)
                  endif
                  u (ival)%val = ibits (jwrk, 0, nbts)
               endif
            enddo
         end subroutine u32a_dtu
!
! _______explicit double to unsigned conversion________________________
!
         function u32a_unsd (r)
            real (kind = kind (1.0d0)), intent (in), dimension (:) :: r
            type (uint_32), dimension (1:size(r))   :: u32a_unsd
            type (uint_32), dimension (1:size(r))   :: uwrk
!
            call u32a_dtu (uwrk (:), r (:))
            u32a_unsd (:) = uwrk (:)
         end function u32a_unsd
!
! _______unsigned to double assignment______________________________
!
         subroutine u32a_utd (r, u)
            type (uint_32), intent (in), dimension (:)              :: u
            real (kind (1.0d0)), intent (out), dimension (1:size(u)):: r
            real (selected_real_kind (nd, nd))   :: xwrk
            integer            :: ival
            intrinsic real, ibits
!
            do ival = 1, size (u)
               xwrk = real (ibits (u (ival)%val, 0, lwm1),              &
     &                      kind = kind (xwrk))
               if (u32_tbt (u (ival), lwm1)) then
                  xwrk = x231 + xwrk
               endif
               r (ival) = xwrk
            enddo
         end subroutine u32a_utd
!
! _______explicit unsigned to double conversion_____________________
!
         function u32a_dbl (a)
            type (uint_32), intent (in), dimension (:)  :: a
            real (kind (1.0d0)), dimension (1:size (a)) :: u32a_dbl
            real (kind (1.0d0)), dimension (1:size (a)) :: rwrk
!
            call u32a_utd (rwrk (:), a (:))
            u32a_dbl (:) = rwrk (:)
         end function u32a_dbl
!
! _______complex to unsigned assignment_____________________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32a_ctu (u, c)
            complex, intent (in), dimension (:)                  :: c
            type (uint_32), intent (out), dimension (1:size (c)) :: u
            real                                 :: r
            integer                              :: jwrk, ival, nbts
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, real, modulo, ibset, ibits, floor, min,      &
     &                bit_size
!
            nbts = min (lw, int (bit_size (i), kind = kind (lw)))
            do ival = 1, size (c)
               r = real (c (ival), kind = kind (r))
               if (abs (r) < one) then
                  u (ival) = uint_32 (0)
               else
                  xwrk = r
                  xwrk = modulo (xwrk, x232)
                  if (xwrk == x231) then
                     jwrk = ibset (0, lwm1)
                  elseif (xwrk > x231) then
                     xwrk = xwrk - x232
                     jwrk = floor (xwrk)
                  else
                     jwrk = floor (xwrk)
                  endif
                  u (ival)%val = ibits (jwrk, 0, nbts)
               endif
            enddo
         end subroutine u32a_ctu
!
! _______explicit complex to unsigned conversion_______________________
!
         function u32a_unsc (c)
            complex, intent (in), dimension (:)     :: c
            type (uint_32), dimension (1:size (c))  :: u32a_unsc
            type (uint_32), dimension (1:size (c))  :: uwrk
!
            call u32a_ctu (uwrk (:), c (:))
            u32a_unsc (:) = uwrk (:)
         end function u32a_unsc
!
! _______unsigned to complex assignment______________________________
!
         subroutine u32a_utc (c, u)
            type (uint_32), intent (in), dimension (:)    :: u
            complex, intent (out), dimension (1:size (u)) :: c
            integer                                       :: ival
            real (selected_real_kind (nd, nd))            :: r
            intrinsic cmplx, ibits, real
!
            do ival = 1, size (u)
               r = real (ibits (u (ival)%val, 0, lwm1), kind = kind (r))
               if (u32_tbt (u (ival), lwm1)) then
                  xwrk = x231 + xwrk
               endif
               c (ival) = cmplx (xwrk, 0.0, kind = kind (c))
            enddo
         end subroutine u32a_utc
!
! _______double complex to unsigned assignment______________________
! ******* Assumes radix is power of 2 ******************************
!
         subroutine u32a_ktu (u, c)
            complex (kind (1.0d0)), intent (in), dimension (:)   :: c
            type (uint_32), intent (out), dimension (1:size (c)) :: u
            real (kind (1.0d0))                                  :: r
            integer                              :: jwrk, ival, nbts
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic abs, floor, ibits, ibset, real, modulo, min,      &
     &                bit_size
!
            nbts = min (lw, int (bit_size (jwrk), kind = kind (lw)))
            do ival = 1, size (u)
               r = real (c (ival), kind = kind (r))
               if (abs (r) < one) then
                  u (ival) = uint_32 (0)
               else
                  xwrk = r
                  xwrk = modulo (xwrk, x232)
                  if (xwrk == x231) then
                     jwrk = ibset (0, lwm1)
                  elseif (xwrk > x231) then
                     xwrk = xwrk - x232
                     jwrk = floor (xwrk)
                  else
                     jwrk = floor (xwrk)
                  endif
                  u (ival)%val = ibits (jwrk, 0, nbts)
               endif
            enddo
         end subroutine u32a_ktu
!
! _______explicit double complex to unsigned conversion________________
!
         function u32a_unsk (c)
            complex (kind (1.0d0)), intent (in), dimension (:) :: c
            type (uint_32), dimension (1:size (c))     :: u32a_unsk
            type (uint_32), dimension (1:size (c))     :: uwrk
!
            call u32a_ktu (uwrk (:), c (:))
            u32a_unsk (:) = uwrk (:)
         end function u32a_unsk
!
! _______unsigned to double complex assignment_________________________
!
         subroutine u32a_utk (c, u)
            type (uint_32), intent (in), dimension (:)         :: u
            complex (kind = kind (1.0d0)), intent (out),                &
     &                                  dimension (1:size (u)) :: c
            integer            :: ival
            real (selected_real_kind (nd, nd))   :: xwrk
            intrinsic cmplx, ibits, real
!
            do ival = 1, size (u)
               xwrk = real(ibits(u(ival)%val, 0, lwm1), kind=kind(xwrk))
               if (u32_tbt (u (ival), lwm1)) then
                  xwrk = x231 + xwrk
               endif
               c (ival) = cmplx (xwrk, 0.0, kind = kind (c))
            enddo
         end subroutine u32a_utk
!
! _______explicit unsigned to complex conversions______________________
!
         function u32a_cpl (x, y)
            type (uint_32), intent (in), dimension (:)            :: x
            type (uint_32), intent (in), optional, dimension (:)  :: y
            complex, dimension (1:size (x))        :: u32a_cpl
            integer                                :: ival
            real (kind = lrk)                      :: rx, ry
            intrinsic cmplx
!
            if (present (y)) then
               do ival = 1, size (x)
                  rx = u32_dbl (x (ival))
                  ry = u32_dbl (y (ival))
                  u32a_cpl (ival) = cmplx (rx, ry)
               enddo
            else
               do ival = 1, size (x)
                  rx = u32_dbl (x (ival))
                  u32a_cpl (ival) = cmplx (rx, 0.0)
               enddo
            endif
         end function u32a_cpl
!
! _______convert string to unsigned_____________________________________
!
         subroutine u32a_stu (u, s)
            character (len=*), intent (in), dimension (:)       :: s
            type (uint_32), intent (out), dimension (1:size(s)) :: u
            type (uint_32)     :: uwrk1, uwrk2
            integer            :: i, i0, i1, j, ival
            intrinsic len_trim, index
!
            do ival = 1, size (s)
               i0 = 1
               i1 = len_trim (s (ival))
               uwrk1 = uint_32 (0)
               do i = i0, i1
                  if (s(ival)(i:i) /= ' ' .and. s(ival)(i:i) /= '+')then
                     i0 = i
                     exit
                  endif
               enddo
               do i = i0, i1
                  j = index (zdgt, s (ival)(i:i)) - 1
                  if (j < 0) then
                     exit
                  else
                     uwrk1 = u32_mul (uwrk1, uten)
                     if (j > 0) then
                        uwrk2 = uint_32 (j)
                        uwrk1 = u32_add (uwrk1, uwrk2)
                     endif
                  endif
               enddo
               u (ival) = uwrk1
            enddo  
         end subroutine u32a_stu
!
! _______explicit string to unsigned conversion________________________
!
         function u32a_unss (s)
            character (len = *), intent (in), dimension (:) :: s
            type (uint_32), dimension (1:size(s))           :: u32a_unss
            type (uint_32), dimension (1:size(s))           :: uwrk
!
            call u32a_stu (uwrk (:), s (:))
            u32a_unss (:) = uwrk (:)
         end function u32a_unss
!
! _______convert unsigned to string____________________________________
!
         subroutine u32a_uts (s, u)
            type (uint_32), intent (in), dimension (:)             :: u
            character (len=*), intent (out), dimension (1:size(u)) :: s
            type (uint_32)     :: uwrk, uwrk1, uwrk2
            integer            :: i, i0, i1, lens, ival
            intrinsic repeat, len
!
            lens = len (s)
            do ival = 1, size (u)
               i0 = 1
               i1 = nd
               uwrk2 = uten
               do i = i0, i1
                  if (u32_tlt (u (ival), uwrk2)) then
                     i1 = i
                     exit
                  endif
                  uwrk2 = u32_mul (uten, uwrk2)
               enddo
               if (i1 > lens) then
                  s (ival) = repeat ('*', lens)
               else
                  s = ' '
                  uwrk1 = u (ival)
                  uwrk2 = u32_exp (uten, i1 - 1)
                  i0 = lens - i1 + 1
                  i1 = lens - 1
                  do i = i0, i1
                     uwrk = u32_div (uwrk1, uwrk2)
                     s (ival)(i:i) = zdgt (uwrk%val+1:uwrk%val+1)
                     uwrk = u32_mul (uwrk, uwrk2)
                     uwrk1 = u32_sub (uwrk1, uwrk)
                     uwrk2 = u32_div (uwrk2, uten)
                  enddo
                  s (ival)(lens:lens) = zdgt (uwrk1%val+1:uwrk1%val+1)
               endif
            enddo
         end subroutine u32a_uts
!
! _______convert unsigned to char______________________________________
!
         function u32a_chr (i)
            type (uint_32), intent (in), dimension (:) :: i
            character (len=1), dimension (1:size (i))  :: u32a_chr
            intrinsic  char
            integer              :: ival
!
            do ival = 1, size (i)
               u32a_chr (ival) = char (i (ival)%val)
            enddo
         end function u32a_chr
!
! _______convert char to unsigned______________________________________
!
         function u32a_uch (c)
            character (len=1), intent (in), dimension (:) :: c
            type (uint_32), dimension (1:size (c))        :: u32a_uch
            intrinsic ichar
            integer              :: ival
!
            do ival = 1, size (c)
               u32a_uch (ival)%val = ichar (c (ival))
            enddo
         end function u32a_uch
!
! _______convert unsigned to ascii_____________________________________
!
         function u32a_ach (i)
            type (uint_32), intent (in), dimension (:) :: i
            character (len=1), dimension (1:size (i))  :: u32a_ach
            intrinsic achar
            integer              :: ival
!
            do ival = 1, size (i)
               u32a_ach (ival) = achar (i (ival)%val)
            enddo
         end function u32a_ach
!
! _______convert ascii to unsigned_____________________________________
!
         function u32a_uac (c)
            character (len=1), intent (in), dimension (:) :: c
            type (uint_32), dimension (1:size (c))        :: u32a_uac
            intrinsic iachar
            integer              :: ival
!
            do ival = 1, size (c)
               u32a_uac (ival)%val = iachar (c (ival))
            enddo
         end function u32a_uac
!
! _______test for equality______________________________________________
!
         function u32a_teq (u32a, u32b)
            type (uint_32), intent (in), dimension (:)  :: u32a, u32b
            logical, dimension (1:size(u32a))           :: u32a_teq
            integer                                     :: i
            intrinsic ibits
!
            do i = 1, size (u32a)
                u32a_teq (i) =                                          &
     &      (ibits (u32a(i)%val,  0, lh) == ibits (u32b(i)%val,  0, lh))&
     & .and.(ibits (u32a(i)%val, lh, lh) == ibits (u32b(i)%val, lh, lh))
             enddo
         end function u32a_teq
!
! _______test for inequality____________________________________________
!
         function u32a_tne (u32a, u32b)
            type (uint_32), intent (in), dimension (:)  :: u32a, u32b
            logical, dimension (1:size(u32a))           :: u32a_tne
            integer                                     :: i
            intrinsic ibits
!
            do i = 1, size (u32a)
!
                u32a_tne (i) =                                          &
     &      (ibits (u32a(i)%val,  0, lh) /= ibits (u32b(i)%val,  0, lh))&
     & .or. (ibits (u32a(i)%val, lh, lh) /= ibits (u32b(i)%val, lh, lh))
            enddo
         end function u32a_tne
!
! _______test for greater than__________________________________________
!
         function u32a_tgt (u32a, u32b)
            type (uint_32), intent (in), dimension (:)  :: u32a, u32b
            logical, dimension (1:size(u32a))           :: u32a_tgt
            integer            :: jwrk1, jwrk2, ival
            intrinsic ibits
!
            do ival = 1, size (u32a)
               jwrk1 = ibits (u32a (ival)%val, lh, lh)
               jwrk2 = ibits (u32b (ival)%val, lh, lh)
               if (jwrk1 == jwrk2) then
                  u32a_tgt (ival) = ibits (u32a (ival)%val,  0, lh) >   &
     &                                ibits (u32b (ival)%val,  0, lh) 
               else
                  u32a_tgt (ival) = jwrk1 > jwrk2
               endif
            enddo
         end function u32a_tgt
!
! _______test for lower than____________________________________________
!
         function u32a_tlt (u32a, u32b)
            type (uint_32), intent (in), dimension (:)  :: u32a, u32b
            logical, dimension (1:size(u32a))           :: u32a_tlt
            integer            :: jwrk1, jwrk2, ival
            intrinsic ibits
!
            do ival = 1, size (u32a)
               jwrk1 = ibits (u32a (ival)%val, lh, lh)
               jwrk2 = ibits (u32b (ival)%val, lh, lh)
               if (jwrk1 == jwrk2) then
                  u32a_tlt (ival) = ibits (u32a (ival)%val,  0, lh) <   &
     &                                ibits (u32b (ival)%val,  0, lh) 
               else
                  u32a_tlt (ival) = jwrk1 < jwrk2
               endif
            enddo
         end function u32a_tlt
!
! _______test for greater or equal_____________________________________
!
         function u32a_tge (u32a, u32b)
            type (uint_32), intent (in), dimension (:)  :: u32a, u32b
            logical, dimension (1:size(u32a))           :: u32a_tge
            integer            :: jwrk1, jwrk2, ival
            intrinsic ibits
!
            do ival = 1, size (u32a)
               jwrk1 = ibits (u32a (ival)%val, lh, lh)
               jwrk2 = ibits (u32b (ival)%val, lh, lh)
               if (jwrk1 == jwrk2) then
                  u32a_tge (ival) = ibits (u32a (ival)%val,  0, lh) >=  &
     &                                 ibits (u32b (ival)%val,  0, lh) 
               else
                  u32a_tge (ival) = jwrk1 > jwrk2
               endif
            enddo
         end function u32a_tge
!
! _______test for lower or equal_______________________________________
!
         function u32a_tle (u32a, u32b)
            type (uint_32), intent (in), dimension (:)  :: u32a, u32b
            logical, dimension (1:size(u32a))           :: u32a_tle
            integer            :: jwrk1, jwrk2, ival
            intrinsic ibits
!
            do ival = 1, size (u32a)
               jwrk1 = ibits (u32a (ival)%val, lh, lh)
               jwrk2 = ibits (u32b (ival)%val, lh, lh)
               if (jwrk1 == jwrk2) then
                  u32a_tle (ival) = ibits (u32a (ival)%val,  0, lh) <=  &
     &                                 ibits (u32b (ival)%val,  0, lh) 
               else
                  u32a_tle (ival) = jwrk1 < jwrk2
               endif
            enddo
         end function u32a_tle
!
! _______addition of two unsigned numbers arrays________________________
!
         function u32a_add (u32a, u32b)
            type (uint_32), intent (in), dimension (:)   :: u32a, u32b
            type (uint_32), dimension (1:size(u32a))     :: u32a_add
            integer            :: jwrk1, jwrk2, ival
            intrinsic ibits, ishft
!
            do ival = 1, size (u32a)
               jwrk1 = ibits (u32a (ival)%val,  0, lh) +                &
     &                 ibits (u32b (ival)%val,  0, lh)
               jwrk2 = ibits (u32a (ival)%val, lh, lh) +                &
     &                 ibits (u32b (ival)%val, lh, lh) +                &
     &                 ibits (jwrk1, lh, lh)
               jwrk1 = ibits (jwrk1,  0, lh)
               jwrk2 = ishft (jwrk2, lh)
               u32a_add (ival)%val = ior (jwrk1, jwrk2)
            enddo
         end function u32a_add
!
! _______substraction of two unsigned numbers arrays____________________
!
         function u32a_sub (u32a, u32b)
            type (uint_32), intent (in), dimension (:)   :: u32a, u32b
            type (uint_32), dimension (1:size(u32a))     :: u32a_sub
            integer            :: jwrk1, jwrk2, jwrkb, ival
            intrinsic ibits, ishft, not, ior
!
            do ival = 1, size (u32a)
               jwrkb = not (ibits (u32b (ival)%val,  0, lw))
               jwrk1 = ibits (u32a (ival)%val,  0, lh) +                &
     &                 ibits (jwrkb,  0, lh) + 1
               jwrk2 = ibits (u32a (ival)%val, lh, lh) +                &
     &                 ibits (jwrkb, lh, lh) + ibits (jwrk1, lh, lh)
               jwrk1 = ibits (jwrk1,  0, lh)
               jwrk2 = ishft (jwrk2, lh)
               u32a_sub (ival)%val = ior (jwrk1, jwrk2)
            enddo
         end function u32a_sub
!
! _______multiplication of two unsigned numbers arrays__________________
!
         function u32a_mul (u32a, u32b)
            type (uint_32), intent (in), dimension (:)   :: u32a, u32b
            type (uint_32), dimension (1:size(u32a))     :: u32a_mul
            integer            :: jwrka1, jwrka2, jwrkb1, jwrkb2
            integer            :: jwrk1, jwrk2, ival
            intrinsic ibits, ishft, ior
!
            do ival = 1, size (u32a)
               jwrka1 = ibits (u32a (ival)%val,  0, lh)
               jwrkb1 = ibits (u32b (ival)%val,  0, lh)
               jwrka2 = ibits (u32a (ival)%val, lh, lh)
               jwrkb2 = ibits (u32b (ival)%val, lh, lh)
               jwrk1 = jwrka1 * jwrkb1
               jwrk2 = jwrka1 * jwrkb2 + jwrka2 * jwrkb1 +              &
     &                 ibits (jwrk1, lh, lh)
               jwrk1 = ibits (jwrk1,  0, lh)
               jwrk2 = ishft (jwrk2, lh)
               u32a_mul (ival)%val = ior (jwrk1, jwrk2)
            enddo
         end function u32a_mul
!
! _______division of two unsigned numbers arrays______________________
!
         function u32a_div (u32a, u32b)
            type (uint_32), intent (in), dimension (:)   :: u32a, u32b
            type (uint_32), dimension (1:size(u32a))     :: u32a_div
            type (uint_32) :: uwrk1, uwrk2, uwrk
            integer        :: ival
            intrinsic btest, ibits

            do ival = 1, size (u32a)
               uwrk = uint_32 (0)
               if (u32_teq (u32b (ival), uwrk)) then
                  u32a_div (ival) = uint_32 (0) 
                                 ! Should raise exception, but ...
               elseif (u32_tlt (u32a (ival), u32b (ival))) then
                  u32a_div (ival) = uint_32 (0)
               else
                  if (btest (u32a (ival)%val, lwm1)) then
!
! ... a large, b large
                     if (btest (u32b (ival)%val, lwm1)) then
                        u32a_div (ival) = uint_32 (1)
!
! ... a large, b small, first divide a/2
                     else
                        uwrk1%val = ibits (u32a (ival)%val, 1, lwm1) /  &
     &                              ibits (u32b (ival)%val, 0, lwm1)
                        uwrk2 = uint_32 (2)
                        uwrk2 = u32_mul (uwrk1, uwrk2)
                        uwrk1 = u32_mul (u32b (ival), uwrk2)
                        uwrk  = u32_sub (u32a (ival),  uwrk1)
                        if (u32_tlt (uwrk, u32b (ival))) then
                           u32a_div (ival) = uwrk2
                        else
                           uwrk = u32_sub (uwrk, u32b (ival))
                           uwrk1%val = 1 + ibits (uwrk%val, 0, lwm1) /  &
     &                                  ibits (u32b (ival)%val, 0, lwm1)
                           u32a_div (ival) = u32_add (uwrk1, uwrk2)
                        endif
                     endif
!
! ... a small, b small
                  else
                    u32a_div(ival)%val = ibits(u32a(ival)%val, 0, lwm1)/&
     &                                   ibits(u32b(ival)%val, 0, lwm1)
                  endif
               endif
            enddo
         end function u32a_div
!
! _______exponentiation, to integer power, of array_____________________
!
         function u32a_exp (u32a, jpow)
            type (uint_32), intent (in), dimension (:)   :: u32a
            integer, intent (in)                         :: jpow
            type (uint_32), dimension (1:size(u32a))     :: u32a_exp
            type (uint_32) :: uwrk, uwrk1
            integer        :: jwrk, ival
            intrinsic mod
!
            if (jpow < 0) then
               do ival = 1, size (u32a)
                  if (u32_teq (u32a (ival), uint_32 (0))) then
                     u32a_exp (ival) = uint_32 (0) 
                                     ! Should raise exception, but ...
                  else
                     if (u32_teq (u32a (ival), uint_32 (1))) then
                        u32a_exp (ival) = uint_32 (1)
                     else
                        u32a_exp (ival) = uint_32 (0)
                     endif
                  endif
               enddo
            elseif (jpow == 0) then
               u32a_exp (:) = uint_32 (1)
            else
               do ival = 1, size (u32a)
                  uwrk = u32a (ival)
                  uwrk1 = uwrk
                  jwrk = jpow - 1
                  if (jwrk /= 0) then
                     do
                        if (mod (jwrk, 2) /= 0) then
                           uwrk = u32_mul (uwrk, uwrk1)
                        endif
                        uwrk1 = u32_mul (uwrk1, uwrk1)
                        jwrk = jwrk / 2
                        if (jwrk == 0) exit
                     enddo
                  endif
                  u32a_exp (ival) = uwrk
               enddo
            endif
         end function u32a_exp
!
! _______bit functions_________________________________________________
! _____________________________________________________________________
!
!
! _______btest_________________________________________________________
!
         function u32a_tbt (i, pos)
            type (uint_32), intent (in), dimension (:) :: i
            integer, intent (in)                       :: pos
            logical, dimension (1:size(i))             :: u32a_tbt
            intrinsic btest
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_tbt (ival) = btest (i (ival)%val, pos)
            enddo
         end function u32a_tbt
!
! _______iand__________________________________________________________
!
         function u32a_and (i, j)
            type (uint_32), intent (in), dimension (:) :: i, j
            type (uint_32), dimension (1:size(i))      :: u32a_and
            intrinsic iand
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_and (ival)%val = iand (i (ival)%val, j (ival)%val)
            enddo
         end function u32a_and
!
! _______ibclr_________________________________________________________
!
         function u32a_bcl (i, pos)
            type (uint_32), intent (in), dimension (:) :: i
            integer, intent (in)                       :: pos
            type (uint_32), dimension (1:size(i))      :: u32a_bcl
            intrinsic ibclr
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_bcl (ival)%val = ibclr (i (ival)%val, pos)
            enddo
         end function u32a_bcl
!
! _______ibits_________________________________________________________
!
         function u32a_bts (i, pos, len)
            type (uint_32), intent (in), dimension (:) :: i
            integer, intent (in)                       :: pos, len
            type (uint_32), dimension (1:size(i))      :: u32a_bts
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_bts (ival)%val = ibits (i (ival)%val, pos, len)
            enddo
         end function u32a_bts
!
! _______ibset_________________________________________________________
!
         function u32a_bst (i, pos)
            type (uint_32), intent (in), dimension (:) :: i
            integer, intent (in)                       :: pos
            type (uint_32), dimension (1:size(i))      :: u32a_bst
            intrinsic ibset
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_bst (ival)%val = ibset (i (ival)%val, pos)
            enddo
         end function u32a_bst
!
! _______ieor__________________________________________________________
!
         function u32a_xor (i, j)
            type (uint_32), intent (in), dimension (:) :: i, j
            type (uint_32), dimension (1:size(i))      :: u32a_xor
            intrinsic ieor
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_xor (ival)%val = ieor (i (ival)%val, j (ival)%val)
            enddo
         end function u32a_xor
!
! _______ior___________________________________________________________
!
         function u32a_ior (i, j)
            type (uint_32), intent (in), dimension (:) :: i, j
            type (uint_32), dimension (1:size(i))      :: u32a_ior
            intrinsic ior
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_ior (ival)%val = ior (i (ival)%val, j (ival)%val)
            enddo
         end function u32a_ior
!
! _______ishft_________________________________________________________
!
         function u32a_sft (i, shift)
            type (uint_32), intent (in), dimension (:) :: i
            integer, intent (in)                       :: shift
            type (uint_32), dimension (1:size(i))      :: u32a_sft
            intrinsic ishft
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_sft (ival)%val = ishft (i (ival)%val, shift)
            enddo
         end function u32a_sft
!
! _______ishftc________________________________________________________
!
         function u32a_sfc (i, shift)
            type (uint_32), intent (in), dimension (:) :: i
            integer, intent (in)                       :: shift
            type (uint_32), dimension (1:size(i))      :: u32a_sfc
            intrinsic ishftc
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_sfc(ival)%val = ishftc (i(ival)%val, shift, size=lw)
            enddo
         end function u32a_sfc
!
! _______mvbits________________________________________________________
!
         subroutine u32a_mvb (from, frompos, len, to, topos)
            type (uint_32), intent (in), dimension (:)  :: from
            integer, intent (in)                  :: frompos, len, topos
            type (uint_32), intent (out), dimension (:) :: to
            integer                      :: ival
            intrinsic mvbits
!
            do ival = 1, size (from)
               call mvbits (from (ival)%val, frompos, len,              &
     &                      to (ival)%val, topos)
            enddo
         end subroutine u32a_mvb
!
! _______not___________________________________________________________
!
         function u32a_not (i)
            type (uint_32), intent (in), dimension (:) :: i
            type (uint_32), dimension (1:size(i))      :: u32a_not
            intrinsic not
            integer                                    :: ival
!
            do ival = 1, size (i)
               u32a_not (ival)%val = not (i (ival)%val)
            enddo
         end function u32a_not
!
! _______random_numbers________________________________________________
! _____________________________________________________________________
!
         subroutine u32a_rnd (harvest)
            type (uint_32), intent (out), dimension (:) :: harvest
            type (uint_32)                              :: utmp
            integer                                     :: ival
!
! ... build the value from twice 16 highest bits of random 1-2^31
!
            do ival = 1, size (harvest)
               usee = u32_add (u32_mul (usee, umul), uadd)
               usee = u32_bts (usee, 0, lwm1)
               utmp = u32_bts (usee, lh-1, lh)
               usee = u32_add (u32_mul (usee, umul), uadd)
               usee = u32_bts (usee, 0, lwm1)
               harvest (ival) = u32_ior (u32_sft (utmp, lh),            &
     &                                   u32_bts (usee, lh-1, lh))
            enddo
         end subroutine u32a_rnd
!
! _______max, min, mod, dim____________________________________________
! _____________________________________________________________________
!
         function u32a_max (a1, a2)
            type (uint_32), intent (in), dimension (:) :: a1, a2
            type (uint_32), dimension (1:size (a1))    :: u32a_max
!
            where (u32a_tge (a1, a2))
               u32a_max = a1
            elsewhere
               u32a_max = a2
            endwhere
         end function u32a_max
!
         function u32a_min (a1, a2)
            type (uint_32), intent (in), dimension (:) :: a1, a2
            type (uint_32), dimension (1:size (a1))    :: u32a_min
!
            where (u32a_tle (a1, a2))
               u32a_min = a1
            elsewhere
               u32a_min = a2
            endwhere
         end function u32a_min
!
         function u32a_mod (a, p)
            type (uint_32), intent (in), dimension (:) :: a, p
            type (uint_32), dimension (1:size (a))     :: u32a_mod
            integer                                    :: ival
!
            do ival = 1, size (a)
               u32a_mod (ival) = u32_sub (a (ival), u32_mul (           &
     &                       u32_div (a (ival), p (ival)), p (ival)))
            enddo
         end function u32a_mod
!
         function u32a_dim (x, y)
            type (uint_32), intent (in), dimension (:) :: x, y
            type (uint_32), dimension (1:size (x))     :: u32a_dim
!
            where (u32a_tge (y, x))
               u32a_dim = uint_32 (0)
            elsewhere
               u32a_dim = u32a_sub (x, y)
            endwhere
         end function u32a_dim
!
! Section 3. Array functions
! _____________________________________________________________________
!
!
! _______maxloc, minloc________________________________________________
! _____________________________________________________________________
!
         function u32_mal (array, mask)
            type (uint_32), intent (in), dimension (:)    :: array
            logical, intent (in), optional, dimension (:) :: mask
            integer, dimension (1:1)                      :: u32_mal
            type (uint_32)                                :: uwrk
            integer                                       :: iarr, iwrk
!
            uwrk = uint_32 (0)
            iwrk = 0
            if (present (mask)) then
               do iarr = size (array), 1, -1
                  if (mask (iarr)) then
                     if (u32_tge (array (iarr), uwrk)) then
                        uwrk = array (iarr)
                        iwrk = iarr
                     endif
                  endif
               enddo
            else
               do iarr = size (array), 1, -1
                  if (u32_tge (array (iarr), uwrk)) then
                     uwrk = array (iarr)
                     iwrk = iarr
                  endif
               enddo
            endif
            u32_mal (1) = iwrk
         end function u32_mal
!
         function u32_mil (array, mask)
            type (uint_32), intent (in), dimension (:)    :: array
            logical, intent (in), optional, dimension (:) :: mask
            integer, dimension (1:1)                      :: u32_mil
            type (uint_32)                                :: uwrk
            integer                                       :: iarr, iwrk
!
            uwrk = u32_hge (uwrk)
            iwrk = 0
            if (present (mask)) then
               do iarr = size (array), 1, -1
                  if (mask (iarr)) then
                     if (u32_tle (array (iarr), uwrk)) then
                        uwrk = array (iarr)
                        iwrk = iarr
                     endif
                  endif
               enddo
            else
               do iarr = size (array), 1, -1
                  if (u32_tle (array (iarr), uwrk)) then
                     uwrk = array (iarr)
                     iwrk = iarr
                  endif
               enddo
            endif
            u32_mil (1) = iwrk
         end function u32_mil
!
! _______maxval, minval________________________________________________
! _____________________________________________________________________
!
         function u32_mav (array, mask)
            type (uint_32), intent (in), dimension (:)    :: array
            logical, intent (in), optional, dimension (:) :: mask
            type (uint_32)                                :: u32_mav
            type (uint_32)                                :: uwrk
            integer                                       :: iarr
!
            uwrk = uint_32 (0)
            if (present (mask)) then
               do iarr = 1, size (array)
                  if (mask (iarr)) then
                     if (u32_tgt (array (iarr), uwrk)) then
                        uwrk = array (iarr)
                     endif
                  endif
               enddo
            else
               do iarr = 1, size (array)
                  if (u32_tgt (array (iarr), uwrk)) then
                     uwrk = array (iarr)
                  endif
               enddo
            endif
            u32_mav = uwrk
         end function u32_mav
!
         function u32_miv (array, mask)
            type (uint_32), intent (in), dimension (:)    :: array
            logical, intent (in), optional, dimension (:) :: mask
            type (uint_32)                                :: u32_miv
            type (uint_32)                                :: uwrk
            integer                                       :: iarr
!
            uwrk = u32_hge (uwrk)
            if (present (mask)) then
               do iarr = 1, size (array)
                  if (mask (iarr)) then
                     if (u32_tlt (array (iarr), uwrk)) then
                        uwrk = array (iarr)
                     endif
                  endif
               enddo
            else
               do iarr = 1, size (array)
                  if (u32_tlt (array (iarr), uwrk)) then
                     uwrk = array (iarr)
                  endif
               enddo
            endif
            u32_miv = uwrk
         end function u32_miv
!
! _______dot_product___________________________________________________
!
         function u32_dot (vector_a, vector_b)
            type (uint_32), intent (in), dimension (:) :: vector_a
            type (uint_32), intent (in), dimension (:) :: vector_b
            type (uint_32)                             :: u32_dot
            type (uint_32)                             :: uwrk
            integer                                    :: iarr, larr
!
            larr = size (vector_a)
            if (larr <= 0) then
               u32_dot = uint_32 (0)
            else
               uwrk = u32_mul (vector_a (1), vector_b (1))
               do iarr = 2, larr
                  uwrk = u32_add (uwrk, u32_mul (vector_a(iarr),        &
     &                                           vector_b(iarr)))
               enddo
               u32_dot = uwrk
            endif
         end function u32_dot
!
! _______product_______________________________________________________
!
         function u32_prd (array)
            type (uint_32), intent (in), dimension (:) :: array
            type (uint_32)                             :: u32_prd
            type (uint_32)                             :: uwrk
            integer                                    :: iarr, larr
!
            larr = size (array)
            if (larr <= 0) then
               u32_prd = uint_32 (1)
            else
               uwrk = array (1)
               do iarr = 2, larr
                  uwrk = u32_mul (uwrk, array (iarr))
               enddo
               u32_prd = uwrk
            endif
         end function u32_prd
!
! _______sum___________________________________________________________
!
         function u32_sum (array)
            type (uint_32), intent (in), dimension (:) :: array
            type (uint_32)                             :: u32_sum
            type (uint_32)                             :: uwrk
            integer                                    :: iarr, larr
!
            larr = size (array)
            if (larr <= 0) then
               u32_sum = uint_32 (0)
            else
               uwrk = array (1)
               do iarr = 2, larr
                  uwrk = u32_add (uwrk, array (iarr))
               enddo
               u32_sum = uwrk
            endif
         end function u32_sum
!
! _______matmul________________________________________________________
!
         function u32_mmm (matrix_a, matrix_b)
            type (uint_32), intent (in), dimension (:, :) :: matrix_a
            type (uint_32), intent (in), dimension (:, :) :: matrix_b
            type (uint_32), dimension (1:size (matrix_a, 1),            &
     &                              1:size (matrix_b, 2)) :: u32_mmm
            integer                                       :: i, j
!
            do j = 1, size (matrix_b, 2)
               do i = 1, size (matrix_a, 1)
                  u32_mmm (i, j) = u32_dot (matrix_a (i, :),            &
     &                                      matrix_b (:, j))
               enddo
            enddo
         end function u32_mmm
!
         function u32_mmv (matrix_a, matrix_b)
            type (uint_32), intent (in), dimension (:, :)  :: matrix_a
            type (uint_32), intent (in), dimension (:)     :: matrix_b
            type (uint_32), dimension (size (matrix_a, 1)) :: u32_mmv
            integer                                        :: i
!
            do i = 1, size (matrix_a, 1)
               u32_mmv (i) = u32_dot (matrix_a (i, :), matrix_b (:))
            enddo
         end function u32_mmv
!
         function u32_vmm (matrix_a, matrix_b)
            type (uint_32), intent (in), dimension (:)     :: matrix_a
            type (uint_32), intent (in), dimension (:, :)  :: matrix_b
            type (uint_32), dimension (size (matrix_b, 2)) :: u32_vmm
            integer                                        :: j
!
            do j = 1, size (matrix_b, 2)
               u32_vmm (j) = u32_dot (matrix_a (:), matrix_b (:, j))
            enddo
         end function u32_vmm
!
      end module unsigned_32
      program tstu32
      use unsigned_32
      type (uint_32), dimension (1:2)    :: u, v
      character (len=11) :: zstr
      integer            :: i (2)
      i = (/ -1, 1 /)
      write (*,*)
      u (:) = i (:)
      v (:) = u (1)
      v = dim (v, u)
      v = abs (v)
      zstr = v (2)
      write (*,*) zstr
      zstr = u (2)
      write (*,*) zstr
      u (2)=uint (2.0d0)**31 + (uint (2) * uint (2 ** 30) - uint_32 (1))
      zstr = u (2)
      write (*,*) zstr
      stop
      end program tstu32
