bzip2: port bzip2 1.0.4 to busybox. note: bzip2 code resides
[people/mcb30/busybox.git] / archival / bz / blocksort.c
1 /*
2  * bzip2 is written by Julian Seward <jseward@bzip.org>.
3  * Adapted for busybox by Denys Vlasenko <vda.linux@googlemail.com>.
4  * See README and LICENSE files in this directory for more information.
5  */
6
7 /*-------------------------------------------------------------*/
8 /*--- Block sorting machinery                               ---*/
9 /*---                                           blocksort.c ---*/
10 /*-------------------------------------------------------------*/
11
12 /* ------------------------------------------------------------------
13 This file is part of bzip2/libbzip2, a program and library for
14 lossless, block-sorting data compression.
15
16 bzip2/libbzip2 version 1.0.4 of 20 December 2006
17 Copyright (C) 1996-2006 Julian Seward <jseward@bzip.org>
18
19 Please read the WARNING, DISCLAIMER and PATENTS sections in the
20 README file.
21
22 This program is released under the terms of the license contained
23 in the file LICENSE.
24 ------------------------------------------------------------------ */
25
26 /* #include "bzlib_private.h" */
27
28 /*---------------------------------------------*/
29 /*--- Fallback O(N log(N)^2) sorting        ---*/
30 /*--- algorithm, for repetitive blocks      ---*/
31 /*---------------------------------------------*/
32
33 /*---------------------------------------------*/
34 static
35 inline
36 void fallbackSimpleSort(uint32_t* fmap,
37                 uint32_t* eclass,
38                 int32_t   lo,
39                 int32_t   hi)
40 {
41         int32_t i, j, tmp;
42         uint32_t ec_tmp;
43
44         if (lo == hi) return;
45
46         if (hi - lo > 3) {
47                 for (i = hi-4; i >= lo; i--) {
48                         tmp = fmap[i];
49                         ec_tmp = eclass[tmp];
50                         for (j = i+4; j <= hi && ec_tmp > eclass[fmap[j]]; j += 4)
51                                 fmap[j-4] = fmap[j];
52                         fmap[j-4] = tmp;
53                 }
54         }
55
56         for (i = hi-1; i >= lo; i--) {
57                 tmp = fmap[i];
58                 ec_tmp = eclass[tmp];
59                 for (j = i+1; j <= hi && ec_tmp > eclass[fmap[j]]; j++)
60                         fmap[j-1] = fmap[j];
61                 fmap[j-1] = tmp;
62         }
63 }
64
65
66 /*---------------------------------------------*/
67 #define fswap(zz1, zz2) \
68 { \
69         int32_t zztmp = zz1; \
70         zz1 = zz2; \
71         zz2 = zztmp; \
72 }
73
74 #define fvswap(zzp1, zzp2, zzn) \
75 { \
76         int32_t yyp1 = (zzp1); \
77         int32_t yyp2 = (zzp2); \
78         int32_t yyn  = (zzn); \
79         while (yyn > 0) { \
80                 fswap(fmap[yyp1], fmap[yyp2]); \
81                 yyp1++; \
82                 yyp2++; \
83                 yyn--; \
84         } \
85 }
86
87
88 #define fmin(a,b) ((a) < (b)) ? (a) : (b)
89
90 #define fpush(lz,hz) { \
91         stackLo[sp] = lz; \
92         stackHi[sp] = hz; \
93         sp++; \
94 }
95
96 #define fpop(lz,hz) { \
97         sp--; \
98         lz = stackLo[sp]; \
99         hz = stackHi[sp]; \
100 }
101
102 #define FALLBACK_QSORT_SMALL_THRESH 10
103 #define FALLBACK_QSORT_STACK_SIZE   100
104
105
106 static
107 void fallbackQSort3(uint32_t* fmap,
108                 uint32_t* eclass,
109                 int32_t   loSt,
110                 int32_t   hiSt)
111 {
112         int32_t unLo, unHi, ltLo, gtHi, n, m;
113         int32_t sp, lo, hi;
114         uint32_t med, r, r3;
115         int32_t stackLo[FALLBACK_QSORT_STACK_SIZE];
116         int32_t stackHi[FALLBACK_QSORT_STACK_SIZE];
117
118         r = 0;
119
120         sp = 0;
121         fpush(loSt, hiSt);
122
123         while (sp > 0) {
124                 AssertH(sp < FALLBACK_QSORT_STACK_SIZE - 1, 1004);
125
126                 fpop(lo, hi);
127                 if (hi - lo < FALLBACK_QSORT_SMALL_THRESH) {
128                         fallbackSimpleSort(fmap, eclass, lo, hi);
129                         continue;
130                 }
131
132                 /* Random partitioning.  Median of 3 sometimes fails to
133                  * avoid bad cases.  Median of 9 seems to help but
134                  * looks rather expensive.  This too seems to work but
135                  * is cheaper.  Guidance for the magic constants
136                  * 7621 and 32768 is taken from Sedgewick's algorithms
137                  * book, chapter 35.
138                  */
139                 r = ((r * 7621) + 1) % 32768;
140                 r3 = r % 3;
141                 if (r3 == 0)
142                         med = eclass[fmap[lo]];
143                 else if (r3 == 1)
144                         med = eclass[fmap[(lo+hi)>>1]];
145                 else
146                         med = eclass[fmap[hi]];
147
148                 unLo = ltLo = lo;
149                 unHi = gtHi = hi;
150
151                 while (1) {
152                         while (1) {
153                                 if (unLo > unHi) break;
154                                 n = (int32_t)eclass[fmap[unLo]] - (int32_t)med;
155                                 if (n == 0) {
156                                         fswap(fmap[unLo], fmap[ltLo]);
157                                         ltLo++;
158                                         unLo++;
159                                         continue;
160                                 };
161                                 if (n > 0) break;
162                                 unLo++;
163                         }
164                         while (1) {
165                                 if (unLo > unHi) break;
166                                 n = (int32_t)eclass[fmap[unHi]] - (int32_t)med;
167                                 if (n == 0) {
168                                         fswap(fmap[unHi], fmap[gtHi]);
169                                         gtHi--; unHi--;
170                                         continue;
171                                 };
172                                 if (n < 0) break;
173                                 unHi--;
174                         }
175                         if (unLo > unHi) break;
176                         fswap(fmap[unLo], fmap[unHi]); unLo++; unHi--;
177                 }
178
179                 AssertD(unHi == unLo-1, "fallbackQSort3(2)");
180
181                 if (gtHi < ltLo) continue;
182
183                 n = fmin(ltLo-lo, unLo-ltLo); fvswap(lo, unLo-n, n);
184                 m = fmin(hi-gtHi, gtHi-unHi); fvswap(unLo, hi-m+1, m);
185
186                 n = lo + unLo - ltLo - 1;
187                 m = hi - (gtHi - unHi) + 1;
188
189                 if (n - lo > hi - m) {
190                         fpush(lo, n);
191                         fpush(m, hi);
192                 } else {
193                         fpush(m, hi);
194                         fpush(lo, n);
195                 }
196         }
197 }
198
199 #undef fmin
200 #undef fpush
201 #undef fpop
202 #undef fswap
203 #undef fvswap
204 #undef FALLBACK_QSORT_SMALL_THRESH
205 #undef FALLBACK_QSORT_STACK_SIZE
206
207
208 /*---------------------------------------------*/
209 /* Pre:
210  *      nblock > 0
211  *      eclass exists for [0 .. nblock-1]
212  *      ((UChar*)eclass) [0 .. nblock-1] holds block
213  *      ptr exists for [0 .. nblock-1]
214  *
215  * Post:
216  *      ((UChar*)eclass) [0 .. nblock-1] holds block
217  *      All other areas of eclass destroyed
218  *      fmap [0 .. nblock-1] holds sorted order
219  *      bhtab[0 .. 2+(nblock/32)] destroyed
220 */
221
222 #define       SET_BH(zz)  bhtab[(zz) >> 5] |= (1 << ((zz) & 31))
223 #define     CLEAR_BH(zz)  bhtab[(zz) >> 5] &= ~(1 << ((zz) & 31))
224 #define     ISSET_BH(zz)  (bhtab[(zz) >> 5] & (1 << ((zz) & 31)))
225 #define      WORD_BH(zz)  bhtab[(zz) >> 5]
226 #define UNALIGNED_BH(zz)  ((zz) & 0x01f)
227
228 static
229 void fallbackSort(uint32_t* fmap,
230                 uint32_t* eclass,
231                 uint32_t* bhtab,
232                 int32_t   nblock)
233 {
234         int32_t ftab[257];
235         int32_t ftabCopy[256];
236         int32_t H, i, j, k, l, r, cc, cc1;
237         int32_t nNotDone;
238         int32_t nBhtab;
239         UChar* eclass8 = (UChar*)eclass;
240
241         /*
242          * Initial 1-char radix sort to generate
243          * initial fmap and initial BH bits.
244          */
245         for (i = 0; i < 257;    i++) ftab[i] = 0;
246         for (i = 0; i < nblock; i++) ftab[eclass8[i]]++;
247         for (i = 0; i < 256;    i++) ftabCopy[i] = ftab[i];
248         for (i = 1; i < 257;    i++) ftab[i] += ftab[i-1];
249
250         for (i = 0; i < nblock; i++) {
251                 j = eclass8[i];
252                 k = ftab[j] - 1;
253                 ftab[j] = k;
254                 fmap[k] = i;
255         }
256
257         nBhtab = 2 + (nblock / 32);
258         for (i = 0; i < nBhtab; i++) bhtab[i] = 0;
259         for (i = 0; i < 256; i++) SET_BH(ftab[i]);
260
261         /*
262          * Inductively refine the buckets.  Kind-of an
263          * "exponential radix sort" (!), inspired by the
264          * Manber-Myers suffix array construction algorithm.
265          */
266
267         /*-- set sentinel bits for block-end detection --*/
268         for (i = 0; i < 32; i++) {
269                 SET_BH(nblock + 2*i);
270                 CLEAR_BH(nblock + 2*i + 1);
271         }
272
273         /*-- the log(N) loop --*/
274         H = 1;
275         while (1) {
276                 j = 0;
277                 for (i = 0; i < nblock; i++) {
278                         if (ISSET_BH(i))
279                                 j = i;
280                         k = fmap[i] - H;
281                         if (k < 0)
282                                 k += nblock;
283                         eclass[k] = j;
284                 }
285
286                 nNotDone = 0;
287                 r = -1;
288                 while (1) {
289
290          /*-- find the next non-singleton bucket --*/
291                         k = r + 1;
292                         while (ISSET_BH(k) && UNALIGNED_BH(k))
293                                 k++;
294                         if (ISSET_BH(k)) {
295                                 while (WORD_BH(k) == 0xffffffff) k += 32;
296                                 while (ISSET_BH(k)) k++;
297                         }
298                         l = k - 1;
299                         if (l >= nblock)
300                                 break;
301                         while (!ISSET_BH(k) && UNALIGNED_BH(k))
302                                 k++;
303                         if (!ISSET_BH(k)) {
304                                 while (WORD_BH(k) == 0x00000000) k += 32;
305                                 while (!ISSET_BH(k)) k++;
306                         }
307                         r = k - 1;
308                         if (r >= nblock)
309                                 break;
310
311                         /*-- now [l, r] bracket current bucket --*/
312                         if (r > l) {
313                                 nNotDone += (r - l + 1);
314                                 fallbackQSort3(fmap, eclass, l, r);
315
316                                 /*-- scan bucket and generate header bits-- */
317                                 cc = -1;
318                                 for (i = l; i <= r; i++) {
319                                         cc1 = eclass[fmap[i]];
320                                         if (cc != cc1) {
321                                                 SET_BH(i);
322                                                 cc = cc1;
323                                         };
324                                 }
325                         }
326                 }
327
328                 H *= 2;
329                 if (H > nblock || nNotDone == 0)
330                         break;
331         }
332
333         /*
334          * Reconstruct the original block in
335          * eclass8 [0 .. nblock-1], since the
336          * previous phase destroyed it.
337          */
338         j = 0;
339         for (i = 0; i < nblock; i++) {
340                 while (ftabCopy[j] == 0)
341                         j++;
342                 ftabCopy[j]--;
343                 eclass8[fmap[i]] = (UChar)j;
344         }
345         AssertH(j < 256, 1005);
346 }
347
348 #undef       SET_BH
349 #undef     CLEAR_BH
350 #undef     ISSET_BH
351 #undef      WORD_BH
352 #undef UNALIGNED_BH
353
354
355 /*---------------------------------------------*/
356 /*--- The main, O(N^2 log(N)) sorting       ---*/
357 /*--- algorithm.  Faster for "normal"       ---*/
358 /*--- non-repetitive blocks.                ---*/
359 /*---------------------------------------------*/
360
361 /*---------------------------------------------*/
362 static
363 inline
364 Bool mainGtU(
365                 uint32_t  i1,
366                 uint32_t  i2,
367                 UChar*    block,
368                 uint16_t* quadrant,
369                 uint32_t  nblock,
370                 int32_t*  budget)
371 {
372         int32_t  k;
373         UChar  c1, c2;
374         uint16_t s1, s2;
375
376 ///unrolling
377         AssertD(i1 != i2, "mainGtU");
378         /* 1 */
379         c1 = block[i1]; c2 = block[i2];
380         if (c1 != c2) return (c1 > c2);
381         i1++; i2++;
382         /* 2 */
383         c1 = block[i1]; c2 = block[i2];
384         if (c1 != c2) return (c1 > c2);
385         i1++; i2++;
386         /* 3 */
387         c1 = block[i1]; c2 = block[i2];
388         if (c1 != c2) return (c1 > c2);
389         i1++; i2++;
390         /* 4 */
391         c1 = block[i1]; c2 = block[i2];
392         if (c1 != c2) return (c1 > c2);
393         i1++; i2++;
394         /* 5 */
395         c1 = block[i1]; c2 = block[i2];
396         if (c1 != c2) return (c1 > c2);
397         i1++; i2++;
398         /* 6 */
399         c1 = block[i1]; c2 = block[i2];
400         if (c1 != c2) return (c1 > c2);
401         i1++; i2++;
402         /* 7 */
403         c1 = block[i1]; c2 = block[i2];
404         if (c1 != c2) return (c1 > c2);
405         i1++; i2++;
406         /* 8 */
407         c1 = block[i1]; c2 = block[i2];
408         if (c1 != c2) return (c1 > c2);
409         i1++; i2++;
410         /* 9 */
411         c1 = block[i1]; c2 = block[i2];
412         if (c1 != c2) return (c1 > c2);
413         i1++; i2++;
414         /* 10 */
415         c1 = block[i1]; c2 = block[i2];
416         if (c1 != c2) return (c1 > c2);
417         i1++; i2++;
418         /* 11 */
419         c1 = block[i1]; c2 = block[i2];
420         if (c1 != c2) return (c1 > c2);
421         i1++; i2++;
422         /* 12 */
423         c1 = block[i1]; c2 = block[i2];
424         if (c1 != c2) return (c1 > c2);
425         i1++; i2++;
426
427         k = nblock + 8;
428
429 ///unrolling
430         do {
431                 /* 1 */
432                 c1 = block[i1]; c2 = block[i2];
433                 if (c1 != c2) return (c1 > c2);
434                 s1 = quadrant[i1]; s2 = quadrant[i2];
435                 if (s1 != s2) return (s1 > s2);
436                 i1++; i2++;
437                 /* 2 */
438                 c1 = block[i1]; c2 = block[i2];
439                 if (c1 != c2) return (c1 > c2);
440                 s1 = quadrant[i1]; s2 = quadrant[i2];
441                 if (s1 != s2) return (s1 > s2);
442                 i1++; i2++;
443                 /* 3 */
444                 c1 = block[i1]; c2 = block[i2];
445                 if (c1 != c2) return (c1 > c2);
446                 s1 = quadrant[i1]; s2 = quadrant[i2];
447                 if (s1 != s2) return (s1 > s2);
448                 i1++; i2++;
449                 /* 4 */
450                 c1 = block[i1]; c2 = block[i2];
451                 if (c1 != c2) return (c1 > c2);
452                 s1 = quadrant[i1]; s2 = quadrant[i2];
453                 if (s1 != s2) return (s1 > s2);
454                 i1++; i2++;
455                 /* 5 */
456                 c1 = block[i1]; c2 = block[i2];
457                 if (c1 != c2) return (c1 > c2);
458                 s1 = quadrant[i1]; s2 = quadrant[i2];
459                 if (s1 != s2) return (s1 > s2);
460                 i1++; i2++;
461                 /* 6 */
462                 c1 = block[i1]; c2 = block[i2];
463                 if (c1 != c2) return (c1 > c2);
464                 s1 = quadrant[i1]; s2 = quadrant[i2];
465                 if (s1 != s2) return (s1 > s2);
466                 i1++; i2++;
467                 /* 7 */
468                 c1 = block[i1]; c2 = block[i2];
469                 if (c1 != c2) return (c1 > c2);
470                 s1 = quadrant[i1]; s2 = quadrant[i2];
471                 if (s1 != s2) return (s1 > s2);
472                 i1++; i2++;
473                 /* 8 */
474                 c1 = block[i1]; c2 = block[i2];
475                 if (c1 != c2) return (c1 > c2);
476                 s1 = quadrant[i1]; s2 = quadrant[i2];
477                 if (s1 != s2) return (s1 > s2);
478                 i1++; i2++;
479
480                 if (i1 >= nblock) i1 -= nblock;
481                 if (i2 >= nblock) i2 -= nblock;
482
483                 k -= 8;
484                 (*budget)--;
485         } while (k >= 0);
486
487         return False;
488 }
489
490
491 /*---------------------------------------------*/
492 /*
493  * Knuth's increments seem to work better
494  * than Incerpi-Sedgewick here.  Possibly
495  * because the number of elems to sort is
496  * usually small, typically <= 20.
497  */
498 static
499 const int32_t incs[14] = {
500         1, 4, 13, 40, 121, 364, 1093, 3280,
501         9841, 29524, 88573, 265720,
502         797161, 2391484
503 };
504
505 static
506 void mainSimpleSort(uint32_t* ptr,
507                 UChar*  block,
508                 uint16_t* quadrant,
509                 int32_t   nblock,
510                 int32_t   lo,
511                 int32_t   hi,
512                 int32_t   d,
513                 int32_t*  budget)
514 {
515         int32_t i, j, h, bigN, hp;
516         uint32_t v;
517
518         bigN = hi - lo + 1;
519         if (bigN < 2) return;
520
521         hp = 0;
522         while (incs[hp] < bigN) hp++;
523         hp--;
524
525         for (; hp >= 0; hp--) {
526                 h = incs[hp];
527
528                 i = lo + h;
529                 while (1) {
530
531 ///unrolling
532                         /*-- copy 1 --*/
533                         if (i > hi) break;
534                         v = ptr[i];
535                         j = i;
536                         while (mainGtU(ptr[j-h]+d, v+d, block, quadrant, nblock, budget)) {
537                                 ptr[j] = ptr[j-h];
538                                 j = j - h;
539                                 if (j <= (lo + h - 1)) break;
540                         }
541                         ptr[j] = v;
542                         i++;
543
544                         /*-- copy 2 --*/
545                         if (i > hi) break;
546                         v = ptr[i];
547                         j = i;
548                         while (mainGtU(ptr[j-h]+d, v+d, block, quadrant, nblock, budget)) {
549                                 ptr[j] = ptr[j-h];
550                                 j = j - h;
551                                 if (j <= (lo + h - 1)) break;
552                         }
553                         ptr[j] = v;
554                         i++;
555
556                         /*-- copy 3 --*/
557                         if (i > hi) break;
558                         v = ptr[i];
559                         j = i;
560                         while (mainGtU (ptr[j-h]+d, v+d, block, quadrant, nblock, budget)) {
561                                 ptr[j] = ptr[j-h];
562                                 j = j - h;
563                                 if (j <= (lo + h - 1)) break;
564                         }
565                         ptr[j] = v;
566                         i++;
567
568                         if (*budget < 0) return;
569                 }
570         }
571 }
572
573
574 /*---------------------------------------------*/
575 /*
576  * The following is an implementation of
577  * an elegant 3-way quicksort for strings,
578  * described in a paper "Fast Algorithms for
579  * Sorting and Searching Strings", by Robert
580  * Sedgewick and Jon L. Bentley.
581  */
582
583 #define mswap(zz1, zz2) \
584 { \
585         int32_t zztmp = zz1; \
586         zz1 = zz2; \
587         zz2 = zztmp; \
588 }
589
590 #define mvswap(zzp1, zzp2, zzn) \
591 { \
592         int32_t yyp1 = (zzp1); \
593         int32_t yyp2 = (zzp2); \
594         int32_t yyn  = (zzn); \
595         while (yyn > 0) { \
596                 mswap(ptr[yyp1], ptr[yyp2]); \
597                 yyp1++; \
598                 yyp2++; \
599                 yyn--; \
600         } \
601 }
602
603 static
604 inline
605 UChar mmed3(UChar a, UChar b, UChar c)
606 {
607         UChar t;
608         if (a > b) {
609                 t = a;
610                 a = b;
611                 b = t;
612         };
613         if (b > c) {
614                 b = c;
615                 if (a > b)
616                         b = a;
617         }
618         return b;
619 }
620
621 #define mmin(a,b) ((a) < (b)) ? (a) : (b)
622
623 #define mpush(lz,hz,dz) \
624 { \
625         stackLo[sp] = lz; \
626         stackHi[sp] = hz; \
627         stackD [sp] = dz; \
628         sp++; \
629 }
630
631 #define mpop(lz,hz,dz) \
632 { \
633         sp--; \
634         lz = stackLo[sp]; \
635         hz = stackHi[sp]; \
636         dz = stackD [sp]; \
637 }
638
639
640 #define mnextsize(az) (nextHi[az]-nextLo[az])
641
642 #define mnextswap(az,bz) \
643 { \
644         int32_t tz; \
645         tz = nextLo[az]; nextLo[az] = nextLo[bz]; nextLo[bz] = tz; \
646         tz = nextHi[az]; nextHi[az] = nextHi[bz]; nextHi[bz] = tz; \
647         tz = nextD [az]; nextD [az] = nextD [bz]; nextD [bz] = tz; \
648 }
649
650 #define MAIN_QSORT_SMALL_THRESH 20
651 #define MAIN_QSORT_DEPTH_THRESH (BZ_N_RADIX + BZ_N_QSORT)
652 #define MAIN_QSORT_STACK_SIZE   100
653
654 static
655 void mainQSort3(uint32_t* ptr,
656                 UChar*    block,
657                 uint16_t* quadrant,
658                 int32_t   nblock,
659                 int32_t   loSt,
660                 int32_t   hiSt,
661                 int32_t   dSt,
662                 int32_t*  budget)
663 {
664         int32_t unLo, unHi, ltLo, gtHi, n, m, med;
665         int32_t sp, lo, hi, d;
666
667         int32_t stackLo[MAIN_QSORT_STACK_SIZE];
668         int32_t stackHi[MAIN_QSORT_STACK_SIZE];
669         int32_t stackD [MAIN_QSORT_STACK_SIZE];
670
671         int32_t nextLo[3];
672         int32_t nextHi[3];
673         int32_t nextD [3];
674
675         sp = 0;
676         mpush(loSt, hiSt, dSt);
677
678         while (sp > 0) {
679                 AssertH(sp < MAIN_QSORT_STACK_SIZE - 2, 1001);
680
681                 mpop(lo, hi, d);
682                 if (hi - lo < MAIN_QSORT_SMALL_THRESH
683                  || d > MAIN_QSORT_DEPTH_THRESH
684                 ) {
685                         mainSimpleSort(ptr, block, quadrant, nblock, lo, hi, d, budget);
686                         if (*budget < 0)
687                                 return;
688                         continue;
689                 }
690
691                 med = (int32_t) mmed3(block[ptr[lo          ] + d],
692                                       block[ptr[hi          ] + d],
693                                       block[ptr[(lo+hi) >> 1] + d]);
694
695                 unLo = ltLo = lo;
696                 unHi = gtHi = hi;
697
698                 while (1) {
699                         while (1) {
700                                 if (unLo > unHi)
701                                         break;
702                                 n = ((int32_t)block[ptr[unLo]+d]) - med;
703                                 if (n == 0) {
704                                         mswap(ptr[unLo], ptr[ltLo]);
705                                         ltLo++;
706                                         unLo++;
707                                         continue;
708                                 };
709                                 if (n >  0) break;
710                                 unLo++;
711                         }
712                         while (1) {
713                                 if (unLo > unHi)
714                                         break;
715                                 n = ((int32_t)block[ptr[unHi]+d]) - med;
716                                 if (n == 0) {
717                                         mswap(ptr[unHi], ptr[gtHi]);
718                                         gtHi--;
719                                         unHi--;
720                                         continue;
721                                 };
722                                 if (n <  0) break;
723                                 unHi--;
724                         }
725                         if (unLo > unHi)
726                                 break;
727                         mswap(ptr[unLo], ptr[unHi]);
728                         unLo++;
729                         unHi--;
730                 }
731
732                 AssertD(unHi == unLo-1, "mainQSort3(2)");
733
734                 if (gtHi < ltLo) {
735                         mpush(lo, hi, d + 1);
736                         continue;
737                 }
738
739                 n = mmin(ltLo-lo, unLo-ltLo); mvswap(lo, unLo-n, n);
740                 m = mmin(hi-gtHi, gtHi-unHi); mvswap(unLo, hi-m+1, m);
741
742                 n = lo + unLo - ltLo - 1;
743                 m = hi - (gtHi - unHi) + 1;
744
745                 nextLo[0] = lo;  nextHi[0] = n;   nextD[0] = d;
746                 nextLo[1] = m;   nextHi[1] = hi;  nextD[1] = d;
747                 nextLo[2] = n+1; nextHi[2] = m-1; nextD[2] = d+1;
748
749                 if (mnextsize(0) < mnextsize(1)) mnextswap(0,1);
750                 if (mnextsize(1) < mnextsize(2)) mnextswap(1,2);
751                 if (mnextsize(0) < mnextsize(1)) mnextswap(0,1);
752
753                 AssertD (mnextsize(0) >= mnextsize(1), "mainQSort3(8)");
754                 AssertD (mnextsize(1) >= mnextsize(2), "mainQSort3(9)");
755
756                 mpush (nextLo[0], nextHi[0], nextD[0]);
757                 mpush (nextLo[1], nextHi[1], nextD[1]);
758                 mpush (nextLo[2], nextHi[2], nextD[2]);
759         }
760 }
761
762 #undef mswap
763 #undef mvswap
764 #undef mpush
765 #undef mpop
766 #undef mmin
767 #undef mnextsize
768 #undef mnextswap
769 #undef MAIN_QSORT_SMALL_THRESH
770 #undef MAIN_QSORT_DEPTH_THRESH
771 #undef MAIN_QSORT_STACK_SIZE
772
773
774 /*---------------------------------------------*/
775 /* Pre:
776  *      nblock > N_OVERSHOOT
777  *      block32 exists for [0 .. nblock-1 +N_OVERSHOOT]
778  *      ((UChar*)block32) [0 .. nblock-1] holds block
779  *      ptr exists for [0 .. nblock-1]
780  *
781  * Post:
782  *      ((UChar*)block32) [0 .. nblock-1] holds block
783  *      All other areas of block32 destroyed
784  *      ftab[0 .. 65536] destroyed
785  *      ptr [0 .. nblock-1] holds sorted order
786  *      if (*budget < 0), sorting was abandoned
787  */
788
789 #define BIGFREQ(b) (ftab[((b)+1) << 8] - ftab[(b) << 8])
790 #define SETMASK (1 << 21)
791 #define CLEARMASK (~(SETMASK))
792
793 static NOINLINE
794 void mainSort(uint32_t*   ptr,
795                 UChar*    block,
796                 uint16_t* quadrant,
797                 uint32_t* ftab,
798                 int32_t   nblock,
799                 int32_t*  budget)
800 {
801         int32_t  i, j, k, ss, sb;
802         int32_t  runningOrder[256];
803         Bool   bigDone[256];
804         int32_t  copyStart[256];
805         int32_t  copyEnd  [256];
806         UChar  c1;
807         int32_t  numQSorted;
808         uint16_t s;
809
810         /*-- set up the 2-byte frequency table --*/
811         /* was: for (i = 65536; i >= 0; i--) ftab[i] = 0; */
812         memset(ftab, 0, 65537 * sizeof(ftab[0]));
813
814         j = block[0] << 8;
815         i = nblock-1;
816 #if 0
817         for (; i >= 3; i -= 4) {
818                 quadrant[i] = 0;
819                 j = (j >> 8) |(((uint16_t)block[i]) << 8);
820                 ftab[j]++;
821                 quadrant[i-1] = 0;
822                 j = (j >> 8) |(((uint16_t)block[i-1]) << 8);
823                 ftab[j]++;
824                 quadrant[i-2] = 0;
825                 j = (j >> 8) |(((uint16_t)block[i-2]) << 8);
826                 ftab[j]++;
827                 quadrant[i-3] = 0;
828                 j = (j >> 8) |(((uint16_t)block[i-3]) << 8);
829                 ftab[j]++;
830         }
831 #endif
832         for (; i >= 0; i--) {
833                 quadrant[i] = 0;
834                 j = (j >> 8) |(((uint16_t)block[i]) << 8);
835                 ftab[j]++;
836         }
837
838         /*-- (emphasises close relationship of block & quadrant) --*/
839         for (i = 0; i < BZ_N_OVERSHOOT; i++) {
840                 block   [nblock+i] = block[i];
841                 quadrant[nblock+i] = 0;
842         }
843
844         /*-- Complete the initial radix sort --*/
845         for (i = 1; i <= 65536; i++) ftab[i] += ftab[i-1];
846
847         s = block[0] << 8;
848         i = nblock-1;
849 #if 0
850         for (; i >= 3; i -= 4) {
851                 s = (s >> 8) | (block[i] << 8);
852                 j = ftab[s] -1;
853                 ftab[s] = j;
854                 ptr[j] = i;
855                 s = (s >> 8) | (block[i-1] << 8);
856                 j = ftab[s] -1;
857                 ftab[s] = j;
858                 ptr[j] = i-1;
859                 s = (s >> 8) | (block[i-2] << 8);
860                 j = ftab[s] -1;
861                 ftab[s] = j;
862                 ptr[j] = i-2;
863                 s = (s >> 8) | (block[i-3] << 8);
864                 j = ftab[s] -1;
865                 ftab[s] = j;
866                 ptr[j] = i-3;
867         }
868 #endif
869         for (; i >= 0; i--) {
870                 s = (s >> 8) | (block[i] << 8);
871                 j = ftab[s] -1;
872                 ftab[s] = j;
873                 ptr[j] = i;
874         }
875
876         /*
877          * Now ftab contains the first loc of every small bucket.
878          * Calculate the running order, from smallest to largest
879          * big bucket.
880          */
881         for (i = 0; i <= 255; i++) {
882                 bigDone     [i] = False;
883                 runningOrder[i] = i;
884         }
885
886         {
887                 int32_t vv;
888                 /* was: int32_t h = 1; */
889                 /* do h = 3 * h + 1; while (h <= 256); */
890                 int32_t h = 364;
891
892                 do {
893                         h = h / 3;
894                         for (i = h; i <= 255; i++) {
895                                 vv = runningOrder[i];
896                                 j = i;
897                                 while (BIGFREQ(runningOrder[j-h]) > BIGFREQ(vv)) {
898                                         runningOrder[j] = runningOrder[j-h];
899                                         j = j - h;
900                                         if (j <= (h - 1)) goto zero;
901                                 }
902                                 zero:
903                                 runningOrder[j] = vv;
904                         }
905                 } while (h != 1);
906         }
907
908         /*
909          * The main sorting loop.
910          */
911
912         numQSorted = 0;
913
914         for (i = 0; i <= 255; i++) {
915
916                 /*
917                  * Process big buckets, starting with the least full.
918                  * Basically this is a 3-step process in which we call
919                  * mainQSort3 to sort the small buckets [ss, j], but
920                  * also make a big effort to avoid the calls if we can.
921                  */
922                 ss = runningOrder[i];
923
924                 /*
925                  * Step 1:
926                  * Complete the big bucket [ss] by quicksorting
927                  * any unsorted small buckets [ss, j], for j != ss.
928                  * Hopefully previous pointer-scanning phases have already
929                  * completed many of the small buckets [ss, j], so
930                  * we don't have to sort them at all.
931                  */
932                 for (j = 0; j <= 255; j++) {
933                         if (j != ss) {
934                                 sb = (ss << 8) + j;
935                                 if (!(ftab[sb] & SETMASK)) {
936                                         int32_t lo = ftab[sb]   & CLEARMASK;
937                                         int32_t hi = (ftab[sb+1] & CLEARMASK) - 1;
938                                         if (hi > lo) {
939                                                 mainQSort3 (
940                                                         ptr, block, quadrant, nblock,
941                                                         lo, hi, BZ_N_RADIX, budget
942                                                 );
943                                                 numQSorted += (hi - lo + 1);
944                                                 if (*budget < 0) return;
945                                         }
946                                 }
947                                 ftab[sb] |= SETMASK;
948                         }
949                 }
950
951                 AssertH(!bigDone[ss], 1006);
952
953                 /*
954                  * Step 2:
955                  * Now scan this big bucket [ss] so as to synthesise the
956                  * sorted order for small buckets [t, ss] for all t,
957                  * including, magically, the bucket [ss,ss] too.
958                  * This will avoid doing Real Work in subsequent Step 1's.
959                  */
960                 {
961                         for (j = 0; j <= 255; j++) {
962                                 copyStart[j] =  ftab[(j << 8) + ss]     & CLEARMASK;
963                                 copyEnd  [j] = (ftab[(j << 8) + ss + 1] & CLEARMASK) - 1;
964                         }
965                         for (j = ftab[ss << 8] & CLEARMASK; j < copyStart[ss]; j++) {
966                                 k = ptr[j] - 1;
967                                 if (k < 0)
968                                         k += nblock;
969                                 c1 = block[k];
970                                 if (!bigDone[c1])
971                                         ptr[copyStart[c1]++] = k;
972                         }
973                         for (j = (ftab[(ss+1) << 8] & CLEARMASK) - 1; j > copyEnd[ss]; j--) {
974                                 k = ptr[j]-1;
975                                 if (k < 0)
976                                         k += nblock;
977                                 c1 = block[k];
978                                 if (!bigDone[c1])
979                                         ptr[copyEnd[c1]--] = k;
980                         }
981                 }
982
983                 AssertH((copyStart[ss]-1 == copyEnd[ss])
984                                          ||
985                                         /* Extremely rare case missing in bzip2-1.0.0 and 1.0.1.
986                                          * Necessity for this case is demonstrated by compressing
987                                          * a sequence of approximately 48.5 million of character
988                                          * 251; 1.0.0/1.0.1 will then die here. */
989                                         (copyStart[ss] == 0 && copyEnd[ss] == nblock-1),
990                                         1007)
991
992                 for (j = 0; j <= 255; j++)
993                         ftab[(j << 8) + ss] |= SETMASK;
994
995                 /*
996                  * Step 3:
997                  * The [ss] big bucket is now done.  Record this fact,
998                  * and update the quadrant descriptors.  Remember to
999                  * update quadrants in the overshoot area too, if
1000                  * necessary.  The "if (i < 255)" test merely skips
1001                  * this updating for the last bucket processed, since
1002                  * updating for the last bucket is pointless.
1003                  *
1004                  * The quadrant array provides a way to incrementally
1005                  * cache sort orderings, as they appear, so as to
1006                  * make subsequent comparisons in fullGtU() complete
1007                  * faster.  For repetitive blocks this makes a big
1008                  * difference (but not big enough to be able to avoid
1009                  * the fallback sorting mechanism, exponential radix sort).
1010                  *
1011                  * The precise meaning is: at all times:
1012                  *
1013                  *      for 0 <= i < nblock and 0 <= j <= nblock
1014                  *
1015                  *      if block[i] != block[j],
1016                  *
1017                  *              then the relative values of quadrant[i] and
1018                  *                        quadrant[j] are meaningless.
1019                  *
1020                  *              else {
1021                  *                      if quadrant[i] < quadrant[j]
1022                  *                              then the string starting at i lexicographically
1023                  *                              precedes the string starting at j
1024                  *
1025                  *                      else if quadrant[i] > quadrant[j]
1026                  *                              then the string starting at j lexicographically
1027                  *                              precedes the string starting at i
1028                  *
1029                  *                      else
1030                  *                              the relative ordering of the strings starting
1031                  *                              at i and j has not yet been determined.
1032                  *              }
1033                  */
1034                 bigDone[ss] = True;
1035
1036                 if (i < 255) {
1037                         int32_t bbStart = ftab[ss << 8] & CLEARMASK;
1038                         int32_t bbSize  = (ftab[(ss+1) << 8] & CLEARMASK) - bbStart;
1039                         int32_t shifts  = 0;
1040
1041                         while ((bbSize >> shifts) > 65534) shifts++;
1042
1043                         for (j = bbSize-1; j >= 0; j--) {
1044                                 int32_t a2update     = ptr[bbStart + j];
1045                                 uint16_t qVal        = (uint16_t)(j >> shifts);
1046                                 quadrant[a2update] = qVal;
1047                                 if (a2update < BZ_N_OVERSHOOT)
1048                                         quadrant[a2update + nblock] = qVal;
1049                         }
1050                         AssertH(((bbSize-1) >> shifts) <= 65535, 1002);
1051                 }
1052
1053         }
1054 }
1055
1056 #undef BIGFREQ
1057 #undef SETMASK
1058 #undef CLEARMASK
1059
1060
1061 /*---------------------------------------------*/
1062 /* Pre:
1063  *      nblock > 0
1064  *      arr2 exists for [0 .. nblock-1 +N_OVERSHOOT]
1065  *        ((UChar*)arr2)[0 .. nblock-1] holds block
1066  *      arr1 exists for [0 .. nblock-1]
1067  *
1068  * Post:
1069  *      ((UChar*)arr2) [0 .. nblock-1] holds block
1070  *      All other areas of block destroyed
1071  *      ftab[0 .. 65536] destroyed
1072  *      arr1[0 .. nblock-1] holds sorted order
1073  */
1074 static NOINLINE
1075 void BZ2_blockSort(EState* s)
1076 {
1077         /* In original bzip2 1.0.4, it's a parameter, but 30
1078          * should work ok. */
1079         enum { wfact = 30 };
1080
1081         uint32_t* ptr    = s->ptr;
1082         UChar*    block  = s->block;
1083         uint32_t* ftab   = s->ftab;
1084         int32_t   nblock = s->nblock;
1085         uint16_t* quadrant;
1086         int32_t   budget;
1087         int32_t   i;
1088
1089         if (nblock < 10000) {
1090                 fallbackSort(s->arr1, s->arr2, ftab, nblock);
1091         } else {
1092                 /* Calculate the location for quadrant, remembering to get
1093                  * the alignment right.  Assumes that &(block[0]) is at least
1094                  * 2-byte aligned -- this should be ok since block is really
1095                  * the first section of arr2.
1096                  */
1097                 i = nblock + BZ_N_OVERSHOOT;
1098                 if (i & 1) i++;
1099                 quadrant = (uint16_t*)(&(block[i]));
1100
1101                 /* (wfact-1) / 3 puts the default-factor-30
1102                  * transition point at very roughly the same place as
1103                  * with v0.1 and v0.9.0.
1104                  * Not that it particularly matters any more, since the
1105                  * resulting compressed stream is now the same regardless
1106                  * of whether or not we use the main sort or fallback sort.
1107                  */
1108                 budget = nblock * ((wfact-1) / 3);
1109
1110                 mainSort(ptr, block, quadrant, ftab, nblock, &budget);
1111                 if (budget < 0) {
1112                         fallbackSort(s->arr1, s->arr2, ftab, nblock);
1113                 }
1114         }
1115
1116         s->origPtr = -1;
1117         for (i = 0; i < s->nblock; i++)
1118                 if (ptr[i] == 0) {
1119                         s->origPtr = i; break;
1120                 };
1121
1122         AssertH(s->origPtr != -1, 1003);
1123 }
1124
1125
1126 /*-------------------------------------------------------------*/
1127 /*--- end                                       blocksort.c ---*/
1128 /*-------------------------------------------------------------*/