radixsort.f90

! ===================================================================
! | ©2026 Z. Grof (grofz@vscht.cz)                                  |
! | Licensed under MIT licence. Provided "as is", without warranty. |
! ===================================================================

module radixsort_module
  implicit none (type, external)
  private
  public radixsort, radixsort_msb

  integer, parameter :: BASE = 2**8

contains

  pure subroutine radixsort(A)
    integer, intent(inout) :: A(:)
!
! LSB radix sort algorithm for an array of integers.
!
    integer, allocatable, dimension(:) :: new_nums, nums, new_rems, rems, buckets
    integer, dimension(0:BASE-1) :: bucket_ids, bucket_counts ! 0-based indexing
    integer :: n, i, ib, ipass

    n = size(A)
    allocate(new_nums(n), nums(n), new_rems(n), rems(n), buckets(n))
    new_nums = A
    new_rems = abs(A) ! negative numbers will be dealt at the end

    MAINLOOP: do ipass = 1, max_passes(new_rems)
      call swap_alloc(new_nums, nums)
      call swap_alloc(new_rems, rems)
      ! strip last "digit" and assign item to a particular bucket
      call strip_last_digit(rems, buckets)
      ! count number of items in every bucket
      bucket_counts = 0
      do i=1, n
        bucket_counts(buckets(i)) = bucket_counts(buckets(i)) + 1
      end do
      ! starting position for each bucket
      bucket_ids(0) = 1
      do ib=1, BASE-1
        bucket_ids(ib) = bucket_counts(ib-1) + bucket_ids(ib-1)
      end do
      ! copy items to their new position in the array
      do i=1, n
        associate(b => bucket_ids(buckets(i)))
          b = bucket_ids(buckets(i))
          new_nums(b) = nums(i)
          new_rems(b) = rems(i)
          b = b + 1
        end associate
      end do
    end do MAINLOOP

    ! copy back to original array by unwinding negative numbers 
    block
      integer :: positive_id, negative_id
      negative_id = count(new_nums < 0) 
      positive_id = negative_id + 1
      do i=1, n
        if (new_nums(i) < 0) then
          A(negative_id) = new_nums(i)
          negative_id = negative_id - 1
        else
          A(positive_id) = new_nums(i)
          positive_id = positive_id + 1
        end if
      end do
    end block
  end subroutine radixsort


  pure subroutine swap_alloc(a, b)
    integer, intent(inout), allocatable :: a(:), b(:)
    integer, allocatable :: tmp(:)
    call move_alloc(a, tmp)
    call move_alloc(b, a)
    call move_alloc(tmp, b)
  end subroutine swap_alloc


  elemental subroutine strip_last_digit(number, last_digit)
    integer, intent(inout) :: number
    integer, intent(out) :: last_digit
!
! Number must be non-negative. Remove least significant "digit" that determines
! the assignment to the bucket.
!
    last_digit = mod(number, BASE)
    number = number/BASE
  end subroutine strip_last_digit


  pure integer function max_passes(arr)
    integer, intent(in) :: arr(:)
!
! Using the largest value, determine the required number of passes
!
    integer :: number, digit
    number = maxval(arr)
    max_passes = 0
    do while(number /= 0)
      call strip_last_digit(number, digit)
      max_passes = max_passes + 1
    end do
  end function max_passes


  pure subroutine radixsort_msb(A)
    integer, intent(inout) :: A(:)
!
! MSB radix sort implementation
!
    integer, allocatable :: tmpA(:), buckets(:), rems(:), tmp_rems(:)
    integer :: n, i, maxorder, positive_id, negative_id, negative_count

    ! working space for recursive subroutine
    n = size(A)
    allocate(tmpA(n), buckets(n), rems(n), tmp_rems(n))

    ! invert negative numbers (abs(x)) and move them before non-negative ones
    negative_count = count(A < 0)
    negative_id = 1
    positive_id = negative_count + 1
    do i=1, n
      if (A(i) < 0) then
        tmpA(negative_id) = abs(A(i))
        negative_id = negative_id + 1
      else
        tmpA(positive_id) = A(i)
        positive_id = positive_id + 1
      end if
    end do

    ! remainder as MSB digits are stripped
    rems = tmpA

    ! identify the largest MSB "digit"
    maxorder = max_passes(tmpA)

    ! negatives and positives are sorted separately
    call rxsort(tmpA, A, buckets, rems, tmp_rems, 1, negative_count, maxorder)
    call rxsort(tmpA, A, buckets, rems, tmp_rems, negative_count+1, n, maxorder)

    ! copy to original array, invert order of negative numbers
    A(1:negative_count) = -tmpA(negative_count:1:-1)
    A(negative_count+1:n) = tmpA(negative_count+1:n)
  end subroutine radixsort_msb


  pure recursive subroutine rxsort(nums, tmp_nums, buckets, rems, tmp_rems, &
                                   from, to, order)
    integer, intent(inout) :: nums(:), rems(:)
    integer, intent(out) :: tmp_nums(:), tmp_rems(:), buckets(:) ! temporary working space
    integer, intent(in) :: from, to, order

    integer :: i, ib
    integer, dimension(0:BASE) :: bucket_ids, bucket_counts, bucket_ids_saved
      ! note zero-based indexing

    if (from >= to) return ! nothing to be done with one or zero items

    ! strip the actual digit and count number of items in every bucket
    bucket_counts = 0
    call strip_msb_digit(rems(from:to), buckets(from:to), order)
    do i=from, to
      bucket_counts(buckets(i)) = bucket_counts(buckets(i)) + 1
    end do
    ! range assigned to buckets
    bucket_ids(0) = from
    do ib=1, BASE
      bucket_ids(ib) = bucket_ids(ib-1)+bucket_counts(ib-1)
    end do
    bucket_ids_saved = bucket_ids
    ! move items to its new positions
    do i=from, to
      associate(b => bucket_ids(buckets(i)))
        tmp_nums(b) = nums(i)
        tmp_rems(b) = rems(i)
        b = b + 1
      end associate
    end do
    ! move back from temporary locations
    rems(from:to) = tmp_rems(from:to)
    nums(from:to) = tmp_nums(from:to)

    ! recursively sort items in every bucket
    if (order > 1) then
      do ib=0,BASE-1
        call rxsort(nums, tmp_nums, buckets, rems, tmp_rems, &
                    bucket_ids_saved(ib), bucket_ids_saved(ib+1)-1, order-1)
      end do
    end if
  end subroutine rxsort


  elemental subroutine strip_msb_digit(number, msb_digit, order)
    integer, intent(inout) :: number
    integer, intent(out) :: msb_digit
    integer, intent(in) :: order
!
! Number must be non-negative. Remove most significant "digit" that determines
! the assignment to the bucket.
!
    msb_digit = number / (BASE**(order-1))
    number = number - msb_digit*BASE**(order-1)
  end subroutine strip_msb_digit

end module radixsort_module
!EOF: radixsort.f90