stdlib_sorting_radix_sort.f90

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