*-----------------------------------------------------------------------
* rsorti - radix sort and carry along an integer array
* by Andy Allinger, 2011-2016, released to the public domain
* This program may be used by any person for any purpose
*
* origin:  Herman Hollerith, 1887
*
* NOTE:  This routine is a sort-and-carry.  To obtain a permutation 
* vector to apply to other arrays, you must first initialize IND to 
* 1, 2, 3...N
*
* NOTE ALSO:
*   if B is positive, then the integers will be interpreted as UNSIGNED
*   note that FORTRAN does not support unsigned integers, so make sure
*   you know what you are doing if you use B=32
*   if B is negative, the sign bit will be inspected as a final step
*   So, for the range -256...255 set B = -8
*
*__Name____Type_____In/Out_________Description__________________________
*  D[N]    Integer  Both           Array, returned sorted
*  IND[N]  Integer  Both           Array to carry along, meaning
*                                         permuted in the same fashion as D
*  N       Integer  In             Length of the array
*  B       Integer  In             Number of significant bits
*  W[2@N]  Integer  Neither        Workspace
*-----------------------------------------------------------------------
      SUBROUTINE RSORTI (D, IND, N, B, W)    ! partition workspace
       IMPLICIT NONE
       INTEGER D, IND, N, B, W
       DIMENSION D(N), IND(N), W(2*N)
       
       CALL RSI01 (D, IND, N, B, W(1), W(N+1))
       RETURN
      END  ! of RSORTI


*-----------------------------------------------------------------------
      SUBROUTINE RSI01 (D, D2, N, B, W, W2)
      
       IMPLICIT NONE
       INTEGER D, D2, N, B, W, W2
       DIMENSION D(N), D2(N), W(N), W2(N)

       INTEGER I, J, K, ILIM, P1OLD, P0OLD, P1, P0, S
       LOGICAL INW                         ! data is in W
       INTEGER BITSZ                       ! external function

*          validate
       IF (B .EQ. 0) STOP 'Bit size must not be zero'
       IF (ABS(B) > BITSZ(N)) STOP 'Bit size to large'

*-----------------------------------------------------------------------
*      sort starts on the least significant bit, position 0
*-----------------------------------------------------------------------
       P1 = N+1
       P0 = 0
       DO J = 1, N

*      test for negative values
         IF (D(J) < 0 .AND. B > 0 .AND. B .NE. BITSZ(N)) 
     &       STOP 'Set argument B negative to sort negative numbers'

         IF ( BTEST(D(J), 0) ) THEN
           P1 = P1 - 1
           W(P1) = D(J)
           W2(P1) = D2(J)
         ELSE 
           P0 = P0 + 1
           W(P0) = D(J)
           W2(P0) = D2(J)
         END IF
       END DO
       INW = .TRUE.

*-----------------------------------------------------------------------
*      now alternate between putting data into D and putting data into W
*-----------------------------------------------------------------------
       ILIM = 0
       IF (B < 0) ILIM = MIN( ABS(B)-1, BITSZ(I)-2 )
       IF (B > 0) ILIM = B - 1
       DO I = 1, ILIM
         P1OLD = P1
         P0OLD = P0         ! save the value from previous bit
         P1 = N+1
         P0 = 0                 ! start a fresh count for next bit

*         odd I
         IF (INW) THEN

*          copy data from the zeros
           DO J = 1, P0OLD, 1
             IF ( BTEST(W(J), I) ) THEN
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             ELSE
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             END IF
           END DO
          
*          copy data from the ones
           DO J = N, P1OLD, -1           
             IF ( BTEST(W(J), I) ) THEN
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             ELSE
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             END IF
           END DO
          
*        even I
         ELSE 

*          copy data from the zeros
           DO J = 1, P0OLD, 1
             IF ( BTEST(D(J), I) ) THEN
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             ELSE
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             END IF
           END DO

*            copy data from the ones
           DO J = N, P1OLD, -1
             IF ( BTEST(D(J), I) ) THEN
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             ELSE
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             END IF
           END DO
         END IF ! even or odd i
        
         INW = .NOT. INW
       END DO   ! next i

*-----------------------------------------------------------------------
*        the sign bit
*-----------------------------------------------------------------------
       IF (B < 0) THEN
         P1OLD = P1
         P0OLD = P0
         P1 = N+1
         P0 = 0 

*          if sign bit is set, send to the zero end of the work array
         IF (INW) THEN
           DO J = 1, P0OLD, 1
             IF ( BTEST(W(J), BITSZ(I)-1) ) THEN 
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             ELSE
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             END IF
           END DO          
           DO J = N, P1OLD, -1
             IF ( BTEST(W(J), BITSZ(I)-1) ) THEN
               P0 = P0 + 1
               D(P0) = W(J)
               D2(P0) = W2(J)
             ELSE
               P1 = P1 - 1
               D(P1) = W(J)
               D2(P1) = W2(J)
             END IF
           END DO
          
         ELSE
           DO J = 1, P0OLD, 1
             IF ( BTEST(D(J), BITSZ(I)-1) ) THEN 
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             ELSE
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             END IF
           END DO          
           DO J = N, P1OLD, -1
             IF ( BTEST(D(J), BITSZ(I)-1) ) THEN
               P0 = P0 + 1
               W(P0) = D(J)
               W2(P0) = D2(J)
             ELSE
               P1 = P1 - 1
               W(P1) = D(J)
               W2(P1) = D2(J)
             END IF
           END DO
         END IF

         INW = .NOT. INW  
       END IF
          
*-----------------------------------------------------------------------
*      put the data back into the input array
*-----------------------------------------------------------------------
       K = P1
       IF (INW) THEN     
         DO J = 1, P0, 1
           D(J) = W(J)
           D2(J) = W2(J)
         END DO
         DO J = N, P1, -1        
           D(K) = W(J)
           D2(K) = W2(J)
           K = K + 1
         END DO
       
*      Swap values in the ones portion
       ELSE 
         DO J = N, (N+P1+2)/2, -1
           S = D(J)
           D(J) = D(K)
           D(K) = S
           S = D2(J)
           D2(J) = D2(K)
           D2(K) = S
           K = K + 1
         END DO
       END IF

      END ! of RSI01