--- /dev/null
+Intro\r
+-----\r
+This describes an adaptive, stable, natural mergesort, modestly called\r
+timsort (hey, I earned it <wink>). It has supernatural performance on many\r
+kinds of partially ordered arrays (less than lg(N!) comparisons needed, and\r
+as few as N-1), yet as fast as Python's previous highly tuned samplesort\r
+hybrid on random arrays.\r
+\r
+In a nutshell, the main routine marches over the array once, left to right,\r
+alternately identifying the next run, then merging it into the previous\r
+runs "intelligently". Everything else is complication for speed, and some\r
+hard-won measure of memory efficiency.\r
+\r
+\r
+Comparison with Python's Samplesort Hybrid\r
+------------------------------------------\r
++ timsort can require a temp array containing as many as N//2 pointers,\r
+ which means as many as 2*N extra bytes on 32-bit boxes. It can be\r
+ expected to require a temp array this large when sorting random data; on\r
+ data with significant structure, it may get away without using any extra\r
+ heap memory. This appears to be the strongest argument against it, but\r
+ compared to the size of an object, 2 temp bytes worst-case (also expected-\r
+ case for random data) doesn't scare me much.\r
+\r
+ It turns out that Perl is moving to a stable mergesort, and the code for\r
+ that appears always to require a temp array with room for at least N\r
+ pointers. (Note that I wouldn't want to do that even if space weren't an\r
+ issue; I believe its efforts at memory frugality also save timsort\r
+ significant pointer-copying costs, and allow it to have a smaller working\r
+ set.)\r
+\r
++ Across about four hours of generating random arrays, and sorting them\r
+ under both methods, samplesort required about 1.5% more comparisons\r
+ (the program is at the end of this file).\r
+\r
++ In real life, this may be faster or slower on random arrays than\r
+ samplesort was, depending on platform quirks. Since it does fewer\r
+ comparisons on average, it can be expected to do better the more\r
+ expensive a comparison function is. OTOH, it does more data movement\r
+ (pointer copying) than samplesort, and that may negate its small\r
+ comparison advantage (depending on platform quirks) unless comparison\r
+ is very expensive.\r
+\r
++ On arrays with many kinds of pre-existing order, this blows samplesort out\r
+ of the water. It's significantly faster than samplesort even on some\r
+ cases samplesort was special-casing the snot out of. I believe that lists\r
+ very often do have exploitable partial order in real life, and this is the\r
+ strongest argument in favor of timsort (indeed, samplesort's special cases\r
+ for extreme partial order are appreciated by real users, and timsort goes\r
+ much deeper than those, in particular naturally covering every case where\r
+ someone has suggested "and it would be cool if list.sort() had a special\r
+ case for this too ... and for that ...").\r
+\r
++ Here are exact comparison counts across all the tests in sortperf.py,\r
+ when run with arguments "15 20 1".\r
+\r
+ Column Key:\r
+ *sort: random data\r
+ \sort: descending data\r
+ /sort: ascending data\r
+ 3sort: ascending, then 3 random exchanges\r
+ +sort: ascending, then 10 random at the end\r
+ %sort: ascending, then randomly replace 1% of elements w/ random values\r
+ ~sort: many duplicates\r
+ =sort: all equal\r
+ !sort: worst case scenario\r
+\r
+ First the trivial cases, trivial for samplesort because it special-cased\r
+ them, and trivial for timsort because it naturally works on runs. Within\r
+ an "n" block, the first line gives the # of compares done by samplesort,\r
+ the second line by timsort, and the third line is the percentage by\r
+ which the samplesort count exceeds the timsort count:\r
+\r
+ n \sort /sort =sort\r
+------- ------ ------ ------\r
+ 32768 32768 32767 32767 samplesort\r
+ 32767 32767 32767 timsort\r
+ 0.00% 0.00% 0.00% (samplesort - timsort) / timsort\r
+\r
+ 65536 65536 65535 65535\r
+ 65535 65535 65535\r
+ 0.00% 0.00% 0.00%\r
+\r
+ 131072 131072 131071 131071\r
+ 131071 131071 131071\r
+ 0.00% 0.00% 0.00%\r
+\r
+ 262144 262144 262143 262143\r
+ 262143 262143 262143\r
+ 0.00% 0.00% 0.00%\r
+\r
+ 524288 524288 524287 524287\r
+ 524287 524287 524287\r
+ 0.00% 0.00% 0.00%\r
+\r
+1048576 1048576 1048575 1048575\r
+ 1048575 1048575 1048575\r
+ 0.00% 0.00% 0.00%\r
+\r
+ The algorithms are effectively identical in these cases, except that\r
+ timsort does one less compare in \sort.\r
+\r
+ Now for the more interesting cases. Where lg(x) is the logarithm of x to\r
+ the base 2 (e.g., lg(8)=3), lg(n!) is the information-theoretic limit for\r
+ the best any comparison-based sorting algorithm can do on average (across\r
+ all permutations). When a method gets significantly below that, it's\r
+ either astronomically lucky, or is finding exploitable structure in the\r
+ data.\r
+\r
+\r
+ n lg(n!) *sort 3sort +sort %sort ~sort !sort\r
+------- ------- ------ ------- ------- ------ ------- --------\r
+ 32768 444255 453096 453614 32908 452871 130491 469141 old\r
+ 448885 33016 33007 50426 182083 65534 new\r
+ 0.94% 1273.92% -0.30% 798.09% -28.33% 615.87% %ch from new\r
+\r
+ 65536 954037 972699 981940 65686 973104 260029 1004607\r
+ 962991 65821 65808 101667 364341 131070\r
+ 1.01% 1391.83% -0.19% 857.15% -28.63% 666.47%\r
+\r
+ 131072 2039137 2101881 2091491 131232 2092894 554790 2161379\r
+ 2057533 131410 131361 206193 728871 262142\r
+ 2.16% 1491.58% -0.10% 915.02% -23.88% 724.51%\r
+\r
+ 262144 4340409 4464460 4403233 262314 4445884 1107842 4584560\r
+ 4377402 262437 262459 416347 1457945 524286\r
+ 1.99% 1577.82% -0.06% 967.83% -24.01% 774.44%\r
+\r
+ 524288 9205096 9453356 9408463 524468 9441930 2218577 9692015\r
+ 9278734 524580 524633 837947 2916107 1048574\r
+ 1.88% 1693.52% -0.03% 1026.79% -23.92% 824.30%\r
+\r
+1048576 19458756 19950272 19838588 1048766 19912134 4430649 20434212\r
+ 19606028 1048958 1048941 1694896 5832445 2097150\r
+ 1.76% 1791.27% -0.02% 1074.83% -24.03% 874.38%\r
+\r
+ Discussion of cases:\r
+\r
+ *sort: There's no structure in random data to exploit, so the theoretical\r
+ limit is lg(n!). Both methods get close to that, and timsort is hugging\r
+ it (indeed, in a *marginal* sense, it's a spectacular improvement --\r
+ there's only about 1% left before hitting the wall, and timsort knows\r
+ darned well it's doing compares that won't pay on random data -- but so\r
+ does the samplesort hybrid). For contrast, Hoare's original random-pivot\r
+ quicksort does about 39% more compares than the limit, and the median-of-3\r
+ variant about 19% more.\r
+\r
+ 3sort, %sort, and !sort: No contest; there's structure in this data, but\r
+ not of the specific kinds samplesort special-cases. Note that structure\r
+ in !sort wasn't put there on purpose -- it was crafted as a worst case for\r
+ a previous quicksort implementation. That timsort nails it came as a\r
+ surprise to me (although it's obvious in retrospect).\r
+\r
+ +sort: samplesort special-cases this data, and does a few less compares\r
+ than timsort. However, timsort runs this case significantly faster on all\r
+ boxes we have timings for, because timsort is in the business of merging\r
+ runs efficiently, while samplesort does much more data movement in this\r
+ (for it) special case.\r
+\r
+ ~sort: samplesort's special cases for large masses of equal elements are\r
+ extremely effective on ~sort's specific data pattern, and timsort just\r
+ isn't going to get close to that, despite that it's clearly getting a\r
+ great deal of benefit out of the duplicates (the # of compares is much less\r
+ than lg(n!)). ~sort has a perfectly uniform distribution of just 4\r
+ distinct values, and as the distribution gets more skewed, samplesort's\r
+ equal-element gimmicks become less effective, while timsort's adaptive\r
+ strategies find more to exploit; in a database supplied by Kevin Altis, a\r
+ sort on its highly skewed "on which stock exchange does this company's\r
+ stock trade?" field ran over twice as fast under timsort.\r
+\r
+ However, despite that timsort does many more comparisons on ~sort, and\r
+ that on several platforms ~sort runs highly significantly slower under\r
+ timsort, on other platforms ~sort runs highly significantly faster under\r
+ timsort. No other kind of data has shown this wild x-platform behavior,\r
+ and we don't have an explanation for it. The only thing I can think of\r
+ that could transform what "should be" highly significant slowdowns into\r
+ highly significant speedups on some boxes are catastrophic cache effects\r
+ in samplesort.\r
+\r
+ But timsort "should be" slower than samplesort on ~sort, so it's hard\r
+ to count that it isn't on some boxes as a strike against it <wink>.\r
+\r
++ Here's the highwater mark for the number of heap-based temp slots (4\r
+ bytes each on this box) needed by each test, again with arguments\r
+ "15 20 1":\r
+\r
+ 2**i *sort \sort /sort 3sort +sort %sort ~sort =sort !sort\r
+ 32768 16384 0 0 6256 0 10821 12288 0 16383\r
+ 65536 32766 0 0 21652 0 31276 24576 0 32767\r
+ 131072 65534 0 0 17258 0 58112 49152 0 65535\r
+ 262144 131072 0 0 35660 0 123561 98304 0 131071\r
+ 524288 262142 0 0 31302 0 212057 196608 0 262143\r
+1048576 524286 0 0 312438 0 484942 393216 0 524287\r
+\r
+ Discussion: The tests that end up doing (close to) perfectly balanced\r
+ merges (*sort, !sort) need all N//2 temp slots (or almost all). ~sort\r
+ also ends up doing balanced merges, but systematically benefits a lot from\r
+ the preliminary pre-merge searches described under "Merge Memory" later.\r
+ %sort approaches having a balanced merge at the end because the random\r
+ selection of elements to replace is expected to produce an out-of-order\r
+ element near the midpoint. \sort, /sort, =sort are the trivial one-run\r
+ cases, needing no merging at all. +sort ends up having one very long run\r
+ and one very short, and so gets all the temp space it needs from the small\r
+ temparray member of the MergeState struct (note that the same would be\r
+ true if the new random elements were prefixed to the sorted list instead,\r
+ but not if they appeared "in the middle"). 3sort approaches N//3 temp\r
+ slots twice, but the run lengths that remain after 3 random exchanges\r
+ clearly has very high variance.\r
+\r
+\r
+A detailed description of timsort follows.\r
+\r
+Runs\r
+----\r
+count_run() returns the # of elements in the next run. A run is either\r
+"ascending", which means non-decreasing:\r
+\r
+ a0 <= a1 <= a2 <= ...\r
+\r
+or "descending", which means strictly decreasing:\r
+\r
+ a0 > a1 > a2 > ...\r
+\r
+Note that a run is always at least 2 long, unless we start at the array's\r
+last element.\r
+\r
+The definition of descending is strict, because the main routine reverses\r
+a descending run in-place, transforming a descending run into an ascending\r
+run. Reversal is done via the obvious fast "swap elements starting at each\r
+end, and converge at the middle" method, and that can violate stability if\r
+the slice contains any equal elements. Using a strict definition of\r
+descending ensures that a descending run contains distinct elements.\r
+\r
+If an array is random, it's very unlikely we'll see long runs. If a natural\r
+run contains less than minrun elements (see next section), the main loop\r
+artificially boosts it to minrun elements, via a stable binary insertion sort\r
+applied to the right number of array elements following the short natural\r
+run. In a random array, *all* runs are likely to be minrun long as a\r
+result. This has two primary good effects:\r
+\r
+1. Random data strongly tends then toward perfectly balanced (both runs have\r
+ the same length) merges, which is the most efficient way to proceed when\r
+ data is random.\r
+\r
+2. Because runs are never very short, the rest of the code doesn't make\r
+ heroic efforts to shave a few cycles off per-merge overheads. For\r
+ example, reasonable use of function calls is made, rather than trying to\r
+ inline everything. Since there are no more than N/minrun runs to begin\r
+ with, a few "extra" function calls per merge is barely measurable.\r
+\r
+\r
+Computing minrun\r
+----------------\r
+If N < 64, minrun is N. IOW, binary insertion sort is used for the whole\r
+array then; it's hard to beat that given the overheads of trying something\r
+fancier (see note BINSORT).\r
+\r
+When N is a power of 2, testing on random data showed that minrun values of\r
+16, 32, 64 and 128 worked about equally well. At 256 the data-movement cost\r
+in binary insertion sort clearly hurt, and at 8 the increase in the number\r
+of function calls clearly hurt. Picking *some* power of 2 is important\r
+here, so that the merges end up perfectly balanced (see next section). We\r
+pick 32 as a good value in the sweet range; picking a value at the low end\r
+allows the adaptive gimmicks more opportunity to exploit shorter natural\r
+runs.\r
+\r
+Because sortperf.py only tries powers of 2, it took a long time to notice\r
+that 32 isn't a good choice for the general case! Consider N=2112:\r
+\r
+>>> divmod(2112, 32)\r
+(66, 0)\r
+>>>\r
+\r
+If the data is randomly ordered, we're very likely to end up with 66 runs\r
+each of length 32. The first 64 of these trigger a sequence of perfectly\r
+balanced merges (see next section), leaving runs of lengths 2048 and 64 to\r
+merge at the end. The adaptive gimmicks can do that with fewer than 2048+64\r
+compares, but it's still more compares than necessary, and-- mergesort's\r
+bugaboo relative to samplesort --a lot more data movement (O(N) copies just\r
+to get 64 elements into place).\r
+\r
+If we take minrun=33 in this case, then we're very likely to end up with 64\r
+runs each of length 33, and then all merges are perfectly balanced. Better!\r
+\r
+What we want to avoid is picking minrun such that in\r
+\r
+ q, r = divmod(N, minrun)\r
+\r
+q is a power of 2 and r>0 (then the last merge only gets r elements into\r
+place, and r < minrun is small compared to N), or q a little larger than a\r
+power of 2 regardless of r (then we've got a case similar to "2112", again\r
+leaving too little work for the last merge to do).\r
+\r
+Instead we pick a minrun in range(32, 65) such that N/minrun is exactly a\r
+power of 2, or if that isn't possible, is close to, but strictly less than,\r
+a power of 2. This is easier to do than it may sound: take the first 6\r
+bits of N, and add 1 if any of the remaining bits are set. In fact, that\r
+rule covers every case in this section, including small N and exact powers\r
+of 2; merge_compute_minrun() is a deceptively simple function.\r
+\r
+\r
+The Merge Pattern\r
+-----------------\r
+In order to exploit regularities in the data, we're merging on natural\r
+run lengths, and they can become wildly unbalanced. That's a Good Thing\r
+for this sort! It means we have to find a way to manage an assortment of\r
+potentially very different run lengths, though.\r
+\r
+Stability constrains permissible merging patterns. For example, if we have\r
+3 consecutive runs of lengths\r
+\r
+ A:10000 B:20000 C:10000\r
+\r
+we dare not merge A with C first, because if A, B and C happen to contain\r
+a common element, it would get out of order wrt its occurrence(s) in B. The\r
+merging must be done as (A+B)+C or A+(B+C) instead.\r
+\r
+So merging is always done on two consecutive runs at a time, and in-place,\r
+although this may require some temp memory (more on that later).\r
+\r
+When a run is identified, its base address and length are pushed on a stack\r
+in the MergeState struct. merge_collapse() is then called to see whether it\r
+should merge it with preceding run(s). We would like to delay merging as\r
+long as possible in order to exploit patterns that may come up later, but we\r
+like even more to do merging as soon as possible to exploit that the run just\r
+found is still high in the memory hierarchy. We also can't delay merging\r
+"too long" because it consumes memory to remember the runs that are still\r
+unmerged, and the stack has a fixed size.\r
+\r
+What turned out to be a good compromise maintains two invariants on the\r
+stack entries, where A, B and C are the lengths of the three righmost not-yet\r
+merged slices:\r
+\r
+1. A > B+C\r
+2. B > C\r
+\r
+Note that, by induction, #2 implies the lengths of pending runs form a\r
+decreasing sequence. #1 implies that, reading the lengths right to left,\r
+the pending-run lengths grow at least as fast as the Fibonacci numbers.\r
+Therefore the stack can never grow larger than about log_base_phi(N) entries,\r
+where phi = (1+sqrt(5))/2 ~= 1.618. Thus a small # of stack slots suffice\r
+for very large arrays.\r
+\r
+If A <= B+C, the smaller of A and C is merged with B (ties favor C, for the\r
+freshness-in-cache reason), and the new run replaces the A,B or B,C entries;\r
+e.g., if the last 3 entries are\r
+\r
+ A:30 B:20 C:10\r
+\r
+then B is merged with C, leaving\r
+\r
+ A:30 BC:30\r
+\r
+on the stack. Or if they were\r
+\r
+ A:500 B:400: C:1000\r
+\r
+then A is merged with B, leaving\r
+\r
+ AB:900 C:1000\r
+\r
+on the stack.\r
+\r
+In both examples, the stack configuration after the merge still violates\r
+invariant #2, and merge_collapse() goes on to continue merging runs until\r
+both invariants are satisfied. As an extreme case, suppose we didn't do the\r
+minrun gimmick, and natural runs were of lengths 128, 64, 32, 16, 8, 4, 2,\r
+and 2. Nothing would get merged until the final 2 was seen, and that would\r
+trigger 7 perfectly balanced merges.\r
+\r
+The thrust of these rules when they trigger merging is to balance the run\r
+lengths as closely as possible, while keeping a low bound on the number of\r
+runs we have to remember. This is maximally effective for random data,\r
+where all runs are likely to be of (artificially forced) length minrun, and\r
+then we get a sequence of perfectly balanced merges (with, perhaps, some\r
+oddballs at the end).\r
+\r
+OTOH, one reason this sort is so good for partly ordered data has to do\r
+with wildly unbalanced run lengths.\r
+\r
+\r
+Merge Memory\r
+------------\r
+Merging adjacent runs of lengths A and B in-place, and in linear time, is\r
+difficult. Theoretical constructions are known that can do it, but they're\r
+too difficult and slow for practical use. But if we have temp memory equal\r
+to min(A, B), it's easy.\r
+\r
+If A is smaller (function merge_lo), copy A to a temp array, leave B alone,\r
+and then we can do the obvious merge algorithm left to right, from the temp\r
+area and B, starting the stores into where A used to live. There's always a\r
+free area in the original area comprising a number of elements equal to the\r
+number not yet merged from the temp array (trivially true at the start;\r
+proceed by induction). The only tricky bit is that if a comparison raises an\r
+exception, we have to remember to copy the remaining elements back in from\r
+the temp area, lest the array end up with duplicate entries from B. But\r
+that's exactly the same thing we need to do if we reach the end of B first,\r
+so the exit code is pleasantly common to both the normal and error cases.\r
+\r
+If B is smaller (function merge_hi, which is merge_lo's "mirror image"),\r
+much the same, except that we need to merge right to left, copying B into a\r
+temp array and starting the stores at the right end of where B used to live.\r
+\r
+A refinement: When we're about to merge adjacent runs A and B, we first do\r
+a form of binary search (more on that later) to see where B[0] should end up\r
+in A. Elements in A preceding that point are already in their final\r
+positions, effectively shrinking the size of A. Likewise we also search to\r
+see where A[-1] should end up in B, and elements of B after that point can\r
+also be ignored. This cuts the amount of temp memory needed by the same\r
+amount.\r
+\r
+These preliminary searches may not pay off, and can be expected *not* to\r
+repay their cost if the data is random. But they can win huge in all of\r
+time, copying, and memory savings when they do pay, so this is one of the\r
+"per-merge overheads" mentioned above that we're happy to endure because\r
+there is at most one very short run. It's generally true in this algorithm\r
+that we're willing to gamble a little to win a lot, even though the net\r
+expectation is negative for random data.\r
+\r
+\r
+Merge Algorithms\r
+----------------\r
+merge_lo() and merge_hi() are where the bulk of the time is spent. merge_lo\r
+deals with runs where A <= B, and merge_hi where A > B. They don't know\r
+whether the data is clustered or uniform, but a lovely thing about merging\r
+is that many kinds of clustering "reveal themselves" by how many times in a\r
+row the winning merge element comes from the same run. We'll only discuss\r
+merge_lo here; merge_hi is exactly analogous.\r
+\r
+Merging begins in the usual, obvious way, comparing the first element of A\r
+to the first of B, and moving B[0] to the merge area if it's less than A[0],\r
+else moving A[0] to the merge area. Call that the "one pair at a time"\r
+mode. The only twist here is keeping track of how many times in a row "the\r
+winner" comes from the same run.\r
+\r
+If that count reaches MIN_GALLOP, we switch to "galloping mode". Here\r
+we *search* B for where A[0] belongs, and move over all the B's before\r
+that point in one chunk to the merge area, then move A[0] to the merge\r
+area. Then we search A for where B[0] belongs, and similarly move a\r
+slice of A in one chunk. Then back to searching B for where A[0] belongs,\r
+etc. We stay in galloping mode until both searches find slices to copy\r
+less than MIN_GALLOP elements long, at which point we go back to one-pair-\r
+at-a-time mode.\r
+\r
+A refinement: The MergeState struct contains the value of min_gallop that\r
+controls when we enter galloping mode, initialized to MIN_GALLOP.\r
+merge_lo() and merge_hi() adjust this higher when galloping isn't paying\r
+off, and lower when it is.\r
+\r
+\r
+Galloping\r
+---------\r
+Still without loss of generality, assume A is the shorter run. In galloping\r
+mode, we first look for A[0] in B. We do this via "galloping", comparing\r
+A[0] in turn to B[0], B[1], B[3], B[7], ..., B[2**j - 1], ..., until finding\r
+the k such that B[2**(k-1) - 1] < A[0] <= B[2**k - 1]. This takes at most\r
+roughly lg(B) comparisons, and, unlike a straight binary search, favors\r
+finding the right spot early in B (more on that later).\r
+\r
+After finding such a k, the region of uncertainty is reduced to 2**(k-1) - 1\r
+consecutive elements, and a straight binary search requires exactly k-1\r
+additional comparisons to nail it (see note REGION OF UNCERTAINTY). Then we\r
+copy all the B's up to that point in one chunk, and then copy A[0]. Note\r
+that no matter where A[0] belongs in B, the combination of galloping + binary\r
+search finds it in no more than about 2*lg(B) comparisons.\r
+\r
+If we did a straight binary search, we could find it in no more than\r
+ceiling(lg(B+1)) comparisons -- but straight binary search takes that many\r
+comparisons no matter where A[0] belongs. Straight binary search thus loses\r
+to galloping unless the run is quite long, and we simply can't guess\r
+whether it is in advance.\r
+\r
+If data is random and runs have the same length, A[0] belongs at B[0] half\r
+the time, at B[1] a quarter of the time, and so on: a consecutive winning\r
+sub-run in B of length k occurs with probability 1/2**(k+1). So long\r
+winning sub-runs are extremely unlikely in random data, and guessing that a\r
+winning sub-run is going to be long is a dangerous game.\r
+\r
+OTOH, if data is lopsided or lumpy or contains many duplicates, long\r
+stretches of winning sub-runs are very likely, and cutting the number of\r
+comparisons needed to find one from O(B) to O(log B) is a huge win.\r
+\r
+Galloping compromises by getting out fast if there isn't a long winning\r
+sub-run, yet finding such very efficiently when they exist.\r
+\r
+I first learned about the galloping strategy in a related context; see:\r
+\r
+ "Adaptive Set Intersections, Unions, and Differences" (2000)\r
+ Erik D. Demaine, Alejandro López-Ortiz, J. Ian Munro\r
+\r
+and its followup(s). An earlier paper called the same strategy\r
+"exponential search":\r
+\r
+ "Optimistic Sorting and Information Theoretic Complexity"\r
+ Peter McIlroy\r
+ SODA (Fourth Annual ACM-SIAM Symposium on Discrete Algorithms), pp\r
+ 467-474, Austin, Texas, 25-27 January 1993.\r
+\r
+and it probably dates back to an earlier paper by Bentley and Yao. The\r
+McIlroy paper in particular has good analysis of a mergesort that's\r
+probably strongly related to this one in its galloping strategy.\r
+\r
+\r
+Galloping with a Broken Leg\r
+---------------------------\r
+So why don't we always gallop? Because it can lose, on two counts:\r
+\r
+1. While we're willing to endure small per-merge overheads, per-comparison\r
+ overheads are a different story. Calling Yet Another Function per\r
+ comparison is expensive, and gallop_left() and gallop_right() are\r
+ too long-winded for sane inlining.\r
+\r
+2. Galloping can-- alas --require more comparisons than linear one-at-time\r
+ search, depending on the data.\r
+\r
+#2 requires details. If A[0] belongs before B[0], galloping requires 1\r
+compare to determine that, same as linear search, except it costs more\r
+to call the gallop function. If A[0] belongs right before B[1], galloping\r
+requires 2 compares, again same as linear search. On the third compare,\r
+galloping checks A[0] against B[3], and if it's <=, requires one more\r
+compare to determine whether A[0] belongs at B[2] or B[3]. That's a total\r
+of 4 compares, but if A[0] does belong at B[2], linear search would have\r
+discovered that in only 3 compares, and that's a huge loss! Really. It's\r
+an increase of 33% in the number of compares needed, and comparisons are\r
+expensive in Python.\r
+\r
+index in B where # compares linear # gallop # binary gallop\r
+A[0] belongs search needs compares compares total\r
+---------------- ----------------- -------- -------- ------\r
+ 0 1 1 0 1\r
+\r
+ 1 2 2 0 2\r
+\r
+ 2 3 3 1 4\r
+ 3 4 3 1 4\r
+\r
+ 4 5 4 2 6\r
+ 5 6 4 2 6\r
+ 6 7 4 2 6\r
+ 7 8 4 2 6\r
+\r
+ 8 9 5 3 8\r
+ 9 10 5 3 8\r
+ 10 11 5 3 8\r
+ 11 12 5 3 8\r
+ ...\r
+\r
+In general, if A[0] belongs at B[i], linear search requires i+1 comparisons\r
+to determine that, and galloping a total of 2*floor(lg(i))+2 comparisons.\r
+The advantage of galloping is unbounded as i grows, but it doesn't win at\r
+all until i=6. Before then, it loses twice (at i=2 and i=4), and ties\r
+at the other values. At and after i=6, galloping always wins.\r
+\r
+We can't guess in advance when it's going to win, though, so we do one pair\r
+at a time until the evidence seems strong that galloping may pay. MIN_GALLOP\r
+is 7, and that's pretty strong evidence. However, if the data is random, it\r
+simply will trigger galloping mode purely by luck every now and again, and\r
+it's quite likely to hit one of the losing cases next. On the other hand,\r
+in cases like ~sort, galloping always pays, and MIN_GALLOP is larger than it\r
+"should be" then. So the MergeState struct keeps a min_gallop variable\r
+that merge_lo and merge_hi adjust: the longer we stay in galloping mode,\r
+the smaller min_gallop gets, making it easier to transition back to\r
+galloping mode (if we ever leave it in the current merge, and at the\r
+start of the next merge). But whenever the gallop loop doesn't pay,\r
+min_gallop is increased by one, making it harder to transition back\r
+to galloping mode (and again both within a merge and across merges). For\r
+random data, this all but eliminates the gallop penalty: min_gallop grows\r
+large enough that we almost never get into galloping mode. And for cases\r
+like ~sort, min_gallop can fall to as low as 1. This seems to work well,\r
+but in all it's a minor improvement over using a fixed MIN_GALLOP value.\r
+\r
+\r
+Galloping Complication\r
+----------------------\r
+The description above was for merge_lo. merge_hi has to merge "from the\r
+other end", and really needs to gallop starting at the last element in a run\r
+instead of the first. Galloping from the first still works, but does more\r
+comparisons than it should (this is significant -- I timed it both ways). For\r
+this reason, the gallop_left() and gallop_right() (see note LEFT OR RIGHT)\r
+functions have a "hint" argument, which is the index at which galloping\r
+should begin. So galloping can actually start at any index, and proceed at\r
+offsets of 1, 3, 7, 15, ... or -1, -3, -7, -15, ... from the starting index.\r
+\r
+In the code as I type it's always called with either 0 or n-1 (where n is\r
+the # of elements in a run). It's tempting to try to do something fancier,\r
+melding galloping with some form of interpolation search; for example, if\r
+we're merging a run of length 1 with a run of length 10000, index 5000 is\r
+probably a better guess at the final result than either 0 or 9999. But\r
+it's unclear how to generalize that intuition usefully, and merging of\r
+wildly unbalanced runs already enjoys excellent performance.\r
+\r
+~sort is a good example of when balanced runs could benefit from a better\r
+hint value: to the extent possible, this would like to use a starting\r
+offset equal to the previous value of acount/bcount. Doing so saves about\r
+10% of the compares in ~sort. However, doing so is also a mixed bag,\r
+hurting other cases.\r
+\r
+\r
+Comparing Average # of Compares on Random Arrays\r
+------------------------------------------------\r
+[NOTE: This was done when the new algorithm used about 0.1% more compares\r
+ on random data than does its current incarnation.]\r
+\r
+Here list.sort() is samplesort, and list.msort() this sort:\r
+\r
+"""\r
+import random\r
+from time import clock as now\r
+\r
+def fill(n):\r
+ from random import random\r
+ return [random() for i in xrange(n)]\r
+\r
+def mycmp(x, y):\r
+ global ncmp\r
+ ncmp += 1\r
+ return cmp(x, y)\r
+\r
+def timeit(values, method):\r
+ global ncmp\r
+ X = values[:]\r
+ bound = getattr(X, method)\r
+ ncmp = 0\r
+ t1 = now()\r
+ bound(mycmp)\r
+ t2 = now()\r
+ return t2-t1, ncmp\r
+\r
+format = "%5s %9.2f %11d"\r
+f2 = "%5s %9.2f %11.2f"\r
+\r
+def drive():\r
+ count = sst = sscmp = mst = mscmp = nelts = 0\r
+ while True:\r
+ n = random.randrange(100000)\r
+ nelts += n\r
+ x = fill(n)\r
+\r
+ t, c = timeit(x, 'sort')\r
+ sst += t\r
+ sscmp += c\r
+\r
+ t, c = timeit(x, 'msort')\r
+ mst += t\r
+ mscmp += c\r
+\r
+ count += 1\r
+ if count % 10:\r
+ continue\r
+\r
+ print "count", count, "nelts", nelts\r
+ print format % ("sort", sst, sscmp)\r
+ print format % ("msort", mst, mscmp)\r
+ print f2 % ("", (sst-mst)*1e2/mst, (sscmp-mscmp)*1e2/mscmp)\r
+\r
+drive()\r
+"""\r
+\r
+I ran this on Windows and kept using the computer lightly while it was\r
+running. time.clock() is wall-clock time on Windows, with better than\r
+microsecond resolution. samplesort started with a 1.52% #-of-comparisons\r
+disadvantage, fell quickly to 1.48%, and then fluctuated within that small\r
+range. Here's the last chunk of output before I killed the job:\r
+\r
+count 2630 nelts 130906543\r
+ sort 6110.80 1937887573\r
+msort 6002.78 1909389381\r
+ 1.80 1.49\r
+\r
+We've done nearly 2 billion comparisons apiece at Python speed there, and\r
+that's enough <wink>.\r
+\r
+For random arrays of size 2 (yes, there are only 2 interesting ones),\r
+samplesort has a 50%(!) comparison disadvantage. This is a consequence of\r
+samplesort special-casing at most one ascending run at the start, then\r
+falling back to the general case if it doesn't find an ascending run\r
+immediately. The consequence is that it ends up using two compares to sort\r
+[2, 1]. Gratifyingly, timsort doesn't do any special-casing, so had to be\r
+taught how to deal with mixtures of ascending and descending runs\r
+efficiently in all cases.\r
+\r
+\r
+NOTES\r
+-----\r
+\r
+BINSORT\r
+A "binary insertion sort" is just like a textbook insertion sort, but instead\r
+of locating the correct position of the next item via linear (one at a time)\r
+search, an equivalent to Python's bisect.bisect_right is used to find the\r
+correct position in logarithmic time. Most texts don't mention this\r
+variation, and those that do usually say it's not worth the bother: insertion\r
+sort remains quadratic (expected and worst cases) either way. Speeding the\r
+search doesn't reduce the quadratic data movement costs.\r
+\r
+But in CPython's case, comparisons are extraordinarily expensive compared to\r
+moving data, and the details matter. Moving objects is just copying\r
+pointers. Comparisons can be arbitrarily expensive (can invoke arbitary\r
+user-supplied Python code), but even in simple cases (like 3 < 4) _all_\r
+decisions are made at runtime: what's the type of the left comparand? the\r
+type of the right? do they need to be coerced to a common type? where's the\r
+code to compare these types? And so on. Even the simplest Python comparison\r
+triggers a large pile of C-level pointer dereferences, conditionals, and\r
+function calls.\r
+\r
+So cutting the number of compares is almost always measurably helpful in\r
+CPython, and the savings swamp the quadratic-time data movement costs for\r
+reasonable minrun values.\r
+\r
+\r
+LEFT OR RIGHT\r
+gallop_left() and gallop_right() are akin to the Python bisect module's\r
+bisect_left() and bisect_right(): they're the same unless the slice they're\r
+searching contains a (at least one) value equal to the value being searched\r
+for. In that case, gallop_left() returns the position immediately before the\r
+leftmost equal value, and gallop_right() the position immediately after the\r
+rightmost equal value. The distinction is needed to preserve stability. In\r
+general, when merging adjacent runs A and B, gallop_left is used to search\r
+thru B for where an element from A belongs, and gallop_right to search thru A\r
+for where an element from B belongs.\r
+\r
+\r
+REGION OF UNCERTAINTY\r
+Two kinds of confusion seem to be common about the claim that after finding\r
+a k such that\r
+\r
+ B[2**(k-1) - 1] < A[0] <= B[2**k - 1]\r
+\r
+then a binary search requires exactly k-1 tries to find A[0]'s proper\r
+location. For concreteness, say k=3, so B[3] < A[0] <= B[7].\r
+\r
+The first confusion takes the form "OK, then the region of uncertainty is at\r
+indices 3, 4, 5, 6 and 7: that's 5 elements, not the claimed 2**(k-1) - 1 =\r
+3"; or the region is viewed as a Python slice and the objection is "but that's\r
+the slice B[3:7], so has 7-3 = 4 elements". Resolution: we've already\r
+compared A[0] against B[3] and against B[7], so A[0]'s correct location is\r
+already known wrt _both_ endpoints. What remains is to find A[0]'s correct\r
+location wrt B[4], B[5] and B[6], which spans 3 elements. Or in general, the\r
+slice (leaving off both endpoints) (2**(k-1)-1)+1 through (2**k-1)-1\r
+inclusive = 2**(k-1) through (2**k-1)-1 inclusive, which has\r
+ (2**k-1)-1 - 2**(k-1) + 1 =\r
+ 2**k-1 - 2**(k-1) =\r
+ 2*2**k-1 - 2**(k-1) =\r
+ (2-1)*2**(k-1) - 1 =\r
+ 2**(k-1) - 1\r
+elements.\r
+\r
+The second confusion: "k-1 = 2 binary searches can find the correct location\r
+among 2**(k-1) = 4 elements, but you're only applying it to 3 elements: we\r
+could make this more efficient by arranging for the region of uncertainty to\r
+span 2**(k-1) elements." Resolution: that confuses "elements" with\r
+"locations". In a slice with N elements, there are N+1 _locations_. In the\r
+example, with the region of uncertainty B[4], B[5], B[6], there are 4\r
+locations: before B[4], between B[4] and B[5], between B[5] and B[6], and\r
+after B[6]. In general, across 2**(k-1)-1 elements, there are 2**(k-1)\r
+locations. That's why k-1 binary searches are necessary and sufficient.\r
projects / mirror_edk2.git / blobdiff