count_set_bits Function

    elemental function count_set_bits(a) result(nbits)
        integer(n_int), intent(in) :: a
        integer :: nbits
        integer(n_int) :: tmp

#ifdef INT64_
        integer(n_int), parameter :: m1 = 6148914691236517205_n_int  !Z'5555555555555555'
        integer(n_int), parameter :: m2 = 3689348814741910323_n_int  !Z'3333333333333333'
        integer(n_int), parameter :: m3 = 1085102592571150095_n_int  !Z'0f0f0f0f0f0f0f0f'
        integer(n_int), parameter :: m4 = 72340172838076673_n_int    !Z'0101010101010101'

        ! For 64 bit integers:
        tmp = a - iand(ishft(a, -1), m1)
        tmp = iand(tmp, m2) + iand(ishft(tmp, -2), m2)
        tmp = iand(tmp, m3) + iand(ishft(tmp, -4), m3)
        nbits = int(ishft(tmp * m4, -56))
#else
        integer(n_int), parameter :: m1 = 1431655765_n_int    !Z'55555555'
        integer(n_int), parameter :: m2 = 858993459_n_int     !Z'33333333'
        integer(n_int), parameter :: m3 = 252645135_n_int     !Z'0F0F0F0F'
        integer(n_int), parameter :: m4 = 16843009_n_int      !Z'01010101'

        ! For 32 bit integers:
        tmp = a - iand(ishft(a, -1), m1)
        tmp = iand(tmp, m2) + iand(ishft(tmp, -2), m2)
        tmp = iand((tmp + ishft(tmp, -4)), m3) * m4
        nbits = ishft(tmp, -24)
#endif

    end function count_set_bits