archival/bz/blocksort.c

/*
 * bzip2 is written by Julian Seward <jseward@bzip.org>.
 * Adapted for busybox by Denys Vlasenko <vda.linux@googlemail.com>.
 * See README and LICENSE files in this directory for more information.
 */

/*-------------------------------------------------------------*/
/*--- Block sorting machinery                               ---*/
/*---                                           blocksort.c ---*/
/*-------------------------------------------------------------*/

/* ------------------------------------------------------------------
This file is part of bzip2/libbzip2, a program and library for
lossless, block-sorting data compression.

bzip2/libbzip2 version 1.0.4 of 20 December 2006
Copyright (C) 1996-2006 Julian Seward <jseward@bzip.org>

Please read the WARNING, DISCLAIMER and PATENTS sections in the
README file.

This program is released under the terms of the license contained
in the file LICENSE.
------------------------------------------------------------------ */

/* #include "bzlib_private.h" */

#define mswap(zz1, zz2) \
{ \
	int32_t zztmp = zz1; \
	zz1 = zz2; \
	zz2 = zztmp; \
}

static
/* No measurable speed gain with inlining */
/* ALWAYS_INLINE */
void mvswap(uint32_t* ptr, int32_t zzp1, int32_t zzp2, int32_t zzn)
{
	while (zzn > 0) {
		mswap(ptr[zzp1], ptr[zzp2]);
		zzp1++;
		zzp2++;
		zzn--;
	}
}

static
ALWAYS_INLINE
int32_t mmin(int32_t a, int32_t b)
{
	return (a < b) ? a : b;
}


/*---------------------------------------------*/
/*--- Fallback O(N log(N)^2) sorting        ---*/
/*--- algorithm, for repetitive blocks      ---*/
/*---------------------------------------------*/

/*---------------------------------------------*/
static
inline
void fallbackSimpleSort(uint32_t* fmap,
		uint32_t* eclass,
		int32_t   lo,
		int32_t   hi)
{
	int32_t i, j, tmp;
	uint32_t ec_tmp;

	if (lo == hi) return;

	if (hi - lo > 3) {
		for (i = hi-4; i >= lo; i--) {
			tmp = fmap[i];
			ec_tmp = eclass[tmp];
			for (j = i+4; j <= hi && ec_tmp > eclass[fmap[j]]; j += 4)
				fmap[j-4] = fmap[j];
			fmap[j-4] = tmp;
		}
	}

	for (i = hi-1; i >= lo; i--) {
		tmp = fmap[i];
		ec_tmp = eclass[tmp];
		for (j = i+1; j <= hi && ec_tmp > eclass[fmap[j]]; j++)
			fmap[j-1] = fmap[j];
		fmap[j-1] = tmp;
	}
}


/*---------------------------------------------*/
#define fpush(lz,hz) { \
	stackLo[sp] = lz; \
	stackHi[sp] = hz; \
	sp++; \
}

#define fpop(lz,hz) { \
	sp--; \
	lz = stackLo[sp]; \
	hz = stackHi[sp]; \
}

#define FALLBACK_QSORT_SMALL_THRESH 10
#define FALLBACK_QSORT_STACK_SIZE   100

