C/BwtSort.c

   1 /* BwtSort.c -- BWT block sorting\r
   2 2008-08-17\r
   3 Igor Pavlov\r
   4 Public domain */\r
   5 \r
   6 #include "BwtSort.h"\r
   7 #include "Sort.h"\r
   8 \r
   9 /* #define BLOCK_SORT_USE_HEAP_SORT */\r
  10 \r
  11 #define NO_INLINE MY_FAST_CALL\r
  12 \r
  13 /* Don't change it !!! */\r
  14 #define kNumHashBytes 2\r
  15 #define kNumHashValues (1 << (kNumHashBytes * 8))\r
  16 \r
  17 /* kNumRefBitsMax must be < (kNumHashBytes * 8) = 16 */\r
  18 #define kNumRefBitsMax 12\r
  19 \r
  20 #define BS_TEMP_SIZE kNumHashValues\r
  21 \r
  22 #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
  23 \r
  24 /* 32 Flags in UInt32 word */\r
  25 #define kNumFlagsBits 5\r
  26 #define kNumFlagsInWord (1 << kNumFlagsBits)\r
  27 #define kFlagsMask (kNumFlagsInWord - 1)\r
  28 #define kAllFlags 0xFFFFFFFF\r
  29 \r
  30 #else\r
  31 \r
  32 #define kNumBitsMax 20\r
  33 #define kIndexMask ((1 << kNumBitsMax) - 1)\r
  34 #define kNumExtraBits (32 - kNumBitsMax)\r
  35 #define kNumExtra0Bits (kNumExtraBits - 2)\r
  36 #define kNumExtra0Mask ((1 << kNumExtra0Bits) - 1)\r
  37 \r
  38 #define SetFinishedGroupSize(p, size) \\r
  39   {  *(p) |= ((((size) - 1) & kNumExtra0Mask) << kNumBitsMax); \\r
  40     if ((size) > (1 << kNumExtra0Bits)) { \\r
  41     *(p) |= 0x40000000;  *((p) + 1) |= ((((size) - 1)>> kNumExtra0Bits) << kNumBitsMax); } } \\r
  42 \r
  43 static void SetGroupSize(UInt32 *p, UInt32 size)\r
  44 {\r
  45   if (--size == 0)\r
  46     return;\r
  47   *p |= 0x80000000 | ((size & kNumExtra0Mask) << kNumBitsMax);\r
  48   if (size >= (1 << kNumExtra0Bits))\r
  49   {\r
  50     *p |= 0x40000000;\r
  51     p[1] |= ((size >> kNumExtra0Bits) << kNumBitsMax);\r
  52   }\r
  53 }\r
  54 \r
  55 #endif\r
  56 \r
  57 /*\r
  58 SortGroup - is recursive Range-Sort function with HeapSort optimization for small blocks\r
  59   "range" is not real range. It's only for optimization.\r
  60 returns: 1 - if there are groups, 0 - no more groups\r
  61 */\r
  62 \r
  63 UInt32 NO_INLINE SortGroup(UInt32 BlockSize, UInt32 NumSortedBytes, UInt32 groupOffset, UInt32 groupSize, int NumRefBits, UInt32 *Indices\r
  64   #ifndef BLOCK_SORT_USE_HEAP_SORT\r
  65   , UInt32 left, UInt32 range\r
  66   #endif\r
  67   )\r
  68 {\r
  69   UInt32 *ind2 = Indices + groupOffset;\r
  70   UInt32 *Groups;\r
  71   if (groupSize <= 1)\r
  72   {\r
  73     /*\r
  74     #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
  75     SetFinishedGroupSize(ind2, 1);\r
  76     #endif\r
  77     */\r
  78     return 0;\r
  79   }\r
  80   Groups = Indices + BlockSize + BS_TEMP_SIZE;\r
  81   if (groupSize <= ((UInt32)1 << NumRefBits)\r
  82       #ifndef BLOCK_SORT_USE_HEAP_SORT\r
  83       && groupSize <= range\r
  84       #endif\r
  85       )\r
  86   {\r
  87     UInt32 *temp = Indices + BlockSize;\r
  88     UInt32 j;\r
  89     UInt32 mask, thereAreGroups, group, cg;\r
  90     {\r
  91       UInt32 gPrev;\r
  92       UInt32 gRes = 0;\r
  93       {\r
  94         UInt32 sp = ind2[0] + NumSortedBytes;\r
  95         if (sp >= BlockSize) sp -= BlockSize;\r
  96         gPrev = Groups[sp];\r
  97         temp[0] = (gPrev << NumRefBits);\r
  98       }\r
  99       \r
 100       for (j = 1; j < groupSize; j++)\r
 101       {\r
 102         UInt32 sp = ind2[j] + NumSortedBytes;\r
 103         UInt32 g;\r
 104         if (sp >= BlockSize) sp -= BlockSize;\r
 105         g = Groups[sp];\r
 106         temp[j] = (g << NumRefBits) | j;\r
 107         gRes |= (gPrev ^ g);\r
 108       }\r
 109       if (gRes == 0)\r
 110       {\r
 111         #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 112         SetGroupSize(ind2, groupSize);\r
 113         #endif\r
 114         return 1;\r
 115       }\r
 116     }\r
 117     \r
 118     HeapSort(temp, groupSize);\r
 119     mask = ((1 << NumRefBits) - 1);\r
 120     thereAreGroups = 0;\r
 121     \r
 122     group = groupOffset;\r
 123     cg = (temp[0] >> NumRefBits);\r
 124     temp[0] = ind2[temp[0] & mask];\r
 125 \r
 126     {\r
 127     #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 128     UInt32 *Flags = Groups + BlockSize;\r
 129     #else\r
 130     UInt32 prevGroupStart = 0;\r
 131     #endif\r
 132     \r
 133     for (j = 1; j < groupSize; j++)\r
 134     {\r
 135       UInt32 val = temp[j];\r
 136       UInt32 cgCur = (val >> NumRefBits);\r
 137       \r
 138       if (cgCur != cg)\r
 139       {\r
 140         cg = cgCur;\r
 141         group = groupOffset + j;\r
 142 \r
 143         #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 144         {\r
 145         UInt32 t = group - 1;\r
 146         Flags[t >> kNumFlagsBits] &= ~(1 << (t & kFlagsMask));\r
 147         }\r
 148         #else\r
 149         SetGroupSize(temp + prevGroupStart, j - prevGroupStart);\r
 150         prevGroupStart = j;\r
 151         #endif\r
 152       }\r
 153       else\r
 154         thereAreGroups = 1;\r
 155       {\r
 156       UInt32 ind = ind2[val & mask];\r
 157       temp[j] = ind;\r
 158       Groups[ind] = group;\r
 159       }\r
 160     }\r
 161 \r
 162     #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 163     SetGroupSize(temp + prevGroupStart, j - prevGroupStart);\r
 164     #endif\r
 165     }\r
 166 \r
 167     for (j = 0; j < groupSize; j++)\r
 168       ind2[j] = temp[j];\r
 169     return thereAreGroups;\r
 170   }\r
 171 \r
 172   /* Check that all strings are in one group (cannot sort) */\r
 173   {\r
 174     UInt32 group, j;\r
 175     UInt32 sp = ind2[0] + NumSortedBytes; if (sp >= BlockSize) sp -= BlockSize;\r
 176     group = Groups[sp];\r
 177     for (j = 1; j < groupSize; j++)\r
 178     {\r
 179       sp = ind2[j] + NumSortedBytes; if (sp >= BlockSize) sp -= BlockSize;\r
 180       if (Groups[sp] != group)\r
 181         break;\r
 182     }\r
 183     if (j == groupSize)\r
 184     {\r
 185       #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 186       SetGroupSize(ind2, groupSize);\r
 187       #endif\r
 188       return 1;\r
 189     }\r
 190   }\r
 191 \r
 192   #ifndef BLOCK_SORT_USE_HEAP_SORT\r
 193   {\r
 194   /* ---------- Range Sort ---------- */\r
 195   UInt32 i;\r
 196   UInt32 mid;\r
 197   for (;;)\r
 198   {\r
 199     UInt32 j;\r
 200     if (range <= 1)\r
 201     {\r
 202       #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 203       SetGroupSize(ind2, groupSize);\r
 204       #endif\r
 205       return 1;\r
 206     }\r
 207     mid = left + ((range + 1) >> 1);\r
 208     j = groupSize;\r
 209     i = 0;\r
 210     do\r
 211     {\r
 212       UInt32 sp = ind2[i] + NumSortedBytes; if (sp >= BlockSize) sp -= BlockSize;\r
 213       if (Groups[sp] >= mid)\r
 214       {\r
 215         for (j--; j > i; j--)\r
 216         {\r
 217           sp = ind2[j] + NumSortedBytes; if (sp >= BlockSize) sp -= BlockSize;\r
 218           if (Groups[sp] < mid)\r
 219           {\r
 220             UInt32 temp = ind2[i]; ind2[i] = ind2[j]; ind2[j] = temp;\r
 221             break;\r
 222           }\r
 223         }\r
 224         if (i >= j)\r
 225           break;\r
 226       }\r
 227     }\r
 228     while (++i < j);\r
 229     if (i == 0)\r
 230     {\r
 231       range = range - (mid - left);\r
 232       left = mid;\r
 233     }\r
 234     else if (i == groupSize)\r
 235       range = (mid - left);\r
 236     else\r
 237       break;\r
 238   }\r
 239 \r
 240   #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 241   {\r
 242     UInt32 t = (groupOffset + i - 1);\r
 243     UInt32 *Flags = Groups + BlockSize;\r
 244     Flags[t >> kNumFlagsBits] &= ~(1 << (t & kFlagsMask));\r
 245   }\r
 246   #endif\r
 247 \r
 248   {\r
 249     UInt32 j;\r
 250     for (j = i; j < groupSize; j++)\r
 251       Groups[ind2[j]] = groupOffset + i;\r
 252   }\r
 253 \r
 254   {\r
 255   UInt32 res = SortGroup(BlockSize, NumSortedBytes, groupOffset, i, NumRefBits, Indices, left, mid - left);\r
 256   return res | SortGroup(BlockSize, NumSortedBytes, groupOffset + i, groupSize - i, NumRefBits, Indices, mid, range - (mid - left));\r
 257   }\r
 258 \r
 259   }\r
 260 \r
 261   #else\r
 262 \r
 263   /* ---------- Heap Sort ---------- */\r
 264 \r
 265   {\r
 266     UInt32 j;\r
 267     for (j = 0; j < groupSize; j++)\r
 268     {\r
 269       UInt32 sp = ind2[j] + NumSortedBytes; if (sp >= BlockSize) sp -= BlockSize;\r
 270       ind2[j] = sp;\r
 271     }\r
 272 \r
 273     HeapSortRef(ind2, Groups, groupSize);\r
 274 \r
 275     /* Write Flags */\r
 276     {\r
 277     UInt32 sp = ind2[0];\r
 278     UInt32 group = Groups[sp];\r
 279 \r
 280     #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 281     UInt32 *Flags = Groups + BlockSize;\r
 282     #else\r
 283     UInt32 prevGroupStart = 0;\r
 284     #endif\r
 285 \r
 286     for (j = 1; j < groupSize; j++)\r
 287     {\r
 288       sp = ind2[j];\r
 289       if (Groups[sp] != group)\r
 290       {\r
 291         group = Groups[sp];\r
 292         #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 293         {\r
 294         UInt32 t = groupOffset + j - 1;\r
 295         Flags[t >> kNumFlagsBits] &= ~(1 << (t & kFlagsMask));\r
 296         }\r
 297         #else\r
 298         SetGroupSize(ind2 + prevGroupStart, j - prevGroupStart);\r
 299         prevGroupStart = j;\r
 300         #endif\r
 301       }\r
 302     }\r
 303 \r
 304     #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 305     SetGroupSize(ind2 + prevGroupStart, j - prevGroupStart);\r
 306     #endif\r
 307     }\r
 308     {\r
 309     /* Write new Groups values and Check that there are groups */\r
 310     UInt32 thereAreGroups = 0;\r
 311     for (j = 0; j < groupSize; j++)\r
 312     {\r
 313       UInt32 group = groupOffset + j;\r
 314       #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 315       UInt32 subGroupSize = ((ind2[j] & ~0xC0000000) >> kNumBitsMax);\r
 316       if ((ind2[j] & 0x40000000) != 0)\r
 317         subGroupSize += ((ind2[j + 1] >> kNumBitsMax) << kNumExtra0Bits);\r
 318       subGroupSize++;\r
 319       for (;;)\r
 320       {\r
 321         UInt32 original = ind2[j];\r
 322         UInt32 sp = original & kIndexMask;\r
 323         if (sp < NumSortedBytes) sp += BlockSize; sp -= NumSortedBytes;\r
 324         ind2[j] = sp | (original & ~kIndexMask);\r
 325         Groups[sp] = group;\r
 326         if (--subGroupSize == 0)\r
 327           break;\r
 328         j++;\r
 329         thereAreGroups = 1;\r
 330       }\r
 331       #else\r
 332       UInt32 *Flags = Groups + BlockSize;\r
 333       for (;;)\r
 334       {\r
 335         UInt32 sp = ind2[j]; if (sp < NumSortedBytes) sp += BlockSize; sp -= NumSortedBytes;\r
 336         ind2[j] = sp;\r
 337         Groups[sp] = group;\r
 338         if ((Flags[(groupOffset + j) >> kNumFlagsBits] & (1 << ((groupOffset + j) & kFlagsMask))) == 0)\r
 339           break;\r
 340         j++;\r
 341         thereAreGroups = 1;\r
 342       }\r
 343       #endif\r
 344     }\r
 345     return thereAreGroups;\r
 346     }\r
 347   }\r
 348   #endif\r
 349 }\r
 350 \r
 351 /* conditions: blockSize > 0 */\r
 352 UInt32 BlockSort(UInt32 *Indices, const Byte *data, UInt32 blockSize)\r
 353 {\r
 354   UInt32 *counters = Indices + blockSize;\r
 355   UInt32 i;\r
 356   UInt32 *Groups;\r
 357   #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 358   UInt32 *Flags;\r
 359   #endif\r
 360 \r
 361   /* Radix-Sort for 2 bytes */\r
 362   for (i = 0; i < kNumHashValues; i++)\r
 363     counters[i] = 0;\r
 364   for (i = 0; i < blockSize - 1; i++)\r
 365     counters[((UInt32)data[i] << 8) | data[i + 1]]++;\r
 366   counters[((UInt32)data[i] << 8) | data[0]]++;\r
 367 \r
 368   Groups = counters + BS_TEMP_SIZE;\r
 369   #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 370   Flags = Groups + blockSize;\r
 371     {\r
 372       UInt32 numWords = (blockSize + kFlagsMask) >> kNumFlagsBits;\r
 373       for (i = 0; i < numWords; i++)\r
 374         Flags[i] = kAllFlags;\r
 375     }\r
 376   #endif\r
 377 \r
 378   {\r
 379     UInt32 sum = 0;\r
 380     for (i = 0; i < kNumHashValues; i++)\r
 381     {\r
 382       UInt32 groupSize = counters[i];\r
 383       if (groupSize > 0)\r
 384       {\r
 385         #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 386         UInt32 t = sum + groupSize - 1;\r
 387         Flags[t >> kNumFlagsBits] &= ~(1 << (t & kFlagsMask));\r
 388         #endif\r
 389         sum += groupSize;\r
 390       }\r
 391       counters[i] = sum - groupSize;\r
 392     }\r
 393 \r
 394     for (i = 0; i < blockSize - 1; i++)\r
 395       Groups[i] = counters[((UInt32)data[i] << 8) | data[i + 1]];\r
 396     Groups[i] = counters[((UInt32)data[i] << 8) | data[0]];\r
 397 \r
 398     for (i = 0; i < blockSize - 1; i++)\r
 399       Indices[counters[((UInt32)data[i] << 8) | data[i + 1]]++] = i;\r
 400     Indices[counters[((UInt32)data[i] << 8) | data[0]]++] = i;\r
 401     \r
 402     #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 403     {\r
 404     UInt32 prev = 0;\r
 405     for (i = 0; i < kNumHashValues; i++)\r
 406     {\r
 407       UInt32 prevGroupSize = counters[i] - prev;\r
 408       if (prevGroupSize == 0)\r
 409         continue;\r
 410       SetGroupSize(Indices + prev, prevGroupSize);\r
 411       prev = counters[i];\r
 412     }\r
 413     }\r
 414     #endif\r
 415   }\r
 416 \r
 417   {\r
 418   int NumRefBits;\r
 419   UInt32 NumSortedBytes;\r
 420   for (NumRefBits = 0; ((blockSize - 1) >> NumRefBits) != 0; NumRefBits++);\r
 421   NumRefBits = 32 - NumRefBits;\r
 422   if (NumRefBits > kNumRefBitsMax)\r
 423     NumRefBits = kNumRefBitsMax;\r
 424 \r
 425   for (NumSortedBytes = kNumHashBytes; ; NumSortedBytes <<= 1)\r
 426   {\r
 427     #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 428     UInt32 finishedGroupSize = 0;\r
 429     #endif\r
 430     UInt32 newLimit = 0;\r
 431     for (i = 0; i < blockSize;)\r
 432     {\r
 433       UInt32 groupSize;\r
 434       #ifdef BLOCK_SORT_EXTERNAL_FLAGS\r
 435 \r
 436       if ((Flags[i >> kNumFlagsBits] & (1 << (i & kFlagsMask))) == 0)\r
 437       {\r
 438         i++;\r
 439         continue;\r
 440       }\r
 441       for (groupSize = 1;\r
 442         (Flags[(i + groupSize) >> kNumFlagsBits] & (1 << ((i + groupSize) & kFlagsMask))) != 0;\r
 443         groupSize++);\r
 444       \r
 445       groupSize++;\r
 446 \r
 447       #else\r
 448 \r
 449       groupSize = ((Indices[i] & ~0xC0000000) >> kNumBitsMax);\r
 450       {\r
 451       Bool finishedGroup = ((Indices[i] & 0x80000000) == 0);\r
 452       if ((Indices[i] & 0x40000000) != 0)\r
 453       {\r
 454         groupSize += ((Indices[i + 1] >> kNumBitsMax) << kNumExtra0Bits);\r
 455         Indices[i + 1] &= kIndexMask;\r
 456       }\r
 457       Indices[i] &= kIndexMask;\r
 458       groupSize++;\r
 459       if (finishedGroup || groupSize == 1)\r
 460       {\r
 461         Indices[i - finishedGroupSize] &= kIndexMask;\r
 462         if (finishedGroupSize > 1)\r
 463           Indices[i - finishedGroupSize + 1] &= kIndexMask;\r
 464         {\r
 465         UInt32 newGroupSize = groupSize + finishedGroupSize;\r
 466         SetFinishedGroupSize(Indices + i - finishedGroupSize, newGroupSize);\r
 467         finishedGroupSize = newGroupSize;\r
 468         }\r
 469         i += groupSize;\r
 470         continue;\r
 471       }\r
 472       finishedGroupSize = 0;\r
 473       }\r
 474 \r
 475       #endif\r
 476       \r
 477       if (NumSortedBytes >= blockSize)\r
 478       {\r
 479         UInt32 j;\r
 480         for (j = 0; j < groupSize; j++)\r
 481         {\r
 482           UInt32 t = (i + j);\r
 483           /* Flags[t >> kNumFlagsBits] &= ~(1 << (t & kFlagsMask)); */\r
 484           Groups[Indices[t]] = t;\r
 485         }\r
 486       }\r
 487       else\r
 488         if (SortGroup(blockSize, NumSortedBytes, i, groupSize, NumRefBits, Indices\r
 489           #ifndef BLOCK_SORT_USE_HEAP_SORT\r
 490           , 0, blockSize\r
 491           #endif\r
 492           ) != 0)\r
 493           newLimit = i + groupSize;\r
 494       i += groupSize;\r
 495     }\r
 496     if (newLimit == 0)\r
 497       break;\r
 498   }\r
 499   }\r
 500   #ifndef BLOCK_SORT_EXTERNAL_FLAGS\r
 501   for (i = 0; i < blockSize;)\r
 502   {\r
 503     UInt32 groupSize = ((Indices[i] & ~0xC0000000) >> kNumBitsMax);\r
 504     if ((Indices[i] & 0x40000000) != 0)\r
 505     {\r
 506       groupSize += ((Indices[i + 1] >> kNumBitsMax) << kNumExtra0Bits);\r
 507       Indices[i + 1] &= kIndexMask;\r
 508     }\r
 509     Indices[i] &= kIndexMask;\r
 510     groupSize++;\r
 511     i += groupSize;\r
 512   }\r
 513   #endif\r
 514   return Groups[0];\r
 515 }\r
 516 \r