diff --git a/src/Numerics.Tests/SortingTests.cs b/src/Numerics.Tests/SortingTests.cs
index c60c1136..4e8f4995 100644
--- a/src/Numerics.Tests/SortingTests.cs
+++ b/src/Numerics.Tests/SortingTests.cs
@@ -175,7 +175,8 @@ namespace MathNet.Numerics.UnitTests
             }
 
             // default sorting (Ascending)
-            list.Sort();
+            //list.Sort();
+            Sorting.Sort(list);
 
             // just check that the order is as expected, not that the items are correct
             for (var i = 1; i < list.Count; i++)
diff --git a/src/Numerics/Sorting.cs b/src/Numerics/Sorting.cs
index 13475410..03e68e14 100644
--- a/src/Numerics/Sorting.cs
+++ b/src/Numerics/Sorting.cs
@@ -38,7 +38,7 @@ namespace MathNet.Numerics
     public static class Sorting
     {
         /// <summary>
-        /// Sort a list of keys, in place using the quick sort algorithm using the quick sort algorithm.
+        /// Sort a list of keys, in place using the introsort algorithm.
         /// </summary>
         /// <typeparam name="T">The type of elements in the key list.</typeparam>
         /// <param name="keys">List to sort.</param>
@@ -46,42 +46,12 @@ namespace MathNet.Numerics
         public static void Sort<T>(IList<T> keys, IComparer<T> comparer = null)
         {
             int count = keys.Count;
-            if (count <= 1)
-            {
-                return;
-            }
 
             if (null == comparer)
             {
                 comparer = Comparer<T>.Default;
             }
 
-            if (count == 2)
-            {
-                if (comparer.Compare(keys[0], keys[1]) > 0)
-                {
-                    Swap(keys, 0, 1);
-                }
-                return;
-            }
-
-            // insertion sort
-            if (count <= 10)
-            {
-                for (int i = 1; i < count; i++)
-                {
-                    var key = keys[i];
-                    int j = i - 1;
-                    while (j >= 0 && comparer.Compare(keys[j], key) > 0)
-                    {
-                        keys[j + 1] = keys[j];
-                        j--;
-                    }
-                    keys[j + 1] = key;
-                }
-                return;
-            }
-
             // array case
             if (keys is T[] keysArray)
             {
@@ -97,11 +67,11 @@ namespace MathNet.Numerics
             }
 
             // local sort implementation
-            QuickSort(keys, comparer, 0, count - 1);
+            IntroSort<T>(keys, comparer, 0, count - 1, (int a, int b) => Swap(keys, a, b));
         }
 
         /// <summary>
-        /// Sort a list of keys and items with respect to the keys, in place using the quick sort algorithm.
+        /// Sort a list of keys and items with respect to the keys, in place using the introsort algorithm.
         /// </summary>
         /// <typeparam name="TKey">The type of elements in the key list.</typeparam>
         /// <typeparam name="TItem">The type of elements in the item list.</typeparam>
@@ -111,46 +81,12 @@ namespace MathNet.Numerics
         public static void Sort<TKey, TItem>(IList<TKey> keys, IList<TItem> items, IComparer<TKey> comparer = null)
         {
             int count = keys.Count;
-            if (count <= 1)
-            {
-                return;
-            }
 
             if (null == comparer)
             {
                 comparer = Comparer<TKey>.Default;
             }
 
-            if (count == 2)
-            {
-                if (comparer.Compare(keys[0], keys[1]) > 0)
-                {
-                    Swap(keys, 0, 1);
-                    Swap(items, 0, 1);
-                }
-                return;
-            }
-
-            // insertion sort
-            if (count <= 10)
-            {
-                for (int i = 1; i < count; i++)
-                {
-                    var key = keys[i];
-                    var item = items[i];
-                    int j = i - 1;
-                    while (j >= 0 && comparer.Compare(keys[j], key) > 0)
-                    {
-                        keys[j + 1] = keys[j];
-                        items[j + 1] = items[j];
-                        j--;
-                    }
-                    keys[j + 1] = key;
-                    items[j + 1] = item;
-                }
-                return;
-            }
-
             // array case
             if (keys is TKey[] keysArray && items is TItem[] itemsArray)
             {
@@ -159,11 +95,16 @@ namespace MathNet.Numerics
             }
 
             // local sort implementation
-            QuickSort(keys, items, comparer, 0, count - 1);
+            Action<int, int> swap = (int a, int b) =>
+            {
+                Swap<TKey>(keys, a, b);
+                Swap<TItem>(items, a, b);
+            };
+            IntroSort(keys, comparer, 0, count - 1, swap);
         }
 
         /// <summary>
-        /// Sort a list of keys, items1 and items2 with respect to the keys, in place using the quick sort algorithm.
+        /// Sort a list of keys, items1 and items2 with respect to the keys, in place using the introsort algorithm.
         /// </summary>
         /// <typeparam name="TKey">The type of elements in the key list.</typeparam>
         /// <typeparam name="TItem1">The type of elements in the first item list.</typeparam>
@@ -175,52 +116,20 @@ namespace MathNet.Numerics
         public static void Sort<TKey, TItem1, TItem2>(IList<TKey> keys, IList<TItem1> items1, IList<TItem2> items2, IComparer<TKey> comparer = null)
         {
             int count = keys.Count;
-            if (count <= 1)
-            {
-                return;
-            }
 
             if (null == comparer)
             {
                 comparer = Comparer<TKey>.Default;
             }
 
-            if (count == 2)
-            {
-                if (comparer.Compare(keys[0], keys[1]) > 0)
-                {
-                    Swap(keys, 0, 1);
-                    Swap(items1, 0, 1);
-                    Swap(items2, 0, 1);
-                }
-                return;
-            }
-
-            // insertion sort
-            if (count <= 10)
-            {
-                for (int i = 1; i < count; i++)
-                {
-                    var key = keys[i];
-                    var item1 = items1[i];
-                    var item2 = items2[i];
-                    int j = i - 1;
-                    while (j >= 0 && comparer.Compare(keys[j], key) > 0)
-                    {
-                        keys[j + 1] = keys[j];
-                        items1[j + 1] = items1[j];
-                        items2[j + 1] = items2[j];
-                        j--;
-                    }
-                    keys[j + 1] = key;
-                    items1[j + 1] = item1;
-                    items2[j + 1] = item2;
-                }
-                return;
-            }
-
             // local sort implementation
-            QuickSort(keys, items1, items2, comparer, 0, count - 1);
+            Action<int, int> swap = (int a, int b) =>
+            {
+                Swap<TKey>(keys, a, b);
+                Swap<TItem1>(items1, a, b);
+                Swap<TItem2>(items2, a, b);
+            };
+            IntroSort(keys, comparer, 0, count - 1, swap);
         }
 
         /// <summary>
@@ -470,6 +379,130 @@ namespace MathNet.Numerics
             QuickSortAll(primary, secondary, primaryComparer, secondaryComparer, 0, primary.Count - 1);
         }
 