static
void fallbackQSort3(uint32_t* fmap,
		uint32_t* eclass,
		int32_t   loSt,
		int32_t   hiSt)
{
	int32_t unLo, unHi, ltLo, gtHi, n, m;
	int32_t sp, lo, hi;
	uint32_t med, r, r3;
	int32_t stackLo[FALLBACK_QSORT_STACK_SIZE];
	int32_t stackHi[FALLBACK_QSORT_STACK_SIZE];

	r = 0;

	sp = 0;
	fpush(loSt, hiSt);

	while (sp > 0) {
		AssertH(sp < FALLBACK_QSORT_STACK_SIZE - 1, 1004);

		fpop(lo, hi);
		if (hi - lo < FALLBACK_QSORT_SMALL_THRESH) {
			fallbackSimpleSort(fmap, eclass, lo, hi);
			continue;
		}

		/* Random partitioning.  Median of 3 sometimes fails to
		 * avoid bad cases.  Median of 9 seems to help but
		 * looks rather expensive.  This too seems to work but
		 * is cheaper.  Guidance for the magic constants
		 * 7621 and 32768 is taken from Sedgewick's algorithms
		 * book, chapter 35.
		 */
		r = ((r * 7621) + 1) % 32768;
		r3 = r % 3;
		if (r3 == 0)
			med = eclass[fmap[lo]];
		else if (r3 == 1)
			med = eclass[fmap[(lo+hi)>>1]];
		else
			med = eclass[fmap[hi]];

		unLo = ltLo = lo;
		unHi = gtHi = hi;

		while (1) {
			while (1) {
				if (unLo > unHi) break;
				n = (int32_t)eclass[fmap[unLo]] - (int32_t)med;
				if (n == 0) {
					mswap(fmap[unLo], fmap[ltLo]);
					ltLo++;
					unLo++;
					continue;
				};
				if (n > 0) break;
				unLo++;
			}
			while (1) {
				if (unLo > unHi) break;
				n = (int32_t)eclass[fmap[unHi]] - (int32_t)med;
				if (n == 0) {
					mswap(fmap[unHi], fmap[gtHi]);
					gtHi--; unHi--;
					continue;
				};
				if (n < 0) break;
				unHi--;
			}
			if (unLo > unHi) break;
			mswap(fmap[unLo], fmap[unHi]); unLo++; unHi--;
		}

		AssertD(unHi == unLo-1, "fallbackQSort3(2)");

		if (gtHi < ltLo) continue;

		n = mmin(ltLo-lo, unLo-ltLo); mvswap(fmap, lo, unLo-n, n);
		m = mmin(hi-gtHi, gtHi-unHi); mvswap(fmap, unLo, hi-m+1, m);

		n = lo + unLo - ltLo - 1;
		m = hi - (gtHi - unHi) + 1;

		if (n - lo > hi - m) {
			fpush(lo, n);
			fpush(m, hi);
		} else {
			fpush(m, hi);
			fpush(lo, n);
		}
	}
}

#undef fpush
#undef fpop
#undef FALLBACK_QSORT_SMALL_THRESH
#undef FALLBACK_QSORT_STACK_SIZE


/*---------------------------------------------*/
/* Pre:
 *	nblock > 0
 *	eclass exists for [0 .. nblock-1]
 *	((uint8_t*)eclass) [0 .. nblock-1] holds block
 *	ptr exists for [0 .. nblock-1]
 *
 * Post:
 *	((uint8_t*)eclass) [0 .. nblock-1] holds block
 *	All other areas of eclass destroyed
 *	fmap [0 .. nblock-1] holds sorted order
 *	bhtab[0 .. 2+(nblock/32)] destroyed
*/

#define       SET_BH(zz)  bhtab[(zz) >> 5] |= (1 << ((zz) & 31))
#define     CLEAR_BH(zz)  bhtab[(zz) >> 5] &= ~(1 << ((zz) & 31))
#define     ISSET_BH(zz)  (bhtab[(zz) >> 5] & (1 << ((zz) & 31)))
#define      WORD_BH(zz)  bhtab[(zz) >> 5]
#define UNALIGNED_BH(zz)  ((zz) & 0x01f)

