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.2 of 12/17/93 ! ! 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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 (jbits (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 = jbits (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 = jbits (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 (jbits (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 = jbits (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 = (jbits (u32a%val, 0, lh) == jbits (u32b%val, 0, lh))& & .and. (jbits (u32a%val, lh, lh) == jbits (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 = (jbits (u32a%val, 0, lh) /= jbits (u32b%val, 0, lh))& & .or. (jbits (u32a%val, lh, lh) /= jbits (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 = jbits (u32a%val, lh, lh) jwrk2 = jbits (u32b%val, lh, lh) if (jwrk1 == jwrk2) then u32_tgt = jbits (u32a%val, 0, lh) > & & jbits (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 = jbits (u32a%val, lh, lh) jwrk2 = jbits (u32b%val, lh, lh) if (jwrk1 == jwrk2) then u32_tlt = jbits (u32a%val, 0, lh) < & & jbits (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 = jbits (u32a%val, lh, lh) jwrk2 = jbits (u32b%val, lh, lh) if (jwrk1 == jwrk2) then u32_tge = jbits (u32a%val, 0, lh) >= & & jbits (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 = jbits (u32a%val, lh, lh) jwrk2 = jbits (u32b%val, lh, lh) if (jwrk1 == jwrk2) then u32_tle = jbits (u32a%val, 0, lh) <= & & jbits (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 = jbits (u32a%val, 0, lh) + jbits (u32b%val, 0, lh) jwrk2 = jbits (u32a%val, lh, lh) + jbits (u32b%val, lh, lh) & & + jbits (jwrk1, lh, lh) jwrk1 = jbits (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 (jbits (u32b%val, 0, lw)) jwrk1 = jbits (u32a%val, 0, lh) + jbits (jwrkb, 0, lh) + 1 jwrk2 = jbits (u32a%val, lh, lh) + jbits (jwrkb, lh, lh) + & & jbits (jwrk1, lh, lh) jwrk1 = jbits (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 = jbits (u32a%val, 0, lh) jwrkb1 = jbits (u32b%val, 0, lh) jwrka2 = jbits (u32a%val, lh, lh) jwrkb2 = jbits (u32b%val, lh, lh) jwrk1 = jwrka1 * jwrkb1 jwrk2 = jwrka1 * jwrkb2 + jwrka2 * jwrkb1 + & & jbits (jwrk1, lh, lh) jwrk1 = jbits (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 = jbits (u32a%val, 1, lwm1) / & & jbits (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 + jbits (uwrk%val, 0, lwm1) / & & jbits (u32b%val, 0, lwm1) u32_div = u32_add (uwrk1, uwrk2) endif endif ! ! ... a small, b small else u32_div%val = jbits (u32a%val, 0, lwm1) / & & jbits (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 = jbits (-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 = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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) = jbits (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) = jbits (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) = jbits (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) = jbits (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) = jbits (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 = jbits (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 = jbits (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 = jbits (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 = jbits (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) = & & (jbits (u32a(i)%val, 0, lh) == jbits (u32b(i)%val, 0, lh))& & .and.(jbits (u32a(i)%val, lh, lh) == jbits (u32b(i)%val, lh, lh)) enddo end function u32a_teq ! ! _______test for inequality____________________________________________ ! logical 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) = & & (jbits (u32a(i)%val, 0, lh) /= jbits (u32b(i)%val, 0, lh))& & .or. (jbits (u32a(i)%val, lh, lh) /= jbits (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 = jbits (u32a (ival)%val, lh, lh) jwrk2 = jbits (u32b (ival)%val, lh, lh) if (jwrk1 == jwrk2) then u32a_tgt (ival) = jbits (u32a (ival)%val, 0, lh) > & & jbits (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 = jbits (u32a (ival)%val, lh, lh) jwrk2 = jbits (u32b (ival)%val, lh, lh) if (jwrk1 == jwrk2) then u32a_tlt (ival) = jbits (u32a (ival)%val, 0, lh) < & & jbits (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 = jbits (u32a (ival)%val, lh, lh) jwrk2 = jbits (u32b (ival)%val, lh, lh) if (jwrk1 == jwrk2) then u32a_tge (ival) = jbits (u32a (ival)%val, 0, lh) >= & & jbits (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 = jbits (u32a (ival)%val, lh, lh) jwrk2 = jbits (u32b (ival)%val, lh, lh) if (jwrk1 == jwrk2) then u32a_tle (ival) = jbits (u32a (ival)%val, 0, lh) <= & & jbits (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 = jbits (u32a (ival)%val, 0, lh) + & & jbits (u32b (ival)%val, 0, lh) jwrk2 = jbits (u32a (ival)%val, lh, lh) + & & jbits (u32b (ival)%val, lh, lh) + & & jbits (jwrk1, lh, lh) jwrk1 = jbits (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 (jbits (u32b (ival)%val, 0, lw)) jwrk1 = jbits (u32a (ival)%val, 0, lh) + & & jbits (jwrkb, 0, lh) + 1 jwrk2 = jbits (u32a (ival)%val, lh, lh) + & & jbits (jwrkb, lh, lh) + jbits (jwrk1, lh, lh) jwrk1 = jbits (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 = jbits (u32a (ival)%val, 0, lh) jwrkb1 = jbits (u32b (ival)%val, 0, lh) jwrka2 = jbits (u32a (ival)%val, lh, lh) jwrkb2 = jbits (u32b (ival)%val, lh, lh) jwrk1 = jwrka1 * jwrkb1 jwrk2 = jwrka1 * jwrkb2 + jwrka2 * jwrkb1 + & & jbits (jwrk1, lh, lh) jwrk1 = jbits (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 = jbits (u32a (ival)%val, 1, lwm1) / & & jbits (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 + jbits (uwrk%val, 0, lwm1) / & & jbits (u32b (ival)%val, 0, lwm1) u32a_div (ival) = u32_add (uwrk1, uwrk2) endif endif ! ! ... a small, b small else u32a_div(ival)%val = jbits(u32a(ival)%val, 0, lwm1)/& & jbits(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 = jbits (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 ! ! ****** When bug corrected, remove jbits and substitute jbits by ibits ! everywhere else ********************************************** ! _____________________________________________________________________ function jbits (jval, pos, len) integer, intent (in) :: jval, pos, len integer (kind=kind(jval)) :: jbits ! ! *** avoid bug in Nag f90 intrinsic ******************************** if (len == lw .and. pos == 0) then if (btest (jval, lwm1)) then jbits = ibset (ibits (jval, 0, lwm1), lwm1) else jbits = ibits (jval, 0, lwm1) endif else jbits = ibits (jval, pos, len) endif end function jbits ! ! 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