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