static
void fallbackSort(uint32_t* fmap,
		uint32_t* eclass,
		uint32_t* bhtab,
		int32_t   nblock)
{
	int32_t ftab[257];
	int32_t ftabCopy[256];
	int32_t H, i, j, k, l, r, cc, cc1;
	int32_t nNotDone;
	int32_t nBhtab;
	uint8_t* eclass8 = (uint8_t*)eclass;

	/*
	 * Initial 1-char radix sort to generate
	 * initial fmap and initial BH bits.
	 */
	for (i = 0; i < 257;    i++) ftab[i] = 0;
	for (i = 0; i < nblock; i++) ftab[eclass8[i]]++;
	for (i = 0; i < 256;    i++) ftabCopy[i] = ftab[i];

	j = ftab[0];  /* bbox: optimized */
	for (i = 1; i < 257;    i++) {
		j += ftab[i];
		ftab[i] = j;
	}

	for (i = 0; i < nblock; i++) {
		j = eclass8[i];
		k = ftab[j] - 1;
		ftab[j] = k;
		fmap[k] = i;
	}

	nBhtab = 2 + ((uint32_t)nblock / 32); /* bbox: unsigned div is easier */
	for (i = 0; i < nBhtab; i++) bhtab[i] = 0;
	for (i = 0; i < 256; i++) SET_BH(ftab[i]);

	/*
	 * Inductively refine the buckets.  Kind-of an
	 * "exponential radix sort" (!), inspired by the
	 * Manber-Myers suffix array construction algorithm.
	 */

	/*-- set sentinel bits for block-end detection --*/
	for (i = 0; i < 32; i++) {
		SET_BH(nblock + 2*i);
		CLEAR_BH(nblock + 2*i + 1);
	}

	/*-- the log(N) loop --*/
	H = 1;
	while (1) {
		j = 0;
		for (i = 0; i < nblock; i++) {
			if (ISSET_BH(i))
				j = i;
			k = fmap[i] - H;
			if (k < 0)
				k += nblock;
			eclass[k] = j;
		}

		nNotDone = 0;
		r = -1;
		while (1) {

		/*-- find the next non-singleton bucket --*/
			k = r + 1;
			while (ISSET_BH(k) && UNALIGNED_BH(k))
				k++;
			if (ISSET_BH(k)) {
				while (WORD_BH(k) == 0xffffffff) k += 32;
				while (ISSET_BH(k)) k++;
			}
			l = k - 1;
			if (l >= nblock)
				break;
			while (!ISSET_BH(k) && UNALIGNED_BH(k))
				k++;
			if (!ISSET_BH(k)) {
				while (WORD_BH(k) == 0x00000000) k += 32;
				while (!ISSET_BH(k)) k++;
			}
			r = k - 1;
			if (r >= nblock)
				break;

			/*-- now [l, r] bracket current bucket --*/
			if (r > l) {
				nNotDone += (r - l + 1);
				fallbackQSort3(fmap, eclass, l, r);

				/*-- scan bucket and generate header bits-- */
				cc = -1;
				for (i = l; i <= r; i++) {
					cc1 = eclass[fmap[i]];
					if (cc != cc1) {
						SET_BH(i);
						cc = cc1;
					};
				}
			}
		}

		H *= 2;
		if (H > nblock || nNotDone == 0)
			break;
	}

	/*
	 * Reconstruct the original block in
	 * eclass8 [0 .. nblock-1], since the
	 * previous phase destroyed it.
	 */
	j = 0;
	for (i = 0; i < nblock; i++) {
		while (ftabCopy[j] == 0)
			j++;
		ftabCopy[j]--;
		eclass8[fmap[i]] = (uint8_t)j;
	}
	AssertH(j < 256, 1005);
}

#undef       SET_BH
#undef     CLEAR_BH
#undef     ISSET_BH
#undef      WORD_BH
#undef UNALIGNED_BH


/*---------------------------------------------*/
/*--- The main, O(N^2 log(N)) sorting       ---*/
/*--- algorithm.  Faster for "normal"       ---*/
/*--- non-repetitive blocks.                ---*/
/*---------------------------------------------*/

/*---------------------------------------------*/
static
NOINLINE
int mainGtU(
		uint32_t  i1,
		uint32_t  i2,
		uint8_t*  block,
		uint16_t* quadrant,
		uint32_t  nblock,
		int32_t*  budget)
{
	int32_t  k;
	uint8_t  c1, c2;
	uint16_t s1, s2;

/* Loop unrolling here is actually very useful
 * (generated code is much simpler),
 * code size increase is only 270 bytes (i386)
 * but speeds up compression 10% overall
 */

#if CONFIG_BZIP2_FEATURE_SPEED >= 1

#define TIMES_8(code) \
	code; code; code; code; \
	code; code; code; code;
#define TIMES_12(code) \
	code; code; code; code; \
	code; code; code; code; \
	code; code; code; code;

#else

#define TIMES_8(code) \
{ \
	int nn = 8; \
	do { \
		code; \
	} while (--nn); \
}
#define TIMES_12(code) \
{ \
	int nn = 12; \
	do { \
		code; \
	} while (--nn); \
}

#endif

	AssertD(i1 != i2, "mainGtU");
	TIMES_12(
		c1 = block[i1]; c2 = block[i2];
		if (c1 != c2) return (c1 > c2);
		i1++; i2++;
	)

	k = nblock + 8;

	do {
		TIMES_8(
			c1 = block[i1]; c2 = block[i2];
			if (c1 != c2) return (c1 > c2);
			s1 = quadrant[i1]; s2 = quadrant[i2];
			if (s1 != s2) return (s1 > s2);
			i1++; i2++;
		)

		if (i1 >= nblock) i1 -= nblock;
		if (i2 >= nblock) i2 -= nblock;

		(*budget)--;
		k -= 8;
	} while (k >= 0);

	return False;
}
#undef TIMES_8
#undef TIMES_12

