*-----------------------------------------------------------------------
* rsorti - passive radix sort
* by Andy Allinger, 2011-2019, released to the public domain
* This program may be used by any person for any purpose
*
* origin:  Herman Hollerith, 1887
*
* NOTE:
*   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___________________________
*   IX(N)     Integer   Both      Array to sort
*   IND(N)    Integer   Out       Permutation vector
*   N         Integer   In        Length of the array
*   B         Integer   In        Number of significant bits
*   IW(N)     Integer   Neither   Workspace
*-----------------------------------------------------------------------
      SUBROUTINE RSORTI (IX, IND, N, B, IW)
       IMPLICIT NONE
       INTEGER IX, IND, N, B, IW
       DIMENSION IX(N), IND(N), IW(N)

       INTEGER BS, I, J, K, ILIM, P1OLD, P0OLD, P1, P0
       LOGICAL ODD
       INTEGER BITSZ                       ! external function

*         Validate
       BS = BITSZ()
       IF (B .EQ. 0) THEN
         PRINT *, 'rsorti:  Bit size must not be zero'
         RETURN
       END IF
       IF (ABS(B) > BS) THEN
         PRINT *, 'rsorti:  Bit size too large'
         RETURN
       END IF
       DO J = 1, N                     ! Test for negative values.
         IF (IX(J) < 0 .AND. B > 0 .AND. B .NE. BS) THEN
           PRINT *, 'Set argument B negative to sort negative numbers'
           RETURN
         END IF
       END DO

       DO J = 1, N                   ! initialize
         IND(J) = J
       END DO
       IF (N < 2) RETURN

*-----------------------------------------------------------------------
* Alternate between referencing by IND and referencing by IW
*-----------------------------------------------------------------------
       ILIM = 0
       IF (B < 0) ILIM = MIN( ABS(B)-1, BS-2 )  ! -31 and -32 treated the same
       IF (B > 0) ILIM = B - 1

       P1 = N+1
       P0 = N
       ODD = .FALSE.
       DO I = 0, ILIM
         P1OLD = P1
         P0OLD = P0             ! save the value from previous bit
         P1 = N+1
         P0 = 0                 ! start a fresh count for next bit

         IF (ODD) THEN
           DO J = 1, P0OLD, 1               ! copy data from the zeros
             IF ( BTEST(IX(IW(J)), I) ) THEN
               P1 = P1 - 1
               IND(P1) = IW(J)
             ELSE
               P0 = P0 + 1
               IND(P0) = IW(J)
             END IF
           END DO

           DO J = N, P1OLD, -1              ! copy data from the ones
             IF ( BTEST(IX(IW(J)), I) ) THEN
               P1 = P1 - 1
               IND(P1) = IW(J)
             ELSE
               P0 = P0 + 1
               IND(P0) = IW(J)
             END IF
           END DO

         ELSE
           DO J = 1, P0OLD, 1               ! copy data from the zeros
             IF ( BTEST(IX(IND(J)), I) ) THEN
               P1 = P1 - 1
               IW(P1) = IND(J)
             ELSE
               P0 = P0 + 1
               IW(P0) = IND(J)
             END IF
           END DO

           DO J = N, P1OLD, -1              ! copy data from the ones
             IF ( BTEST(IX(IND(J)), I) ) THEN
               P1 = P1 - 1
               IW(P1) = IND(J)
             ELSE
               P0 = P0 + 1
               IW(P0) = IND(J)
             END IF
           END DO
         END IF ! even or odd i

         ODD = .NOT. ODD
       END DO   ! next i

*-----------------------------------------------------------------------
*        the sign bit
*-----------------------------------------------------------------------
       IF (B < 0) THEN
         P1OLD = P1
         P0OLD = P0
         P1 = N+1
         P0 = 0
         IF (ODD) THEN
           DO J = 1, P0OLD, 1
             IF ( BTEST(IX(IW(J)), BS-1) ) THEN
               P0 = P0 + 1
               IND(P0) = IW(J)
             ELSE
               P1 = P1 - 1
               IND(P1) = IW(J)
             END IF
           END DO
           DO J = N, P1OLD, -1
             IF ( BTEST(IX(IW(J)), BS-1) ) THEN
               P0 = P0 + 1
               IND(P0) = IW(J)
             ELSE
               P1 = P1 - 1
               IND(P1) = IW(J)
             END IF
           END DO

         ELSE
           DO J = 1, P0OLD, 1
             IF ( BTEST(IX(IND(J)), BS-1) ) THEN
               P0 = P0 + 1
               IW(P0) = IND(J)
             ELSE
               P1 = P1 - 1
               IW(P1) = IND(J)
             END IF
           END DO
           DO J = N, P1OLD, -1
             IF ( BTEST(IX(IND(J)), BS-1) ) THEN
               P0 = P0 + 1
               IW(P0) = IND(J)
             ELSE
               P1 = P1 - 1
               IW(P1) = IND(J)
             END IF
           END DO
         END IF

         ODD = .NOT. ODD
       END IF

*-----------------------------------------------------------------------
* Put indices into IND
*-----------------------------------------------------------------------
       K = P1
       IF (ODD) THEN
         DO J = 1, P0, 1
           IND(J) = IW(J)
         END DO
         DO J = N, P1, -1
           IND(K) = IW(J)
           K = K + 1
         END DO

*      Swap values in the ones portion
       ELSE
         DO J = N, (N+P1+2)/2, -1
           I = IND(J)
           IND(J) = IND(K)
           IND(K) = I
           K = K + 1
         END DO
       END IF
       RETURN
      END  ! of rsorti


*-----------------------------------------------------------------------
* BIT_SIZE replacement for old compilers.
*-----------------------------------------------------------------------
      FUNCTION BITSZ ()
       IMPLICIT NONE
       INTEGER BITSZ
       INTEGER I, J

       J = 0
       J = NOT(J)
       BITSZ = 0
       DO I = 1, 1024
         J = ISHFT(J, -1)
         IF (.NOT. BTEST(J, 0)) THEN
           BITSZ = I
           RETURN
         END IF
       END DO
       RETURN      ! error!
      END  ! of bitsz