1 files changed, 114 insertions, 0 deletions
diff --git a/gsl-1.9/sort/sort.c b/gsl-1.9/sort/sort.c
new file mode 100644
index 0000000..9e86ab0
--- /dev/null
+++ b/gsl-1.9/sort/sort.c
@@ -0,0 +1,114 @@
+/*
+ * Implement Heap sort -- direct and indirect sorting
+ * Based on descriptions in Sedgewick "Algorithms in C"
+ *
+ * Copyright (C) 1999  Thomas Walter
+ *
+ * 18 February 2000: Modified for GSL by Brian Gough
+ *
+ * This is free software; you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the
+ * Free Software Foundation; either version 2, or (at your option) any
+ * later version.
+ *
+ * This source is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ * for more details.
+ */
+
+#include <config.h>
+#include <stdlib.h>
+#include <gsl/gsl_heapsort.h>
+
+static inline void swap (void *base, size_t size, size_t i, size_t j);
+static inline void downheap (void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare);
+
+/* Inline swap function for moving objects around */
+
+static inline void
+swap (void *base, size_t size, size_t i, size_t j)
+{
+  register char *a = size * i + (char *) base;
+  register char *b = size * j + (char *) base;
+  register size_t s = size;
+
+  if (i == j)
+    return;
+
+  do
+    {
+      char tmp = *a;
+      *a++ = *b;
+      *b++ = tmp;
+    }
+  while (--s > 0);
+}
+
+#define CMP(data,size,j,k) (compare((char *)(data) + (size) * (j), (char *)(data) + (size) * (k)))
+
+static inline void
+downheap (void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare)
+{
+  while (k <= N / 2)
+    {
+      size_t j = 2 * k;
+
+      if (j < N && CMP (data, size, j, j + 1) < 0)
+        {
+          j++;
+        }
+
+      if (CMP (data, size, k, j) < 0)
+        {
+          swap (data, size, j, k);
+        }
+      else
+        {
+          break;
+        }
+
+      k = j;
+    }
+}
+
+void
+gsl_heapsort (void *data, size_t count, size_t size, gsl_comparison_fn_t compare)
+{
+  /* Sort the array in ascending order. This is a true inplace
+     algorithm with N log N operations. Worst case (an already sorted
+     array) is something like 20% slower */
+
+  size_t N;
+  size_t k;
+
+  if (count == 0)
+    {
+      return;                   /* No data to sort */
+    }
+
+  /* We have n_data elements, last element is at 'n_data-1', first at
+     '0' Set N to the last element number. */
+
+  N = count - 1;
+
+  k = N / 2;
+  k++;                          /* Compensate the first use of 'k--' */
+  do
+    {
+      k--;
+      downheap (data, size, N, k, compare);
+    }
+  while (k > 0);
+
+  while (N > 0)
+    {
+      /* first swap the elements */
+      swap (data, size, 0, N);
+
+      /* then process the heap */
+      N--;
+
+      downheap (data, size, N, 0, compare);
+    }
+}