submodule(stdlib_sorting) stdlib_sorting_radix_sort implicit none integer, parameter :: radix_bits = 8 integer, parameter :: radix_mask = 255 integer(kind=int16), parameter :: radix_bits_i16 = 8_int16 integer(kind=int16), parameter :: radix_mask_i16 = 255_int16 integer(kind=int32), parameter :: radix_bits_i32 = 8_int32 integer(kind=int32), parameter :: radix_mask_i32 = 255_int32 integer(kind=int64), parameter :: radix_bits_i64 = 8_int64 integer(kind=int64), parameter :: radix_mask_i64 = 255_int64 contains ! For int8, radix sort becomes counting sort, so buffer is not needed pure subroutine radix_sort_u8_helper(N, arr) integer(kind=int_index), intent(in) :: N integer(kind=int8), dimension(N), intent(inout) :: arr integer(kind=int_index) :: i integer :: bin_idx integer(kind=int_index), dimension(-128:127) :: counts counts(:) = 0 do i = 1, N bin_idx = arr(i) counts(bin_idx) = counts(bin_idx) + 1 end do i = 1 do bin_idx = -128, 127 do while (counts(bin_idx) > 0) arr(i) = int(bin_idx, kind=int8) i = i + 1 counts(bin_idx) = counts(bin_idx) - 1 end do end do end subroutine pure subroutine radix_sort_u16_helper(N, arr, buf) integer(kind=int_index), intent(in) :: N integer(kind=int16), dimension(N), intent(inout) :: arr integer(kind=int16), dimension(N), intent(inout) :: buf integer(kind=int_index) :: i integer :: b, b0, b1 integer(kind=int_index), dimension(0:radix_mask) :: c0, c1 c0(:) = 0 c1(:) = 0 do i = 1, N b0 = iand(arr(i), radix_mask_i16) b1 = ishft(arr(i), -radix_bits_i16) c0(b0) = c0(b0) + 1 c1(b1) = c1(b1) + 1 end do do b = 1, radix_mask c0(b) = c0(b) + c0(b - 1) c1(b) = c1(b) + c1(b - 1) end do do i = N, 1, -1 b0 = iand(arr(i), radix_mask_i16) buf(c0(b0)) = arr(i) c0(b0) = c0(b0) - 1 end do do i = N, 1, -1 b1 = ishft(buf(i), -radix_bits_i16) arr(c1(b1)) = buf(i) c1(b1) = c1(b1) - 1 end do end subroutine pure subroutine radix_sort_u32_helper(N, arr, buf) integer(kind=int_index), intent(in) :: N integer(kind=int32), dimension(N), intent(inout) :: arr integer(kind=int32), dimension(N), intent(inout) :: buf integer(kind=int_index) :: i integer :: b, b0, b1, b2, b3 integer(kind=int_index), dimension(0:radix_mask) :: c0, c1, c2, c3 c0(:) = 0 c1(:) = 0 c2(:) = 0 c3(:) = 0 do i = 1, N b0 = iand(arr(i), radix_mask_i32) b1 = iand(ishft(arr(i), -radix_bits_i32), radix_mask_i32) b2 = iand(ishft(arr(i), -2*radix_bits_i32), radix_mask_i32) b3 = ishft(arr(i), -3*radix_bits_i32) c0(b0) = c0(b0) + 1 c1(b1) = c1(b1) + 1 c2(b2) = c2(b2) + 1 c3(b3) = c3(b3) + 1 end do do b = 1, radix_mask c0(b) = c0(b) + c0(b - 1) c1(b) = c1(b) + c1(b - 1) c2(b) = c2(b) + c2(b - 1) c3(b) = c3(b) + c3(b - 1) end do do i = N, 1, -1 b0 = iand(arr(i), radix_mask_i32) buf(c0(b0)) = arr(i) c0(b0) = c0(b0) - 1 end do do i = N, 1, -1 b1 = iand(ishft(buf(i), -radix_bits_i32), radix_mask_i32) arr(c1(b1)) = buf(i) c1(b1) = c1(b1) - 1 end do do i = N, 1, -1 b2 = iand(ishft(arr(i), -2*radix_bits_i32), radix_mask_i32) buf(c2(b2)) = arr(i) c2(b2) = c2(b2) - 1 end do do i = N, 1, -1 b3 = ishft(buf(i), -3*radix_bits_i32) arr(c3(b3)) = buf(i) c3(b3) = c3(b3) - 1 end do end subroutine radix_sort_u32_helper pure subroutine radix_sort_u64_helper(N, arr, buffer) integer(kind=int_index), intent(in) :: N integer(kind=int64), dimension(N), intent(inout) :: arr integer(kind=int64), dimension(N), intent(inout) :: buffer integer(kind=int_index) :: i integer(kind=int64) :: b, b0, b1, b2, b3, b4, b5, b6, b7 integer(kind=int_index), dimension(0:radix_mask) :: c0, c1, c2, c3, c4, c5, c6, c7 c0(:) = 0 c1(:) = 0 c2(:) = 0 c3(:) = 0 c4(:) = 0 c5(:) = 0 c6(:) = 0 c7(:) = 0 do i = 1, N b0 = iand(arr(i), radix_mask_i64) b1 = iand(ishft(arr(i), -radix_bits_i64), radix_mask_i64) b2 = iand(ishft(arr(i), -2*radix_bits_i64), radix_mask_i64) b3 = iand(ishft(arr(i), -3*radix_bits_i64), radix_mask_i64) b4 = iand(ishft(arr(i), -4*radix_bits_i64), radix_mask_i64) b5 = iand(ishft(arr(i), -5*radix_bits_i64), radix_mask_i64) b6 = iand(ishft(arr(i), -6*radix_bits_i64), radix_mask_i64) b7 = ishft(arr(i), -7*radix_bits_i64) c0(b0) = c0(b0) + 1 c1(b1) = c1(b1) + 1 c2(b2) = c2(b2) + 1 c3(b3) = c3(b3) + 1 c4(b4) = c4(b4) + 1 c5(b5) = c5(b5) + 1 c6(b6) = c6(b6) + 1 c7(b7) = c7(b7) + 1 end do do b = 1, radix_mask c0(b) = c0(b) + c0(b - 1) c1(b) = c1(b) + c1(b - 1) c2(b) = c2(b) + c2(b - 1) c3(b) = c3(b) + c3(b - 1) c4(b) = c4(b) + c4(b - 1) c5(b) = c5(b) + c5(b - 1) c6(b) = c6(b) + c6(b - 1) c7(b) = c7(b) + c7(b - 1) end do do i = N, 1, -1 b0 = iand(arr(i), radix_mask_i64) buffer(c0(b0)) = arr(i) c0(b0) = c0(b0) - 1 end do do i = N, 1, -1 b1 = iand(ishft(buffer(i), -radix_bits_i64), radix_mask_i64) arr(c1(b1)) = buffer(i) c1(b1) = c1(b1) - 1 end do do i = N, 1, -1 b2 = iand(ishft(arr(i), -2*radix_bits_i64), radix_mask_i64) buffer(c2(b2)) = arr(i) c2(b2) = c2(b2) - 1 end do do i = N, 1, -1 b3 = iand(ishft(buffer(i), -3*radix_bits_i64), radix_mask_i64) arr(c3(b3)) = buffer(i) c3(b3) = c3(b3) - 1 end do do i = N, 1, -1 b4 = iand(ishft(arr(i), -4*radix_bits_i64), radix_mask_i64) buffer(c4(b4)) = arr(i) c4(b4) = c4(b4) - 1 end do do i = N, 1, -1 b5 = iand(ishft(buffer(i), -5*radix_bits_i64), radix_mask_i64) arr(c5(b5)) = buffer(i) c5(b5) = c5(b5) - 1 end do do i = N, 1, -1 b6 = iand(ishft(arr(i), -6*radix_bits_i64), radix_mask_i64) buffer(c6(b6)) = arr(i) c6(b6) = c6(b6) - 1 end do do i = N, 1, -1 b7 = ishft(buffer(i), -7*radix_bits_i64) arr(c7(b7)) = buffer(i) c7(b7) = c7(b7) - 1 end do end subroutine radix_sort_u64_helper pure module subroutine int8_radix_sort(array, reverse) integer(kind=int8), dimension(:), intent(inout) :: array logical, intent(in), optional :: reverse integer(kind=int8) :: item integer(kind=int_index) :: i, N N = size(array, kind=int_index) call radix_sort_u8_helper(N, array) if (optval(reverse, .false.)) then do i = 1, N/2 item = array(i) array(i) = array(N - i + 1) array(N - i + 1) = item end do end if end subroutine int8_radix_sort pure module subroutine int16_radix_sort(array, work, reverse) integer(kind=int16), dimension(:), intent(inout) :: array integer(kind=int16), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse integer(kind=int_index) :: i, N, start, middle, end integer(kind=int16), dimension(:), pointer :: buffer integer(kind=int16) :: item logical :: use_internal_buffer N = size(array, kind=int_index) if (present(work)) then if (size(work, kind=int_index) < N) then error stop "int16_radix_sort: work array is too small." end if use_internal_buffer = .false. buffer => work else use_internal_buffer = .true. allocate (buffer(N)) end if call radix_sort_u16_helper(N, array, buffer) ! After calling `radix_sort_u<width>_helper. The array is sorted as unsigned integers. ! In the view of signed array, the negative numbers are sorted but in the tail of the array. ! Use binary search to find the boundary, and move them to correct position. if (array(1) >= 0 .and. array(N) < 0) then start = 1 end = N middle = (1 + N)/2 do while (.true.) if (array(middle) >= 0) then start = middle + 1 else end = middle end if middle = (start + end)/2 if (start == end) exit end do buffer(1:(N - middle + 1)) = array(middle:N) buffer(N - middle + 2:N) = array(1:middle - 1) array(:) = buffer(:) end if if (optval(reverse, .false.)) then do i = 1, N/2 item = array(i) array(i) = array(N - i + 1) array(N - i + 1) = item end do end if if (use_internal_buffer) then deallocate (buffer) end if end subroutine int16_radix_sort pure module subroutine int32_radix_sort(array, work, reverse) integer(kind=int32), dimension(:), intent(inout) :: array integer(kind=int32), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse integer(kind=int_index) :: i, N, start, middle, end integer(kind=int32), dimension(:), pointer :: buffer integer(kind=int32) :: item logical :: use_internal_buffer N = size(array, kind=int_index) if (present(work)) then if (size(work, kind=int_index) < N) then error stop "int32_radix_sort: work array is too small." end if use_internal_buffer = .false. buffer => work else use_internal_buffer = .true. allocate (buffer(N)) end if call radix_sort_u32_helper(N, array, buffer) if (array(1) >= 0 .and. array(N) < 0) then start = 1 end = N middle = (1 + N)/2 do while (.true.) if (array(middle) >= 0) then start = middle + 1 else end = middle end if middle = (start + end)/2 if (start == end) exit end do buffer(1:(N - middle + 1)) = array(middle:N) buffer(N - middle + 2:N) = array(1:middle - 1) array(:) = buffer(:) end if if (optval(reverse, .false.)) then do i = 1, N/2 item = array(i) array(i) = array(N - i + 1) array(N - i + 1) = item end do end if if (use_internal_buffer) then deallocate (buffer) end if end subroutine int32_radix_sort module subroutine sp_radix_sort(array, work, reverse) use iso_c_binding real(kind=sp), dimension(:), intent(inout), target :: array real(kind=sp), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse integer(kind=int_index) :: i, N, pos, rev_pos integer(kind=int32), dimension(:), pointer :: arri32 integer(kind=int32), dimension(:), pointer :: buffer real(kind=sp) :: item logical :: use_internal_buffer N = size(array, kind=int_index) if (present(work)) then if (size(work, kind=int_index) < N) then error stop "sp_radix_sort: work array is too small." end if use_internal_buffer = .false. call c_f_pointer(c_loc(work), buffer, [N]) else use_internal_buffer = .true. allocate (buffer(N)) end if call c_f_pointer(c_loc(array), arri32, [N]) call radix_sort_u32_helper(N, arri32, buffer) ! After calling `radix_sort_u<width>_helper. The array is sorted as unsigned integers. ! The positive real number is sorted, guaranteed by IEEE-754 standard. ! But the negative real number is sorted in a reversed order, and also in the tail of array. ! Remark that -0.0 is the minimum nagative integer, so using radix sort, -0.0 is naturally lesser than 0.0. ! In IEEE-754 standard, the bit representation of `Inf` is greater than all positive real numbers, ! and the `-Inf` is lesser than all negative real numbers. So the order of `Inf, -Inf` is ok. ! The bit representation of `NaN` may be positive or negative integers in different machine, ! thus if the array contains `NaN`, the result is undefined. if (arri32(1) >= 0 .and. arri32(N) < 0) then pos = 1 rev_pos = N do while (arri32(rev_pos) < 0) buffer(pos) = arri32(rev_pos) pos = pos + 1 rev_pos = rev_pos - 1 end do buffer(pos:N) = arri32(1:rev_pos) arri32(:) = buffer(:) end if if (optval(reverse, .false.)) then do i = 1, N/2 item = array(i) array(i) = array(N - i + 1) array(N - i + 1) = item end do end if if (use_internal_buffer) then deallocate (buffer) end if end subroutine sp_radix_sort pure module subroutine int64_radix_sort(array, work, reverse) integer(kind=int64), dimension(:), intent(inout) :: array integer(kind=int64), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse integer(kind=int_index) :: i, N, start, middle, end integer(kind=int64), dimension(:), pointer :: buffer integer(kind=int64) :: item logical :: use_internal_buffer N = size(array, kind=int_index) if (present(work)) then if (size(work, kind=int_index) < N) then error stop "int64_radix_sort: work array is too small." end if use_internal_buffer = .false. buffer => work else use_internal_buffer = .true. allocate (buffer(N)) end if call radix_sort_u64_helper(N, array, buffer) if (array(1) >= 0 .and. array(N) < 0) then start = 1 end = N middle = (1 + N)/2 do while (.true.) if (array(middle) >= 0) then start = middle + 1 else end = middle end if middle = (start + end)/2 if (start == end) exit end do buffer(1:(N - middle + 1)) = array(middle:N) buffer(N - middle + 2:N) = array(1:middle - 1) array(:) = buffer(:) end if if (optval(reverse, .false.)) then do i = 1, N/2 item = array(i) array(i) = array(N - i + 1) array(N - i + 1) = item end do end if if (use_internal_buffer) then deallocate (buffer) end if end subroutine int64_radix_sort module subroutine dp_radix_sort(array, work, reverse) use iso_c_binding real(kind=dp), dimension(:), intent(inout), target :: array real(kind=dp), dimension(:), intent(inout), target, optional :: work logical, intent(in), optional :: reverse integer(kind=int_index) :: i, N, pos, rev_pos integer(kind=int64), dimension(:), pointer :: arri64 integer(kind=int64), dimension(:), pointer :: buffer real(kind=dp) :: item logical :: use_internal_buffer N = size(array, kind=int_index) if (present(work)) then if (size(work, kind=int_index) < N) then error stop "sp_radix_sort: work array is too small." end if use_internal_buffer = .false. call c_f_pointer(c_loc(work), buffer, [N]) else use_internal_buffer = .true. allocate (buffer(N)) end if call c_f_pointer(c_loc(array), arri64, [N]) call radix_sort_u64_helper(N, arri64, buffer) if (arri64(1) >= 0 .and. arri64(N) < 0) then pos = 1 rev_pos = N do while (arri64(rev_pos) < 0) buffer(pos) = arri64(rev_pos) pos = pos + 1 rev_pos = rev_pos - 1 end do buffer(pos:N) = arri64(1:rev_pos) arri64(:) = buffer(:) end if if (optval(reverse, .false.)) then do i = 1, N/2 item = array(i) array(i) = array(N - i + 1) array(N - i + 1) = item end do end if if (use_internal_buffer) then deallocate (buffer) end if end subroutine dp_radix_sort end submodule stdlib_sorting_radix_sort