/*---------------------------------------------*/
/*
 * Knuth's increments seem to work better
 * than Incerpi-Sedgewick here.  Possibly
 * because the number of elems to sort is
 * usually small, typically <= 20.
 */
static
const int32_t incs[14] = {
	1, 4, 13, 40, 121, 364, 1093, 3280,
	9841, 29524, 88573, 265720,
	797161, 2391484
};

static
void mainSimpleSort(uint32_t* ptr,
		uint8_t*  block,
		uint16_t* quadrant,
		int32_t   nblock,
		int32_t   lo,
		int32_t   hi,
		int32_t   d,
		int32_t*  budget)
{
	int32_t i, j, h, bigN, hp;
	uint32_t v;

	bigN = hi - lo + 1;
	if (bigN < 2) return;

	hp = 0;
	while (incs[hp] < bigN) hp++;
	hp--;

	for (; hp >= 0; hp--) {
		h = incs[hp];

		i = lo + h;
		while (1) {
			/*-- copy 1 --*/
			if (i > hi) break;
			v = ptr[i];
			j = i;
			while (mainGtU(ptr[j-h]+d, v+d, block, quadrant, nblock, budget)) {
				ptr[j] = ptr[j-h];
				j = j - h;
				if (j <= (lo + h - 1)) break;
			}
			ptr[j] = v;
			i++;

/* 1.5% overall speedup, +290 bytes */
#if CONFIG_BZIP2_FEATURE_SPEED >= 3
			/*-- copy 2 --*/
			if (i > hi) break;
			v = ptr[i];
			j = i;
			while (mainGtU(ptr[j-h]+d, v+d, block, quadrant, nblock, budget)) {
				ptr[j] = ptr[j-h];
				j = j - h;
				if (j <= (lo + h - 1)) break;
			}
			ptr[j] = v;
			i++;

			/*-- copy 3 --*/
			if (i > hi) break;
			v = ptr[i];
			j = i;
			while (mainGtU(ptr[j-h]+d, v+d, block, quadrant, nblock, budget)) {
				ptr[j] = ptr[j-h];
				j = j - h;
				if (j <= (lo + h - 1)) break;
			}
			ptr[j] = v;
			i++;
#endif
			if (*budget < 0) return;
		}
	}
}


/*---------------------------------------------*/
/*
 * The following is an implementation of
 * an elegant 3-way quicksort for strings,
 * described in a paper "Fast Algorithms for
 * Sorting and Searching Strings", by Robert
 * Sedgewick and Jon L. Bentley.
 */

static
ALWAYS_INLINE
uint8_t mmed3(uint8_t a, uint8_t b, uint8_t c)
{
	uint8_t t;
	if (a > b) {
		t = a;
		a = b;
		b = t;
	};
	/* here b >= a */
	if (b > c) {
		b = c;
		if (a > b)
			b = a;
	}
	return b;
}

#define mpush(lz,hz,dz) \
{ \
	stackLo[sp] = lz; \
	stackHi[sp] = hz; \
	stackD [sp] = dz; \
	sp++; \
}

#define mpop(lz,hz,dz) \
{ \
	sp--; \
	lz = stackLo[sp]; \
	hz = stackHi[sp]; \
	dz = stackD [sp]; \
}

#define mnextsize(az) (nextHi[az] - nextLo[az])

#define mnextswap(az,bz) \
{ \
	int32_t tz; \
	tz = nextLo[az]; nextLo[az] = nextLo[bz]; nextLo[bz] = tz; \
	tz = nextHi[az]; nextHi[az] = nextHi[bz]; nextHi[bz] = tz; \
	tz = nextD [az]; nextD [az] = nextD [bz]; nextD [bz] = tz; \
}