+        /// <summary>
+        /// Recursively reorders the given list to satisfy the max heap property.
+        /// </summary>
+        /// <typeparam name="TKey">The type of elements in the key list.</typeparam>
+        /// <param name="keys">The list which is turned into a heap.</param>
+        /// <param name="comparer">The method with which to compare two elements of the heap.</param>
+        /// <param name="i">The index of the heap to heapify.</param>
+        /// <param name="left">The left boundary of the heapify.</param>
+        /// <param name="right">The right boundary of the heapify.</param>
+        /// <param name="swapper">An Action which takes the given indexes and swaps the key and satellite data.</param>
+        static void MaxHeapify<TKey>(IList<TKey> keys, IComparer<TKey> comparer, int i, int left, int right, Action<int, int> swapper)
+        {
+            i -= left;
+            int leftChild = 2 * i + 1;
+            int rightChild = 2 * i + 2;
+
+            int largest = i;
+
+            if (leftChild + left >= left && leftChild + left <= right && comparer.Compare(keys[leftChild + left], keys[largest + left]) > 0)
+            {
+                largest = leftChild;
+            }
+
+            if (rightChild + left >= left && rightChild + left <= right && comparer.Compare(keys[rightChild + left], keys[largest + left]) > 0)
+            {
+                largest = rightChild;
+            }
+
+            if (largest != i)
+            {
+                swapper(largest + left, i + left);
+                MaxHeapify(keys, comparer, largest + left, left, right, swapper);
+            }
+        }
+
+        /// <summary>
+        /// 
+        /// </summary>
+        /// <typeparam name="TKey">The type of elements in the key list.</typeparam>
+        /// <param name="keys">The list which is turned into a heap.</param>
+        /// <param name="comparer">The method with which to compare two elements of the heap.</param>
+        /// <param name="left">The left boundary of the heap.</param>
+        /// <param name="right">The right boundary of the heap.</param>
+        /// <param name="swapper">An Action which takes the given indexes and swaps the key and satellite data.</param>
+
+        static void BuildMaxHeap<TKey>(IList<TKey> keys, IComparer<TKey> comparer, int left, int right, Action<int, int> swapper)
+        {
+            for (int i = (right - left) / 2; i >= 0; i--)
+            {
+                MaxHeapify(keys, comparer, i + left, left, right, swapper);
+            }
+        }
+
+        /// <summary>
+        /// Recursive implementation for an in place introspective sort on a list.
+        /// </summary>
+        /// <typeparam name="TKey">The type of elements in the key list.</typeparam>
+        /// <param name="keys">The list which is sorted using intro sort.</param>
+        /// <param name="comparer">The method with which to compare two elements of the intro sort.</param>
+        /// <param name="left">The left boundary of the intro sort.</param>
+        /// <param name="right">The right boundary of the intro sort.</param>
+        /// <param name="swapper">An Action which takes the given indexes and swaps the key and satellite data.</param>
+        /// <param name="recursions">Tracks the number of recursions entered.</param>
+        static void IntroSort<TKey>(IList<TKey> keys, IComparer<TKey> comparer, int left, int right, Action<int, int> swapper, int recursions = 0)
+        {
+            const double ln2 = 0.69314718056; // Natural Logarithm of 2
+            if (left >= right)
+            {
+                return;
+            }
+
+            Random.CryptoRandomSource rand = new MathNet.Numerics.Random.CryptoRandomSource();
+            double max_recursion_depth = 2 * Math.Log(keys.Count) / ln2; // This is the cap on recursion depth used by the GNU STL
+
+            if (right - left < 16) // Insertion Sort is faster on very small sequences, 16 is the number that Array.Sort uses
+            {
+                for (int i = left + 1; i <= right; i++)
+                {
+                    int j = i;
+                    while (j > 0 && comparer.Compare(keys[j - 1], keys[j]) > 0)
+                    {
+                        swapper(j, j - 1);
+                        j--;
+                    }
+                }
+            }
+            else if (recursions > max_recursion_depth) // Heapsort is guaranteed O(n log n) 
+            {
+                BuildMaxHeap(keys, comparer, left, right, swapper);
+
+                int heapBoundary = right;
+
+                for (int i = right - left; i > 0; i--)
+                {
+                    swapper(left, left + i);
+                    MaxHeapify(keys, comparer, left, left, --heapBoundary, swapper);
+                }
+            }
+            else // Quicksort
+            {
+                int pivot_index = (rand.Next() % (right - left)) + left; // Don't need to worry about negatives because rand.Next() returns a non-negative number
+                swapper(pivot_index, right);
+
+                TKey pivot = keys[right];
+
+                int i = left - 1;
+                for (int j = left; j < right; j++)
+                {
+                    if (comparer.Compare(keys[j], pivot) <= 0)
+                    {
+                        i++;
+                        swapper(i, j);
+                    }
+                }
+
+                i++;
+                swapper(i, right);
+
+                recursions++;
+                IntroSort(keys, comparer, left, i - 1, swapper, recursions);
+                IntroSort(keys, comparer, i + 1, right, swapper, recursions);
+            }
+
+        }
 
         /// <summary>
         /// Recursive implementation for an in place quick sort on a list.