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 /*-------------------------------------------------------------*/