#define MAIN_QSORT_SMALL_THRESH 20
#define MAIN_QSORT_DEPTH_THRESH (BZ_N_RADIX + BZ_N_QSORT)
#define MAIN_QSORT_STACK_SIZE   100

static
void mainQSort3(uint32_t* ptr,
		uint8_t*  block,
		uint16_t* quadrant,
		int32_t   nblock,
		int32_t   loSt,
		int32_t   hiSt,
		int32_t   dSt,
		int32_t*  budget)
{
	int32_t unLo, unHi, ltLo, gtHi, n, m, med;
	int32_t sp, lo, hi, d;

	int32_t stackLo[MAIN_QSORT_STACK_SIZE];
	int32_t stackHi[MAIN_QSORT_STACK_SIZE];
	int32_t stackD [MAIN_QSORT_STACK_SIZE];

	int32_t nextLo[3];
	int32_t nextHi[3];
	int32_t nextD [3];

	sp = 0;
	mpush(loSt, hiSt, dSt);

	while (sp > 0) {
		AssertH(sp < MAIN_QSORT_STACK_SIZE - 2, 1001);

		mpop(lo, hi, d);
		if (hi - lo < MAIN_QSORT_SMALL_THRESH
		 || d > MAIN_QSORT_DEPTH_THRESH
		) {
			mainSimpleSort(ptr, block, quadrant, nblock, lo, hi, d, budget);
			if (*budget < 0)
				return;
			continue;
		}
		med = (int32_t)	mmed3(block[ptr[lo          ] + d],
		                      block[ptr[hi          ] + d],
		                      block[ptr[(lo+hi) >> 1] + d]);

		unLo = ltLo = lo;
		unHi = gtHi = hi;

		while (1) {
			while (1) {
				if (unLo > unHi)
					break;
				n = ((int32_t)block[ptr[unLo]+d]) - med;
				if (n == 0) {
					mswap(ptr[unLo], ptr[ltLo]);
					ltLo++;
					unLo++;
					continue;
				};
				if (n >  0) break;
				unLo++;
			}
			while (1) {
				if (unLo > unHi)
					break;
				n = ((int32_t)block[ptr[unHi]+d]) - med;
				if (n == 0) {
					mswap(ptr[unHi], ptr[gtHi]);
					gtHi--;
					unHi--;
					continue;
				};
				if (n <  0) break;
				unHi--;
			}
			if (unLo > unHi)
				break;
			mswap(ptr[unLo], ptr[unHi]);
			unLo++;
			unHi--;
		}

		AssertD(unHi == unLo-1, "mainQSort3(2)");

		if (gtHi < ltLo) {
			mpush(lo, hi, d + 1);
			continue;
		}

		n = mmin(ltLo-lo, unLo-ltLo); mvswap(ptr, lo, unLo-n, n);
		m = mmin(hi-gtHi, gtHi-unHi); mvswap(ptr, unLo, hi-m+1, m);

		n = lo + unLo - ltLo - 1;
		m = hi - (gtHi - unHi) + 1;

		nextLo[0] = lo;  nextHi[0] = n;   nextD[0] = d;
		nextLo[1] = m;   nextHi[1] = hi;  nextD[1] = d;
		nextLo[2] = n+1; nextHi[2] = m-1; nextD[2] = d+1;

		if (mnextsize(0) < mnextsize(1)) mnextswap(0, 1);
		if (mnextsize(1) < mnextsize(2)) mnextswap(1, 2);
		if (mnextsize(0) < mnextsize(1)) mnextswap(0, 1);

		AssertD (mnextsize(0) >= mnextsize(1), "mainQSort3(8)");
		AssertD (mnextsize(1) >= mnextsize(2), "mainQSort3(9)");

		mpush(nextLo[0], nextHi[0], nextD[0]);
		mpush(nextLo[1], nextHi[1], nextD[1]);
		mpush(nextLo[2], nextHi[2], nextD[2]);
	}
}

#undef mpush
#undef mpop
#undef mnextsize
#undef mnextswap
#undef MAIN_QSORT_SMALL_THRESH
#undef MAIN_QSORT_DEPTH_THRESH
#undef MAIN_QSORT_STACK_SIZE


