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