intT n, intT m, F extract) {

for (intT i = 0; i < m; i++) counts[i] = 0;

for (intT j = 0; j < n; j++) {

intT k = Tmp[j] = extract(A[j]);

counts[k]++;

}

intT s = 0;

for (intT i = 0; i < m; i++) {

s += counts[i];

counts[i] = s;

}

for (intT j = n-1; j >= 0; j--) {

intT x = --counts[Tmp[j]];

B[x] = A[j];

}

intT _mask; intT _offset; F _f;

eBits(int bits, intT offset, F f): _mask((1<<bits)-1),

_offset(offset), _f(f) {}

intT operator() (E p) {return _mask&(_f(p)>>_offset);}

intT a=0;

while ((1L << a) <= i) a++;

return a;

intT bits = log2Up(m);

// temporary space

E* B = (E*) malloc(sizeof(E)*n);

bIndexT* Tmp = (bIndexT*) malloc(sizeof(bIndexT)*n);

intT* counts = (intT*) malloc(sizeof(intT)*BUCKETS);

intT rounds = 1+(bits-1)/MAX_RADIX;

intT rbits = 1+(bits-1)/rounds;

intT bitOffset = 0;

bool flipped = 0;

while (bitOffset < bits) {

if (bitOffset+rbits > bits) rbits = bits-bitOffset;

if (flipped)

radixStep(B, A, Tmp, counts, n, 1L << rbits,

eBits<E,F>(rbits,bitOffset,f));

else

radixStep(A, B, Tmp, counts, n, 1L << rbits,

eBits<E,F>(rbits,bitOffset,f));

bitOffset += rbits;

flipped = !flipped;

}

if (flipped)

for (intT i=0; i < n; i++)

A[i] = B[i];

free(B); free(Tmp); free(counts);