/*---------------------------------------------*/
/* Pre:
 *	nblock > N_OVERSHOOT
 *	block32 exists for [0 .. nblock-1 +N_OVERSHOOT]
 *	((uint8_t*)block32) [0 .. nblock-1] holds block
 *	ptr exists for [0 .. nblock-1]
 *
 * Post:
 *	((uint8_t*)block32) [0 .. nblock-1] holds block
 *	All other areas of block32 destroyed
 *	ftab[0 .. 65536] destroyed
 *	ptr [0 .. nblock-1] holds sorted order
 *	if (*budget < 0), sorting was abandoned
 */

#define BIGFREQ(b) (ftab[((b)+1) << 8] - ftab[(b) << 8])
#define SETMASK (1 << 21)
#define CLEARMASK (~(SETMASK))

static NOINLINE
void mainSort(EState* state,
		uint32_t* ptr,
		uint8_t*  block,
		uint16_t* quadrant,
		uint32_t* ftab,
		int32_t   nblock,
		int32_t*  budget)
{
	int32_t  i, j, k, ss, sb;
	uint8_t  c1;
	int32_t  numQSorted;
	uint16_t s;
	Bool     bigDone[256];
	/* bbox: moved to EState to save stack
	int32_t  runningOrder[256];
	int32_t  copyStart[256];
	int32_t  copyEnd  [256];
	*/
#define runningOrder (state->mainSort__runningOrder)
#define copyStart    (state->mainSort__copyStart)
#define copyEnd      (state->mainSort__copyEnd)

	/*-- set up the 2-byte frequency table --*/
	/* was: for (i = 65536; i >= 0; i--) ftab[i] = 0; */
	memset(ftab, 0, 65537 * sizeof(ftab[0]));

	j = block[0] << 8;
	i = nblock - 1;
/* 3%, +300 bytes */
#if CONFIG_BZIP2_FEATURE_SPEED >= 2
	for (; i >= 3; i -= 4) {
		quadrant[i] = 0;
		j = (j >> 8) | (((uint16_t)block[i]) << 8);
		ftab[j]++;
		quadrant[i-1] = 0;
		j = (j >> 8) | (((uint16_t)block[i-1]) << 8);
		ftab[j]++;
		quadrant[i-2] = 0;
		j = (j >> 8) | (((uint16_t)block[i-2]) << 8);
		ftab[j]++;
		quadrant[i-3] = 0;
		j = (j >> 8) | (((uint16_t)block[i-3]) << 8);
		ftab[j]++;
	}
#endif
	for (; i >= 0; i--) {
		quadrant[i] = 0;
		j = (j >> 8) | (((uint16_t)block[i]) << 8);
		ftab[j]++;
	}

	/*-- (emphasises close relationship of block & quadrant) --*/
	for (i = 0; i < BZ_N_OVERSHOOT; i++) {
		block   [nblock+i] = block[i];
		quadrant[nblock+i] = 0;
	}

	/*-- Complete the initial radix sort --*/
	j = ftab[0]; /* bbox: optimized */
	for (i = 1; i <= 65536; i++) {
		j += ftab[i];
		ftab[i] = j;
	}

	s = block[0] << 8;
	i = nblock - 1;
#if CONFIG_BZIP2_FEATURE_SPEED >= 2
	for (; i >= 3; i -= 4) {
		s = (s >> 8) | (block[i] << 8);
		j = ftab[s] - 1;
		ftab[s] = j;
		ptr[j] = i;
		s = (s >> 8) | (block[i-1] << 8);
		j = ftab[s] - 1;
		ftab[s] = j;
		ptr[j] = i-1;
		s = (s >> 8) | (block[i-2] << 8);
		j = ftab[s] - 1;
		ftab[s] = j;
		ptr[j] = i-2;
		s = (s >> 8) | (block[i-3] << 8);
		j = ftab[s] - 1;
		ftab[s] = j;
		ptr[j] = i-3;
	}
#endif
	for (; i >= 0; i--) {
		s = (s >> 8) | (block[i] << 8);
		j = ftab[s] - 1;
		ftab[s] = j;
		ptr[j] = i;
	}

	/*
	 * Now ftab contains the first loc of every small bucket.
	 * Calculate the running order, from smallest to largest
	 * big bucket.
	 */
	for (i = 0; i <= 255; i++) {
		bigDone     [i] = False;
		runningOrder[i] = i;
	}

	{
		int32_t vv;
		/* bbox: was: int32_t h = 1; */
		/* do h = 3 * h + 1; while (h <= 256); */
		uint32_t h = 364;

		do {
			/*h = h / 3;*/
			h = (h * 171) >> 9; /* bbox: fast h/3 */
			for (i = h; i <= 255; i++) {
				vv = runningOrder[i];
				j = i;
				while (BIGFREQ(runningOrder[j-h]) > BIGFREQ(vv)) {
					runningOrder[j] = runningOrder[j-h];
					j = j - h;
					if (j <= (h - 1))
						goto zero;
				}
 zero:
				runningOrder[j] = vv;
			}
		} while (h != 1);
	}

	/*
	 * The main sorting loop.
	 */

	numQSorted = 0;

	for (i = 0; i <= 255; i++) {

		/*
		 * Process big buckets, starting with the least full.
		 * Basically this is a 3-step process in which we call
		 * mainQSort3 to sort the small buckets [ss, j], but
		 * also make a big effort to avoid the calls if we can.
		 */
		ss = runningOrder[i];

		/*
		 * Step 1:
		 * Complete the big bucket [ss] by quicksorting
		 * any unsorted small buckets [ss, j], for j != ss.
		 * Hopefully previous pointer-scanning phases have already
		 * completed many of the small buckets [ss, j], so
		 * we don't have to sort them at all.
		 */
		for (j = 0; j <= 255; j++) {
			if (j != ss) {
				sb = (ss << 8) + j;
				if (!(ftab[sb] & SETMASK)) {
					int32_t lo =  ftab[sb]   & CLEARMASK;
					int32_t hi = (ftab[sb+1] & CLEARMASK) - 1;
					if (hi > lo) {
						mainQSort3(
							ptr, block, quadrant, nblock,
							lo, hi, BZ_N_RADIX, budget
						);
						if (*budget < 0) return;
						numQSorted += (hi - lo + 1);
					}
				}
				ftab[sb] |= SETMASK;
			}
		}

		AssertH(!bigDone[ss], 1006);

		/*
		 * Step 2:
		 * Now scan this big bucket [ss] so as to synthesise the
		 * sorted order for small buckets [t, ss] for all t,
		 * including, magically, the bucket [ss,ss] too.
		 * This will avoid doing Real Work in subsequent Step 1's.
		 */
		{
			for (j = 0; j <= 255; j++) {
				copyStart[j] =  ftab[(j << 8) + ss]     & CLEARMASK;
				copyEnd  [j] = (ftab[(j << 8) + ss + 1] & CLEARMASK) - 1;
			}
			for (j = ftab[ss << 8] & CLEARMASK; j < copyStart[ss]; j++) {
				k = ptr[j] - 1;
				if (k < 0)
					k += nblock;
				c1 = block[k];
				if (!bigDone[c1])
					ptr[copyStart[c1]++] = k;
			}
			for (j = (ftab[(ss+1) << 8] & CLEARMASK) - 1; j > copyEnd[ss]; j--) {
				k = ptr[j]-1;
				if (k < 0)
					k += nblock;
				c1 = block[k];
				if (!bigDone[c1])
					ptr[copyEnd[c1]--] = k;
			}
		}

		/* Extremely rare case missing in bzip2-1.0.0 and 1.0.1.
		 * Necessity for this case is demonstrated by compressing
		 * a sequence of approximately 48.5 million of character
		 * 251; 1.0.0/1.0.1 will then die here. */
		AssertH((copyStart[ss]-1 == copyEnd[ss]) \
		     || (copyStart[ss] == 0 && copyEnd[ss] == nblock-1), 1007);

		for (j = 0; j <= 255; j++)
			ftab[(j << 8) + ss] |= SETMASK;

		/*
		 * Step 3:
		 * The [ss] big bucket is now done.  Record this fact,
		 * and update the quadrant descriptors.  Remember to
		 * update quadrants in the overshoot area too, if
		 * necessary.  The "if (i < 255)" test merely skips
		 * this updating for the last bucket processed, since
		 * updating for the last bucket is pointless.
		 *
		 * The quadrant array provides a way to incrementally
		 * cache sort orderings, as they appear, so as to
		 * make subsequent comparisons in fullGtU() complete
		 * faster.  For repetitive blocks this makes a big
		 * difference (but not big enough to be able to avoid
		 * the fallback sorting mechanism, exponential radix sort).
		 *
		 * The precise meaning is: at all times:
		 *
		 *	for 0 <= i < nblock and 0 <= j <= nblock
		 *
		 *	if block[i] != block[j],
		 *
		 *		then the relative values of quadrant[i] and
		 *			  quadrant[j] are meaningless.
		 *
		 *		else {
		 *			if quadrant[i] < quadrant[j]
		 *				then the string starting at i lexicographically
		 *				precedes the string starting at j
		 *
		 *			else if quadrant[i] > quadrant[j]
		 *				then the string starting at j lexicographically
		 *				precedes the string starting at i
		 *
		 *			else
		 *				the relative ordering of the strings starting
		 *				at i and j has not yet been determined.
		 *		}
		 */
		bigDone[ss] = True;

		if (i < 255) {
			int32_t bbStart = ftab[ss << 8] & CLEARMASK;
			int32_t bbSize  = (ftab[(ss+1) << 8] & CLEARMASK) - bbStart;
			int32_t shifts  = 0;

			while ((bbSize >> shifts) > 65534) shifts++;

			for (j = bbSize-1; j >= 0; j--) {
				int32_t a2update   = ptr[bbStart + j];
				uint16_t qVal      = (uint16_t)(j >> shifts);
				quadrant[a2update] = qVal;
				if (a2update < BZ_N_OVERSHOOT)
					quadrant[a2update + nblock] = qVal;
			}
			AssertH(((bbSize-1) >> shifts) <= 65535, 1002);
		}
	}
#undef runningOrder
#undef copyStart
#undef copyEnd
}

