Coverage Report

Created: 2025-10-13 06:36

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