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