contrib/bzip2/blocksort.c

   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.6 of 6 September 2010
  12    Copyright (C) 1996-2010 Julian Seward <jseward@bzip.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 {
  37    Int32 i, j, tmp;
  38    UInt32 ec_tmp;
  39
  40    if (lo == hi) return;
  41
  42    if (hi - lo > 3) {
  43       for ( i = hi-4; i >= lo; i-- ) {
  44          tmp = fmap[i];
  45          ec_tmp = eclass[tmp];
  46          for ( j = i+4; j <= hi && ec_tmp > eclass[fmap[j]]; j += 4 )
  47             fmap[j-4] = fmap[j];
  48          fmap[j-4] = tmp;
  49       }
  50    }
  51
  52    for ( i = hi-1; i >= lo; i-- ) {
  53       tmp = fmap[i];
  54       ec_tmp = eclass[tmp];
  55       for ( j = i+1; j <= hi && ec_tmp > eclass[fmap[j]]; j++ )
  56          fmap[j-1] = fmap[j];
  57       fmap[j-1] = tmp;
  58    }
  59 }
  60
  61
  62 /*---------------------------------------------*/
  63 #define fswap(zz1, zz2) \
  64    { Int32 zztmp = zz1; zz1 = zz2; zz2 = zztmp; }
  65
  66 #define fvswap(zzp1, zzp2, zzn)       \
  67 {                                     \
  68    Int32 yyp1 = (zzp1);               \
  69    Int32 yyp2 = (zzp2);               \
  70    Int32 yyn  = (zzn);                \
  71    while (yyn > 0) {                  \
  72       fswap(fmap[yyp1], fmap[yyp2]);  \
  73       yyp1++; yyp2++; yyn--;          \
  74    }                                  \
  75 }
  76
  77
  78 #define fmin(a,b) ((a) < (b)) ? (a) : (b)
  79
  80 #define fpush(lz,hz) { stackLo[sp] = lz; \
  81                        stackHi[sp] = hz; \
  82                        sp++; }
  83
  84 #define fpop(lz,hz) { sp--;              \
  85                       lz = stackLo[sp];  \
  86                       hz = stackHi[sp]; }
  87
  88 #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 {
  98    Int32 unLo, unHi, ltLo, gtHi, n, m;
  99    Int32 sp, lo, hi;
 100    UInt32 med, r, r3;
 101    Int32 stackLo[FALLBACK_QSORT_STACK_SIZE];
 102    Int32 stackHi[FALLBACK_QSORT_STACK_SIZE];
 103
 104    r = 0;
 105
 106    sp = 0;
 107    fpush ( loSt, hiSt );
 108
 109    while (sp > 0) {
 110
 111       AssertH ( sp < FALLBACK_QSORT_STACK_SIZE - 1, 1004 );
 112
 113       fpop ( lo, hi );
 114       if (hi - lo < FALLBACK_QSORT_SMALL_THRESH) {
 115          fallbackSimpleSort ( fmap, eclass, lo, hi );
 116          continue;
 117       }
 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       r = ((r * 7621) + 1) % 32768;
 127       r3 = r % 3;
 128       if (r3 == 0) med = eclass[fmap[lo]]; else
 129       if (r3 == 1) med = eclass[fmap[(lo+hi)>>1]]; else
 130                    med = eclass[fmap[hi]];
 131
 132       unLo = ltLo = lo;
 133       unHi = gtHi = hi;
 134
 135       while (1) {
 136          while (1) {
 137             if (unLo > unHi) break;
 138             n = (Int32)eclass[fmap[unLo]] - (Int32)med;
 139             if (n == 0) {
 140                fswap(fmap[unLo], fmap[ltLo]);
 141                ltLo++; unLo++;
 142                continue;
 143             };
 144             if (n > 0) break;
 145             unLo++;
 146          }
 147          while (1) {
 148             if (unLo > unHi) break;
 149             n = (Int32)eclass[fmap[unHi]] - (Int32)med;
 150             if (n == 0) {
 151                fswap(fmap[unHi], fmap[gtHi]);
 152                gtHi--; unHi--;
 153                continue;
 154             };
 155             if (n < 0) break;
 156             unHi--;
 157          }
 158          if (unLo > unHi) break;
 159          fswap(fmap[unLo], fmap[unHi]); unLo++; unHi--;
 160       }
 161
 162       AssertD ( unHi == unLo-1, "fallbackQSort3(2)" );
 163
 164       if (gtHi < ltLo) continue;
 165
 166       n = fmin(ltLo-lo, unLo-ltLo); fvswap(lo, unLo-n, n);
 167       m = fmin(hi-gtHi, gtHi-unHi); fvswap(unLo, hi-m+1, m);
 168
 169       n = lo + unLo - ltLo - 1;
 170       m = hi - (gtHi - unHi) + 1;
 171
 172       if (n - lo > hi - m) {
 173          fpush ( lo, n );
 174          fpush ( m, hi );
 175       } else {
 176          fpush ( m, hi );
 177          fpush ( lo, n );
 178       }
 179    }
 180 }
 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 #define       SET_BH(zz)  bhtab[(zz) >> 5] |= (1 << ((zz) & 31))
 206 #define     CLEAR_BH(zz)  bhtab[(zz) >> 5] &= ~(1 << ((zz) & 31))
 207 #define     ISSET_BH(zz)  (bhtab[(zz) >> 5] & (1 << ((zz) & 31)))
 208 #define      WORD_BH(zz)  bhtab[(zz) >> 5]
 209 #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 {
 218    Int32 ftab[257];
 219    Int32 ftabCopy[256];
 220    Int32 H, i, j, k, l, r, cc, cc1;
 221    Int32 nNotDone;
 222    Int32 nBhtab;
 223    UChar* eclass8 = (UChar*)eclass;
 224
 225    /*--
 226       Initial 1-char radix sort to generate
 227       initial fmap and initial BH bits.
 228    --*/
 229    if (verb >= 4)
 230       VPrintf0 ( "        bucket sorting ...\n" );
 231    for (i = 0; i < 257;    i++) ftab[i] = 0;
 232    for (i = 0; i < nblock; i++) ftab[eclass8[i]]++;
 233    for (i = 0; i < 256;    i++) ftabCopy[i] = ftab[i];
 234    for (i = 1; i < 257;    i++) ftab[i] += ftab[i-1];
 235
 236    for (i = 0; i < nblock; i++) {
 237       j = eclass8[i];
 238       k = ftab[j] - 1;
 239       ftab[j] = k;
 240       fmap[k] = i;
 241    }
 242
 243    nBhtab = 2 + (nblock / 32);
 244    for (i = 0; i < nBhtab; i++) bhtab[i] = 0;
 245    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    for (i = 0; i < 32; i++) {
 255       SET_BH(nblock + 2*i);
 256       CLEAR_BH(nblock + 2*i + 1);
 257    }
 258
 259    /*-- the log(N) loop --*/
 260    H = 1;
 261    while (1) {
 262
 263       if (verb >= 4)
 264          VPrintf1 ( "        depth %6d has ", H );
 265
 266       j = 0;
 267       for (i = 0; i < nblock; i++) {
 268          if (ISSET_BH(i)) j = i;
 269          k = fmap[i] - H; if (k < 0) k += nblock;
 270          eclass[k] = j;
 271       }
 272
 273       nNotDone = 0;
 274       r = -1;
 275       while (1) {
 276
 277          /*-- find the next non-singleton bucket --*/
 278          k = r + 1;
 279          while (ISSET_BH(k) && UNALIGNED_BH(k)) k++;
 280          if (ISSET_BH(k)) {
 281             while (WORD_BH(k) == 0xffffffff) k += 32;
 282             while (ISSET_BH(k)) k++;
 283          }
 284          l = k - 1;
 285          if (l >= nblock) break;
 286          while (!ISSET_BH(k) && UNALIGNED_BH(k)) k++;
 287          if (!ISSET_BH(k)) {
 288             while (WORD_BH(k) == 0x00000000) k += 32;
 289             while (!ISSET_BH(k)) k++;
 290          }
 291          r = k - 1;
 292          if (r >= nblock) break;
 293
 294          /*-- now [l, r] bracket current bucket --*/
 295          if (r > l) {
 296             nNotDone += (r - l + 1);
 297             fallbackQSort3 ( fmap, eclass, l, r );
 298
 299             /*-- scan bucket and generate header bits-- */
 300             cc = -1;
 301             for (i = l; i <= r; i++) {
 302                cc1 = eclass[fmap[i]];
 303                if (cc != cc1) { SET_BH(i); cc = cc1; };
 304             }
 305          }
 306       }
 307
 308       if (verb >= 4)
 309          VPrintf1 ( "%6d unresolved strings\n", nNotDone );
 310
 311       H *= 2;
 312       if (H > nblock || nNotDone == 0) break;
 313    }
 314
 315    /*--
 316       Reconstruct the original block in
 317       eclass8 [0 .. nblock-1], since the
 318       previous phase destroyed it.
 319    --*/
 320    if (verb >= 4)
 321       VPrintf0 ( "        reconstructing block ...\n" );
 322    j = 0;
 323    for (i = 0; i < nblock; i++) {
 324       while (ftabCopy[j] == 0) j++;
 325       ftabCopy[j]--;
 326       eclass8[fmap[i]] = (UChar)j;
 327    }
 328    AssertH ( j < 256, 1005 );
 329 }
 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 {
 354    Int32  k;
 355    UChar  c1, c2;
 356    UInt16 s1, s2;
 357
 358    AssertD ( i1 != i2, "mainGtU" );
 359    /* 1 */
 360    c1 = block[i1]; c2 = block[i2];
 361    if (c1 != c2) return (c1 > c2);
 362    i1++; i2++;
 363    /* 2 */
 364    c1 = block[i1]; c2 = block[i2];
 365    if (c1 != c2) return (c1 > c2);
 366    i1++; i2++;
 367    /* 3 */
 368    c1 = block[i1]; c2 = block[i2];
 369    if (c1 != c2) return (c1 > c2);
 370    i1++; i2++;
 371    /* 4 */
 372    c1 = block[i1]; c2 = block[i2];
 373    if (c1 != c2) return (c1 > c2);
 374    i1++; i2++;
 375    /* 5 */
 376    c1 = block[i1]; c2 = block[i2];
 377    if (c1 != c2) return (c1 > c2);
 378    i1++; i2++;
 379    /* 6 */
 380    c1 = block[i1]; c2 = block[i2];
 381    if (c1 != c2) return (c1 > c2);
 382    i1++; i2++;
 383    /* 7 */
 384    c1 = block[i1]; c2 = block[i2];
 385    if (c1 != c2) return (c1 > c2);
 386    i1++; i2++;
 387    /* 8 */
 388    c1 = block[i1]; c2 = block[i2];
 389    if (c1 != c2) return (c1 > c2);
 390    i1++; i2++;
 391    /* 9 */
 392    c1 = block[i1]; c2 = block[i2];
 393    if (c1 != c2) return (c1 > c2);
 394    i1++; i2++;
 395    /* 10 */
 396    c1 = block[i1]; c2 = block[i2];
 397    if (c1 != c2) return (c1 > c2);
 398    i1++; i2++;
 399    /* 11 */
 400    c1 = block[i1]; c2 = block[i2];
 401    if (c1 != c2) return (c1 > c2);
 402    i1++; i2++;
 403    /* 12 */
 404    c1 = block[i1]; c2 = block[i2];
 405    if (c1 != c2) return (c1 > c2);
 406    i1++; i2++;
 407
 408    k = nblock + 8;
 409
 410    do {
 411       /* 1 */
 412       c1 = block[i1]; c2 = block[i2];
 413       if (c1 != c2) return (c1 > c2);
 414       s1 = quadrant[i1]; s2 = quadrant[i2];
 415       if (s1 != s2) return (s1 > s2);
 416       i1++; i2++;
 417       /* 2 */
 418       c1 = block[i1]; c2 = block[i2];
 419       if (c1 != c2) return (c1 > c2);
 420       s1 = quadrant[i1]; s2 = quadrant[i2];
 421       if (s1 != s2) return (s1 > s2);
 422       i1++; i2++;
 423       /* 3 */
 424       c1 = block[i1]; c2 = block[i2];
 425       if (c1 != c2) return (c1 > c2);
 426       s1 = quadrant[i1]; s2 = quadrant[i2];
 427       if (s1 != s2) return (s1 > s2);
 428       i1++; i2++;
 429       /* 4 */
 430       c1 = block[i1]; c2 = block[i2];
 431       if (c1 != c2) return (c1 > c2);
 432       s1 = quadrant[i1]; s2 = quadrant[i2];
 433       if (s1 != s2) return (s1 > s2);
 434       i1++; i2++;
 435       /* 5 */
 436       c1 = block[i1]; c2 = block[i2];
 437       if (c1 != c2) return (c1 > c2);
 438       s1 = quadrant[i1]; s2 = quadrant[i2];
 439       if (s1 != s2) return (s1 > s2);
 440       i1++; i2++;
 441       /* 6 */
 442       c1 = block[i1]; c2 = block[i2];
 443       if (c1 != c2) return (c1 > c2);
 444       s1 = quadrant[i1]; s2 = quadrant[i2];
 445       if (s1 != s2) return (s1 > s2);
 446       i1++; i2++;
 447       /* 7 */
 448       c1 = block[i1]; c2 = block[i2];
 449       if (c1 != c2) return (c1 > c2);
 450       s1 = quadrant[i1]; s2 = quadrant[i2];
 451       if (s1 != s2) return (s1 > s2);
 452       i1++; i2++;
 453       /* 8 */
 454       c1 = block[i1]; c2 = block[i2];
 455       if (c1 != c2) return (c1 > c2);
 456       s1 = quadrant[i1]; s2 = quadrant[i2];
 457       if (s1 != s2) return (s1 > s2);
 458       i1++; i2++;
 459
 460       if (i1 >= nblock) i1 -= nblock;
 461       if (i2 >= nblock) i2 -= nblock;
 462
 463       k -= 8;
 464       (*budget)--;
 465    }
 466       while (k >= 0);
 467
 468    return False;
 469 }
 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 {
 494    Int32 i, j, h, bigN, hp;
 495    UInt32 v;
 496
 497    bigN = hi - lo + 1;
 498    if (bigN < 2) return;
 499
 500    hp = 0;
 501    while (incs[hp] < bigN) hp++;
 502    hp--;
 503
 504    for (; hp >= 0; hp--) {
 505       h = incs[hp];
 506
 507       i = lo + h;
 508       while (True) {
 509
 510          /*-- copy 1 --*/
 511          if (i > hi) break;
 512          v = ptr[i];
 513          j = i;
 514          while ( mainGtU (
 515                     ptr[j-h]+d, v+d, block, quadrant, nblock, budget
 516                  ) ) {
 517             ptr[j] = ptr[j-h];
 518             j = j - h;
 519             if (j <= (lo + h - 1)) break;
 520          }
 521          ptr[j] = v;
 522          i++;
 523
 524          /*-- copy 2 --*/
 525          if (i > hi) break;
 526          v = ptr[i];
 527          j = i;
 528          while ( mainGtU (
 529                     ptr[j-h]+d, v+d, block, quadrant, nblock, budget
 530                  ) ) {
 531             ptr[j] = ptr[j-h];
 532             j = j - h;
 533             if (j <= (lo + h - 1)) break;
 534          }
 535          ptr[j] = v;
 536          i++;
 537
 538          /*-- copy 3 --*/
 539          if (i > hi) break;
 540          v = ptr[i];
 541          j = i;
 542          while ( mainGtU (
 543                     ptr[j-h]+d, v+d, block, quadrant, nblock, budget
 544                  ) ) {
 545             ptr[j] = ptr[j-h];
 546             j = j - h;
 547             if (j <= (lo + h - 1)) break;
 548          }
 549          ptr[j] = v;
 550          i++;
 551
 552          if (*budget < 0) return;
 553       }
 554    }
 555 }
 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    { Int32 zztmp = zz1; zz1 = zz2; zz2 = zztmp; }
 569
 570 #define mvswap(zzp1, zzp2, zzn)       \
 571 {                                     \
 572    Int32 yyp1 = (zzp1);               \
 573    Int32 yyp2 = (zzp2);               \
 574    Int32 yyn  = (zzn);                \
 575    while (yyn > 0) {                  \
 576       mswap(ptr[yyp1], ptr[yyp2]);    \
 577       yyp1++; yyp2++; yyn--;          \
 578    }                                  \
 579 }
 580
 581 static
 582 __inline__
 583 UChar mmed3 ( UChar a, UChar b, UChar c )
 584 {
 585    UChar t;
 586    if (a > b) { t = a; a = b; b = t; };
 587    if (b > c) {
 588       b = c;
 589       if (a > b) b = a;
 590    }
 591    return b;
 592 }
 593
 594 #define mmin(a,b) ((a) < (b)) ? (a) : (b)
 595
 596 #define mpush(lz,hz,dz) { stackLo[sp] = lz; \
 597                           stackHi[sp] = hz; \
 598                           stackD [sp] = dz; \
 599                           sp++; }
 600
 601 #define mpop(lz,hz,dz) { sp--;             \
 602                          lz = stackLo[sp]; \
 603                          hz = stackHi[sp]; \
 604                          dz = stackD [sp]; }
 605
 606
 607 #define mnextsize(az) (nextHi[az]-nextLo[az])
 608
 609 #define mnextswap(az,bz)                                        \
 610    { Int32 tz;                                                  \
 611      tz = nextLo[az]; nextLo[az] = nextLo[bz]; nextLo[bz] = tz; \
 612      tz = nextHi[az]; nextHi[az] = nextHi[bz]; nextHi[bz] = tz; \
 613      tz = nextD [az]; nextD [az] = nextD [bz]; nextD [bz] = tz; }
 614
 615
 616 #define MAIN_QSORT_SMALL_THRESH 20
 617 #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 {
 630    Int32 unLo, unHi, ltLo, gtHi, n, m, med;
 631    Int32 sp, lo, hi, d;
 632
 633    Int32 stackLo[MAIN_QSORT_STACK_SIZE];
 634    Int32 stackHi[MAIN_QSORT_STACK_SIZE];
 635    Int32 stackD [MAIN_QSORT_STACK_SIZE];
 636
 637    Int32 nextLo[3];
 638    Int32 nextHi[3];
 639    Int32 nextD [3];
 640
 641    sp = 0;
 642    mpush ( loSt, hiSt, dSt );
 643
 644    while (sp > 0) {
 645
 646       AssertH ( sp < MAIN_QSORT_STACK_SIZE - 2, 1001 );
 647
 648       mpop ( lo, hi, d );
 649       if (hi - lo < MAIN_QSORT_SMALL_THRESH ||
 650           d > MAIN_QSORT_DEPTH_THRESH) {
 651          mainSimpleSort ( ptr, block, quadrant, nblock, lo, hi, d, budget );
 652          if (*budget < 0) return;
 653          continue;
 654       }
 655
 656       med = (Int32)
 657             mmed3 ( block[ptr[ lo         ]+d],
 658                     block[ptr[ hi         ]+d],
 659                     block[ptr[ (lo+hi)>>1 ]+d] );
 660
 661       unLo = ltLo = lo;
 662       unHi = gtHi = hi;
 663
 664       while (True) {
 665          while (True) {
 666             if (unLo > unHi) break;
 667             n = ((Int32)block[ptr[unLo]+d]) - med;
 668             if (n == 0) {
 669                mswap(ptr[unLo], ptr[ltLo]);
 670                ltLo++; unLo++; continue;
 671             };
 672             if (n >  0) break;
 673             unLo++;
 674          }
 675          while (True) {
 676             if (unLo > unHi) break;
 677             n = ((Int32)block[ptr[unHi]+d]) - med;
 678             if (n == 0) {
 679                mswap(ptr[unHi], ptr[gtHi]);
 680                gtHi--; unHi--; continue;
 681             };
 682             if (n <  0) break;
 683             unHi--;
 684          }
 685          if (unLo > unHi) break;
 686          mswap(ptr[unLo], ptr[unHi]); unLo++; unHi--;
 687       }
 688
 689       AssertD ( unHi == unLo-1, "mainQSort3(2)" );
 690
 691       if (gtHi < ltLo) {
 692          mpush(lo, hi, d+1 );
 693          continue;
 694       }
 695
 696       n = mmin(ltLo-lo, unLo-ltLo); mvswap(lo, unLo-n, n);
 697       m = mmin(hi-gtHi, gtHi-unHi); mvswap(unLo, hi-m+1, m);
 698
 699       n = lo + unLo - ltLo - 1;
 700       m = hi - (gtHi - unHi) + 1;
 701
 702       nextLo[0] = lo;  nextHi[0] = n;   nextD[0] = d;
 703       nextLo[1] = m;   nextHi[1] = hi;  nextD[1] = d;
 704       nextLo[2] = n+1; nextHi[2] = m-1; nextD[2] = d+1;
 705
 706       if (mnextsize(0) < mnextsize(1)) mnextswap(0,1);
 707       if (mnextsize(1) < mnextsize(2)) mnextswap(1,2);
 708       if (mnextsize(0) < mnextsize(1)) mnextswap(0,1);
 709
 710       AssertD (mnextsize(0) >= mnextsize(1), "mainQSort3(8)" );
 711       AssertD (mnextsize(1) >= mnextsize(2), "mainQSort3(9)" );
 712
 713       mpush (nextLo[0], nextHi[0], nextD[0]);
 714       mpush (nextLo[1], nextHi[1], nextD[1]);
 715       mpush (nextLo[2], nextHi[2], nextD[2]);
 716    }
 717 }
 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 #define BIGFREQ(b) (ftab[((b)+1) << 8] - ftab[(b) << 8])
 747 #define SETMASK (1 << 21)
 748 #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 {
 759    Int32  i, j, k, ss, sb;
 760    Int32  runningOrder[256];
 761    Bool   bigDone[256];
 762    Int32  copyStart[256];
 763    Int32  copyEnd  [256];
 764    UChar  c1;
 765    Int32  numQSorted;
 766    UInt16 s;
 767    if (verb >= 4) VPrintf0 ( "        main sort initialise ...\n" );
 768
 769    /*-- set up the 2-byte frequency table --*/
 770    for (i = 65536; i >= 0; i--) ftab[i] = 0;
 771
 772    j = block[0] << 8;
 773    i = nblock-1;
 774    for (; i >= 3; i -= 4) {
 775       quadrant[i] = 0;
 776       j = (j >> 8) | ( ((UInt16)block[i]) << 8);
 777       ftab[j]++;
 778       quadrant[i-1] = 0;
 779       j = (j >> 8) | ( ((UInt16)block[i-1]) << 8);
 780       ftab[j]++;
 781       quadrant[i-2] = 0;
 782       j = (j >> 8) | ( ((UInt16)block[i-2]) << 8);
 783       ftab[j]++;
 784       quadrant[i-3] = 0;
 785       j = (j >> 8) | ( ((UInt16)block[i-3]) << 8);
 786       ftab[j]++;
 787    }
 788    for (; i >= 0; i--) {
 789       quadrant[i] = 0;
 790       j = (j >> 8) | ( ((UInt16)block[i]) << 8);
 791       ftab[j]++;
 792    }
 793
 794    /*-- (emphasises close relationship of block & quadrant) --*/
 795    for (i = 0; i < BZ_N_OVERSHOOT; i++) {
 796       block   [nblock+i] = block[i];
 797       quadrant[nblock+i] = 0;
 798    }
 799
 800    if (verb >= 4) VPrintf0 ( "        bucket sorting ...\n" );
 801
 802    /*-- Complete the initial radix sort --*/
 803    for (i = 1; i <= 65536; i++) ftab[i] += ftab[i-1];
 804
 805    s = block[0] << 8;
 806    i = nblock-1;
 807    for (; i >= 3; i -= 4) {
 808       s = (s >> 8) | (block[i] << 8);
 809       j = ftab[s] -1;
 810       ftab[s] = j;
 811       ptr[j] = i;
 812       s = (s >> 8) | (block[i-1] << 8);
 813       j = ftab[s] -1;
 814       ftab[s] = j;
 815       ptr[j] = i-1;
 816       s = (s >> 8) | (block[i-2] << 8);
 817       j = ftab[s] -1;
 818       ftab[s] = j;
 819       ptr[j] = i-2;
 820       s = (s >> 8) | (block[i-3] << 8);
 821       j = ftab[s] -1;
 822       ftab[s] = j;
 823       ptr[j] = i-3;
 824    }
 825    for (; i >= 0; i--) {
 826       s = (s >> 8) | (block[i] << 8);
 827       j = ftab[s] -1;
 828       ftab[s] = j;
 829       ptr[j] = i;
 830    }
 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    for (i = 0; i <= 255; i++) {
 838       bigDone     [i] = False;
 839       runningOrder[i] = i;
 840    }
 841
 842    {
 843       Int32 vv;
 844       Int32 h = 1;
 845       do h = 3 * h + 1; while (h <= 256);
 846       do {
 847          h = h / 3;
 848          for (i = h; i <= 255; i++) {
 849             vv = runningOrder[i];
 850             j = i;
 851             while ( BIGFREQ(runningOrder[j-h]) > BIGFREQ(vv) ) {
 852                runningOrder[j] = runningOrder[j-h];
 853                j = j - h;
 854                if (j <= (h - 1)) goto zero;
 855             }
 856             zero:
 857             runningOrder[j] = vv;
 858          }
 859       } while (h != 1);
 860    }
 861
 862    /*--
 863       The main sorting loop.
 864    --*/
 865
 866    numQSorted = 0;
 867
 868    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       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       for (j = 0; j <= 255; j++) {
 887          if (j != ss) {
 888             sb = (ss << 8) + j;
 889             if ( ! (ftab[sb] & SETMASK) ) {
 890                Int32 lo = ftab[sb]   & CLEARMASK;
 891                Int32 hi = (ftab[sb+1] & CLEARMASK) - 1;
 892                if (hi > lo) {
 893                   if (verb >= 4)
 894                      VPrintf4 ( "        qsort [0x%x, 0x%x]   "
 895                                 "done %d   this %d\n",
 896                                 ss, j, numQSorted, hi - lo + 1 );
 897                   mainQSort3 (
 898                      ptr, block, quadrant, nblock,
 899                      lo, hi, BZ_N_RADIX, budget
 900                   );
 901                   numQSorted += (hi - lo + 1);
 902                   if (*budget < 0) return;
 903                }
 904             }
 905             ftab[sb] |= SETMASK;
 906          }
 907       }
 908
 909       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       {
 919          for (j = 0; j <= 255; j++) {
 920             copyStart[j] =  ftab[(j << 8) + ss]     & CLEARMASK;
 921             copyEnd  [j] = (ftab[(j << 8) + ss + 1] & CLEARMASK) - 1;
 922          }
 923          for (j = ftab[ss << 8] & CLEARMASK; j < copyStart[ss]; j++) {
 924             k = ptr[j]-1; if (k < 0) k += nblock;
 925             c1 = block[k];
 926             if (!bigDone[c1])
 927                ptr[ copyStart[c1]++ ] = k;
 928          }
 929          for (j = (ftab[(ss+1) << 8] & CLEARMASK) - 1; j > copyEnd[ss]; j--) {
 930             k = ptr[j]-1; if (k < 0) k += nblock;
 931             c1 = block[k];
 932             if (!bigDone[c1])
 933                ptr[ copyEnd[c1]-- ] = k;
 934          }
 935       }
 936
 937       AssertH ( (copyStart[ss]-1 == copyEnd[ss])
 938                 ||
 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                 (copyStart[ss] == 0 && copyEnd[ss] == nblock-1),
 944                 1007 )
 945
 946       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       bigDone[ss] = True;
 988
 989       if (i < 255) {
 990          Int32 bbStart  = ftab[ss << 8] & CLEARMASK;
 991          Int32 bbSize   = (ftab[(ss+1) << 8] & CLEARMASK) - bbStart;
 992          Int32 shifts   = 0;
 993
 994          while ((bbSize >> shifts) > 65534) shifts++;
 995
 996          for (j = bbSize-1; j >= 0; j--) {
 997             Int32 a2update     = ptr[bbStart + j];
 998             UInt16 qVal        = (UInt16)(j >> shifts);
 999             quadrant[a2update] = qVal;
1000             if (a2update < BZ_N_OVERSHOOT)
1001                quadrant[a2update + nblock] = qVal;
1002          }
1003          AssertH ( ((bbSize-1) >> shifts) <= 65535, 1002 );
1004       }
1005
1006    }
1007
1008    if (verb >= 4)
1009       VPrintf3 ( "        %d pointers, %d sorted, %d scanned\n",
1010                  nblock, numQSorted, nblock - numQSorted );
1011 }
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 {
1033    UInt32* ptr    = s->ptr;
1034    UChar*  block  = s->block;
1035    UInt32* ftab   = s->ftab;
1036    Int32   nblock = s->nblock;
1037    Int32   verb   = s->verbosity;
1038    Int32   wfact  = s->workFactor;
1039    UInt16* quadrant;
1040    Int32   budget;
1041    Int32   budgetInit;
1042    Int32   i;
1043
1044    if (nblock < 10000) {
1045       fallbackSort ( s->arr1, s->arr2, ftab, nblock, verb );
1046    } 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       i = nblock+BZ_N_OVERSHOOT;
1053       if (i & 1) i++;
1054       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       if (wfact < 1  ) wfact = 1;
1064       if (wfact > 100) wfact = 100;
1065       budgetInit = nblock * ((wfact-1) / 3);
1066       budget = budgetInit;
1067
1068       mainSort ( ptr, block, quadrant, ftab, nblock, verb, &budget );
1069       if (verb >= 3)
1070          VPrintf3 ( "      %d work, %d block, ratio %5.2f\n",
1071                     budgetInit - budget,
1072                     nblock,
1073                     (float)(budgetInit - budget) /
1074                     (float)(nblock==0 ? 1 : nblock) );
1075       if (budget < 0) {
1076          if (verb >= 2)
1077             VPrintf0 ( "    too repetitive; using fallback"
1078                        " sorting algorithm\n" );
1079          fallbackSort ( s->arr1, s->arr2, ftab, nblock, verb );
1080       }
1081    }
1082
1083    s->origPtr = -1;
1084    for (i = 0; i < s->nblock; i++)
1085       if (ptr[i] == 0)
1086          { s->origPtr = i; break; };
1087
1088    AssertH( s->origPtr != -1, 1003 );
1089 }
1090
1091
1092 /*-------------------------------------------------------------*/
1093 /*--- end                                       blocksort.c ---*/
1094 /*-------------------------------------------------------------*/