#undef BIGFREQ
#undef SETMASK
#undef CLEARMASK


/*---------------------------------------------*/
/* Pre:
 *	nblock > 0
 *	arr2 exists for [0 .. nblock-1 +N_OVERSHOOT]
 *	  ((uint8_t*)arr2)[0 .. nblock-1] holds block
 *	arr1 exists for [0 .. nblock-1]
 *
 * Post:
 *	((uint8_t*)arr2) [0 .. nblock-1] holds block
 *	All other areas of block destroyed
 *	ftab[0 .. 65536] destroyed
 *	arr1[0 .. nblock-1] holds sorted order
 */
static NOINLINE
void BZ2_blockSort(EState* s)
{
	/* In original bzip2 1.0.4, it's a parameter, but 30
	 * (which was the default) should work ok. */
	enum { wfact = 30 };

	uint32_t* ptr    = s->ptr;
	uint8_t*  block  = s->block;
	uint32_t* ftab   = s->ftab;
	int32_t   nblock = s->nblock;
	uint16_t* quadrant;
	int32_t   budget;
	int32_t   i;

	if (nblock < 10000) {
		fallbackSort(s->arr1, s->arr2, ftab, nblock);
	} else {
		/* Calculate the location for quadrant, remembering to get
		 * the alignment right.  Assumes that &(block[0]) is at least
		 * 2-byte aligned -- this should be ok since block is really
		 * the first section of arr2.
		 */
		i = nblock + BZ_N_OVERSHOOT;
		if (i & 1) i++;
		quadrant = (uint16_t*)(&(block[i]));

		/* (wfact-1) / 3 puts the default-factor-30
		 * transition point at very roughly the same place as
		 * with v0.1 and v0.9.0.
		 * Not that it particularly matters any more, since the
		 * resulting compressed stream is now the same regardless
		 * of whether or not we use the main sort or fallback sort.
		 */
		budget = nblock * ((wfact-1) / 3);

		mainSort(s, ptr, block, quadrant, ftab, nblock, &budget);
		if (budget < 0) {
			fallbackSort(s->arr1, s->arr2, ftab, nblock);
		}
	}

	s->origPtr = -1;
	for (i = 0; i < s->nblock; i++)
		if (ptr[i] == 0) {
			s->origPtr = i;
			break;
		};

	AssertH(s->origPtr != -1, 1003);
}


/*-------------------------------------------------------------*/
/*--- end                                       blocksort.c ---*/
/*-------------------------------------------------------------*/