fork download
copy 
  1.  /***************************************************************************
  2.  * File: Smoothsort.hh
  3.  * Author: Keith Schwarz (htiek@cs.stanford.edu)
  4.  *
  5.  * An implementation of Dijkstra's Smoothsort algorithm, a modification of
  6.  * heapsort that runs in O(n lg n) in the worst case, but O(n) if the data
  7.  * are already sorted.  For more information about how this algorithm works
  8.  * and some of the details necessary for its proper operation, please see
  9.  *
  10.  *              http://w...content-available-to-author-only...z.com/smoothsort/
  11.  *
  12.  * This implementation is designed to work on a 32-bit machine and may
  13.  * have portability issues on 64-bit computers.  In particular, I've only
  14.  * precomputed the Leonardo numbers up 2^32, and so if you try to sort a
  15.  * sequence of length greater than that you'll run into trouble.  Similarly,
  16.  * I've used the tricky O(1) optimization to use a constant amount of space
  17.  * given the fact that the machine is 32 bits.
  18.  */
  19. #ifndef Smoothsort_Included
  20. #define Smoothsort_Included
  21.  
  22. #include <iterator>
  23. #include <algorithm>
  24. #include <bitset>
  25.  
  26. /**
  27.  * Function: Smoothsort(RandomIterator begin, RandomIterator end);
  28.  * -----------------------------------------------------------------------
  29.  * Sorts the input range into ascending order using the smoothsort
  30.  * algorithm.
  31.  */
  32. template <typename RandomIterator>
  33. void Smoothsort(RandomIterator begin, RandomIterator end);
  34.  
  35. /**
  36.  * Function: Smoothsort(RandomIterator begin, RandomIterator end,
  37.  *                      Comparator comp);
  38.  * -----------------------------------------------------------------------
  39.  * Sorts the input range into ascending order according to the strict
  40.  * total ordering comp using the smoothsort algorithm.
  41.  */
  42. template <typename RandomIterator, typename Comparator>
  43. void Smoothsort(RandomIterator begin, RandomIterator end, Comparator comp);
  44.  
  45. /* * * * * Implementation Below This Point * * * * */
  46. namespace smoothsort_detail {
  47.   /* A constant containing the number of Leonardo numbers that can fit into
  48.    * 32 bits.  For a 64-bit machine, you'll need to update this value and the
  49.    * list below.
  50.    */
  51.   const size_t kNumLeonardoNumbers = 46;
  52.  
  53.   /* A list of all the Leonardo numbers below 2^32, precomputed for
  54.    * efficiency.
  55.    * Source: http://o...content-available-to-author-only...s.org/classic/b001595.txt
  56.    */
  57.   const size_t kLeonardoNumbers[kNumLeonardoNumbers] = {
  58.     1u, 1u, 3u, 5u, 9u, 15u, 25u, 41u, 67u, 109u, 177u, 287u, 465u, 753u, 
  59.     1219u, 1973u, 3193u, 5167u, 8361u, 13529u, 21891u, 35421u, 57313u, 92735u,
  60.     150049u, 242785u, 392835u, 635621u, 1028457u, 1664079u, 2692537u, 
  61.     4356617u, 7049155u, 11405773u, 18454929u, 29860703u, 48315633u, 78176337u,
  62.     126491971u, 204668309u, 331160281u, 535828591u, 866988873u, 1402817465u,
  63.     2269806339u, 3672623805u
  64.   };
  65.  
  66.   /* A structure containing a bitvector encoding of the trees in a Leonardo
  67.    * heap.  The representation is as a bitvector shifted down so that its
  68.    * first digit is a one, along with the amount that it was shifted.
  69.    */
  70.   struct HeapShape {
  71.     /* A bitvector capable of holding all the Leonardo numbers. */
  72.     std::bitset<kNumLeonardoNumbers> trees;
  73.  
  74.     /* The shift amount, which is also the size of the smallest tree. */
  75.     size_t smallestTreeSize;
  76.   };
  77.  
  78.   /**
  79.    * Function: RandomIterator SecondChild(RandomIterator root)
  80.    * ---------------------------------------------------------------------
  81.    * Given an iterator to the root of Leonardo heap, returns an iterator
  82.    * to the root of that tree's second child.  It's assumed that the heap
  83.    * is well-formed and that size > 1.
  84.    */
  85.   template <typename RandomIterator>
  86.   RandomIterator SecondChild(RandomIterator root) {
  87.     /* The second child root is always one step before the root. */
  88.     return root - 1;
  89.   }
  90.  
  91.   /**
  92.    * Function: RandomIterator FirstChild(RandomIterator root, size_t size)
  93.    * ---------------------------------------------------------------------
  94.    * Given an iterator to the root of Leonardo heap, returns an iterator
  95.    * to the root of that tree's first child.  It's assumed that the heap
  96.    * is well-formed and that size > 1.
  97.    */
  98.   template <typename RandomIterator>
  99.   RandomIterator FirstChild(RandomIterator root, size_t size) {
  100.     /* Go to the second child, then step backwards L(size - 2) steps to
  101.      * skip over it.
  102.      */
  103.     return SecondChild(root) - kLeonardoNumbers[size - 2];
  104.   }
  105.  
  106.   /**
  107.    * Function: RandomIterator LargerChild(RandomIterator root, size_t size,
  108.    *                                      Comparator comp);
  109.    * --------------------------------------------------------------------
  110.    * Given an iterator to the root of a max-heap Leonardo tree, returns
  111.    * an iterator to its larger child.  It's assumed that the heap is
  112.    * well-formatted and that the heap has order > 1.
  113.    */
  114.   template <typename RandomIterator, typename Comparator>
  115.   RandomIterator LargerChild(RandomIterator root, size_t size, Comparator comp) {
  116.     /* Get pointers to the first and second child. */
  117.     RandomIterator first  = FirstChild(root, size);
  118.     RandomIterator second = SecondChild(root);
  119.  
  120.     /* Determine which is greater. */
  121.     return comp(*first, *second)? second : first;
  122.   }
  123.  
  124.   /**
  125.    * Function: RebalanceSingleHeap(RandomIterator root, size_t size, 
  126.    *                               Comparator comp);
  127.    * --------------------------------------------------------------------
  128.    * Given an iterator to the root of a single Leonardo tree that needs
  129.    * rebalancing, rebalances that tree using the standard "bubble-down"
  130.    * approach.
  131.    */
  132.   template <typename RandomIterator, typename Comparator>
  133.   void RebalanceSingleHeap(RandomIterator root, size_t size, Comparator comp) {
  134.     /* Loop until the current node has no children, which happens when the order
  135.      * of the tree is 0 or 1.
  136.      */
  137.     while (size > 1) {
  138.       /* Get pointers to the first and second child. */
  139.       RandomIterator first  = FirstChild(root, size);
  140.       RandomIterator second = SecondChild(root);
  141.  
  142.       /* Determine which child is larger and remember the order of its tree. */
  143.       RandomIterator largerChild;
  144.       size_t childSize;
  145.       if (comp(*first, *second)) { 
  146.         largerChild = second; // Second child is larger...
  147.         childSize = size - 2; // ... and has order k - 2.
  148.       } else {
  149.         largerChild = first;  // First child is larger...
  150.         childSize = size - 1; // ... and has order k - 1.
  151.       }
  152.  
  153.       /* If the root is bigger than this child, we're done. */
  154.       if (!comp(*root, *largerChild))
  155.         return;
  156.  
  157.       /* Otherwise, swap down and update our order. */
  158.       std::iter_swap(root, largerChild);
  159.       root = largerChild;
  160.       size = childSize;
  161.     }
  162.   }
  163.  
  164.   /**
  165.    * Function: LeonardoHeapRectify(RandomIterator begin, RandomIterator end,
  166.    *                               HeapShape shape, Comparator comp);
  167.    * ---------------------------------------------------------------------
  168.    * Given an implicit Leonardo heap spanning [begin, end) that has just
  169.    * had an element inserted into it at the very end, along with the 
  170.    * size list for that heap, rectifies the heap structure by shuffling
  171.    * the new root down to the proper position and rebalancing the target
  172.    * heap.
  173.    */
  174.   template <typename RandomIterator, typename Comparator>
  175.   void LeonardoHeapRectify(RandomIterator begin, RandomIterator end,
  176.                            HeapShape shape, Comparator comp) {
  177.     /* Back up the end iterator one step to get to the root of the rightmost
  178.      * heap.
  179.      */
  180.     RandomIterator itr = end - 1;
  181.  
  182.     /* Keep track of the size of the last heap size that we visited.  We need
  183.      * this so that once we've positioned the new node atop the correct heap
  184.      * we remember how large it is.
  185.      */
  186.     size_t lastHeapSize;
  187.  
  188.     /* Starting at the last heap and working backward, check whether we need
  189.      * to swap the root of the current heap with the previous root.
  190.      */
  191.     while (true) {
  192.       /* Cache the size of the heap we're currently on top of. */
  193.       lastHeapSize = shape.smallestTreeSize;
  194.  
  195.       /* If this is the very first heap in the tree, we're done. */
  196.       if (size_t(std::distance(begin, itr)) == kLeonardoNumbers[lastHeapSize] - 1)
  197.         break;
  198.  
  199.       /* We want to swap the previous root with this one if it's strictly
  200.        * greater than both the root of this tree and both its children.
  201.        * In order to avoid weird edge cases when the current heap has
  202.        * size zero or size one, we'll compute what value will be compared
  203.        * against.
  204.        */
  205.       RandomIterator toCompare = itr;
  206.  
  207.       /* If we aren't an order-0 or order-1 tree, we have two children, and
  208.        * need to check which of the three values is largest.
  209.        */
  210.       if (shape.smallestTreeSize > 1) {
  211.         /* Get the largest child and see if we need to change what we're 
  212.          * comparing against.
  213.          */
  214.         RandomIterator largeChild = LargerChild(itr, shape.smallestTreeSize,
  215.                                                 comp);
  216.  
  217.         /* Update what element is being compared against. */
  218.         if (comp(*toCompare, *largeChild))
  219.           toCompare = largeChild;
  220.       }
  221.  
  222.       /* Get a pointer to the root of the second heap by backing up the size
  223.        * of this heap.
  224.        */
  225.       RandomIterator priorHeap = itr - kLeonardoNumbers[lastHeapSize];
  226.  
  227.       /* If we ran out of trees or the new tree root is less than the element
  228.        * we're comparing, we now have the new node at the top of the correct
  229.        * heap.
  230.        */
  231.       if (!comp(*toCompare, *priorHeap))
  232.         break;
  233.  
  234.       /* Otherwise, do the swap and adjust our location. */
  235.       std::iter_swap(itr, priorHeap);
  236.       itr = priorHeap;
  237.  
  238.       /* Scan down until we find the heap before this one.  We do this by 
  239.        * continously shifting down the tree bitvector and bumping up the size
  240.        * of the smallest tree until we hit a new tree.
  241.        */
  242.       do {
  243.         shape.trees >>= 1;
  244.         ++shape.smallestTreeSize;
  245.       } while (!shape.trees[0]);
  246.     }
  247.  
  248.     /* Finally, rebalance the current heap. */
  249.     RebalanceSingleHeap(itr, lastHeapSize, comp);
  250.   }
  251.  
  252.   /**
  253.    * Function: LeonardoHeapAdd(RandomIterator begin, RandomIterator end,
  254.    *                           RandomIterator heapEnd,
  255.    *                           HeapShape& shape, Comparator comp);
  256.    * ----------------------------------------------------------------------
  257.    * Given an implicit Leonardo heap spanning [begin, end) in a range spanned
  258.    * by [begin, heapEnd], along with the shape and a comparator, increases the
  259.    * size of that heap by one by inserting the element at *end.
  260.    */
  261.   template <typename RandomIterator, typename Comparator>
  262.   void LeonardoHeapAdd(RandomIterator begin, RandomIterator end,
  263.                        RandomIterator heapEnd,
  264.                        HeapShape& shape, Comparator comp) {
  265.     /* There are three cases to consider, which are analogous to the cases
  266.      * in the proof that it is possible to partition the input into heaps
  267.      * of decreasing size:
  268.      *
  269.      * Case 0: If there are no elements in the heap, add a tree of order 1.
  270.      * Case 1: If the last two heaps have sizes that differ by one, we
  271.      *         add the new element by merging the last two heaps.
  272.      * Case 2: Otherwise, if the last heap has Leonardo number 1, add
  273.      *         a singleton heap of Leonardo number 0.
  274.      * Case 3: Otherwise, add a singleton heap of Leonardo number 1.
  275.      */
  276.  
  277.     /* Case 0 represented by the first bit being a zero; it should always be
  278.      * one during normal operation.
  279.      */
  280.     if (!shape.trees[0]) {
  281.       shape.trees[0] = true;
  282.       shape.smallestTreeSize = 1;
  283.     }
  284.     /* Case 1 would be represented by the last two bits of the bitvector both
  285.      * being set.
  286.      */
  287.     else if (shape.trees[1] && shape.trees[0]) {
  288.       /* First, remove those two trees by shifting them off the bitvector. */
  289.       shape.trees >>= 2;
  290.  
  291.       /* Set the last bit of the bitvector; we just added a tree of this
  292.        * size.
  293.        */
  294.       shape.trees[0] = true;
  295.  
  296.       /* Finally, increase the size of the smallest tree by two, since the new
  297.        * Leonardo tree has order one greater than both of them.
  298.        */
  299.       shape.smallestTreeSize += 2;
  300.     } 
  301.     /* Case two is represented by the size of the smallest tree being 1. */
  302.     else if (shape.smallestTreeSize == 1) {
  303.       /* Shift the bits up one spot so that we have room for the zero bit. */
  304.       shape.trees <<= 1;
  305.       shape.smallestTreeSize = 0;
  306.  
  307.       /* Set the bit. */
  308.       shape.trees[0] = true;
  309.     } 
  310.     /* Case three is everything else. */
  311.     else {
  312.       /* We currently have a forest encoded with a format that looks like
  313.        * (W, n) for bitstring W and exponent n.  We want to convert this to
  314.        * (W00...01, 1) by shifting up n - 1 spaces, then setting the last bit.
  315.        */
  316.       shape.trees <<= shape.smallestTreeSize - 1;
  317.       shape.trees[0] = true;
  318.  
  319.       /* Set the smallest tree size to one, since that is the new smallest
  320.        * tree size.
  321.        */
  322.       shape.smallestTreeSize = 1;
  323.     }
  324.  
  325.     /* At this point, we've set up a new tree.  We need to see if this tree is
  326.      * at the final size it's going to take.  If so, we'll do a full rectify
  327.      * on it.  Otherwise, all we need to do is maintain the heap property.
  328.      */
  329.     bool isLast = false;
  330.     switch (shape.smallestTreeSize) {
  331.       /* If this last heap has order 0, then it's in its final position only
  332.        * if it's the very last element of the array.
  333.        */
  334.     case 0:
  335.       if (end + 1 == heapEnd)
  336.         isLast = true;
  337.       break;
  338.  
  339.       /* If this last heap has order 1, then it's in its final position if
  340.        * it's the last element, or it's the penultimate element and it's not
  341.        * about to be merged.  For simplicity
  342.        */
  343.     case 1:
  344.       if (end + 1 == heapEnd || (end + 2 == heapEnd && !shape.trees[1]))
  345.         isLast = true;
  346.       break;
  347.  
  348.       /* Otherwise, this heap is in its final position if there isn't enough
  349.        * room for the next Leonardo number and one extra element.
  350.        */
  351.     default:
  352.       if (size_t(std::distance(end + 1, heapEnd)) < kLeonardoNumbers[shape.smallestTreeSize - 1] + 1)
  353.         isLast = true;
  354.       break;
  355.     }
  356.  
  357.     /* If this isn't a final heap, then just rebalance the current heap. */
  358.     if (!isLast)
  359.       RebalanceSingleHeap(end, shape.smallestTreeSize, comp);
  360.     /* Otherwise do a full rectify to put this node in its place. */
  361.     else
  362.       LeonardoHeapRectify(begin, end + 1, shape, comp);
  363.   }
  364.   /**
  365.    * Function: LeonardoHeapRemove(RandomIterator begin, RandomIterator end,
  366.    *                              HeapShape& shape,  Comparator comp);
  367.    * ----------------------------------------------------------------------
  368.    * Given an implicit Leonardo heap spanning [begin, end), along with the
  369.    * size list and a comparator, dequeues the element at end - 1 and
  370.    * rebalances the heap.  Since the largest element of the heap is already
  371.    * at end, this essentially keeps the max element in its place and does
  372.    * a rebalance if necessary.
  373.    */
  374.   template <typename RandomIterator, typename Comparator>
  375.   void LeonardoHeapRemove(RandomIterator begin, RandomIterator end,
  376.                           HeapShape& shape, Comparator comp) {
  377.     /* There are two cases to consider:
  378.      *
  379.      * Case 1: The last heap is of order zero or one.  In this case,
  380.      *         removing it doesn't expose any new trees and we can just
  381.      *         drop it from the list of trees.
  382.      * Case 2: The last heap is of order two or greater.  In this case,
  383.      *         we exposed two new heaps, which may require rebalancing.
  384.      */
  385.  
  386.     /* Case 1. */
  387.     if (shape.smallestTreeSize <= 1) {
  388.       /* Keep scanning up the list looking for the next tree. */
  389.       do {
  390.         shape.trees >>= 1;
  391.         ++shape.smallestTreeSize;
  392.       } while (shape.trees.any() && !shape.trees[0]);
  393.       return;
  394.     }
  395.  
  396.     /* Break open the last heap to expose two subheaps of order
  397.      * k - 2 and k - 1.  This works by mapping the encoding (W1, n) to the
  398.      * encoding (W011, n - 2).
  399.      */
  400.     const size_t heapOrder = shape.smallestTreeSize;
  401.     shape.trees[0] = false;
  402.     shape.trees <<= 2;
  403.     shape.trees[1] = shape.trees[0] = true;
  404.     shape.smallestTreeSize -= 2;
  405.  
  406.     /* We now do the insertion-sort/rebalance operation on the larger exposed heap to
  407.      * put it in its proper place, then on the smaller of the two.  But first, we need
  408.      * to find where they are.  This can be done by just looking up the first and second
  409.      * children of the former root, which was at end - 1.
  410.      */
  411.     RandomIterator leftHeap  = FirstChild(end - 1, heapOrder);
  412.     RandomIterator rightHeap = SecondChild(end - 1);
  413.  
  414.     /* Rebalance the left heap.  For this step we'll pretend that there is
  415.      * one fewer heap than there actually is, since we're ignoring the
  416.      * rightmost heap.
  417.      */
  418.     HeapShape allButLast = shape;
  419.     ++allButLast.smallestTreeSize;
  420.     allButLast.trees >>= 1;
  421.  
  422.     /* We add one to the position of the left heap because the function
  423.      * assumes an exclusive range, while leftHeap is actually an iterator
  424.      * directly to where the root is.
  425.      */
  426.     LeonardoHeapRectify(begin, leftHeap + 1,  allButLast, comp);
  427.     LeonardoHeapRectify(begin, rightHeap + 1, shape, comp);
  428.   }
  429. }
  430.  
  431. /* Actual smoothsort implementation. */
  432. template <typename RandomIterator, typename Comparator>
  433. void Smoothsort(RandomIterator begin, RandomIterator end, Comparator comp) {
  434.   /* Edge case: Check that the range isn't empty or a singleton. */
  435.   if (begin == end || begin + 1 == end) return;
  436.  
  437.   /* Construct a shape object describing the empty heap. */
  438.   smoothsort_detail::HeapShape shape;
  439.   shape.smallestTreeSize = 0;
  440.  
  441.   /* Convert the input into an implicit Leonardo heap. */
  442.   for (RandomIterator itr = begin; itr != end; ++itr)
  443.     smoothsort_detail::LeonardoHeapAdd(begin, itr, end, shape, comp);
  444.  
  445.   /* Continuously dequeue from the implicit Leonardo heap until we've
  446.    * consumed all the elements.
  447.    */
  448.   for (RandomIterator itr = end; itr != begin; --itr)
  449.     smoothsort_detail::LeonardoHeapRemove(begin, itr, shape, comp);
  450. }
  451.  
  452. /* Non-comparator version just uses the default comparator. */
  453. template <typename RandomIterator>
  454. void Smoothsort(RandomIterator begin, RandomIterator end) {
  455.   Smoothsort(begin, end,
  456.              std::less<typename std::iterator_traits<RandomIterator>::value_type>());
  457. }
  458.  
  459. #endif
  460.  
  461. #include <iostream>
  462. #include <vector>
  463. using namespace std;
  464.  
  465. int main()
  466. {
  467. 	const int init_values[] = { 3, 2, 1 };
  468.     std::vector<int> values(init_values, init_values+3);
  469.     Smoothsort(values.begin(), values.end());
  470. 	return 0;
  471. }
  472.  
Success #stdin #stdout 0s 16064KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/j38MFs
language:
C++ (gcc 8.3)
created:
7 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!