htswanalysis/MACS/lib/gsl/gsl-1.11/sort/sort.c

   1 /*
   2  * Implement Heap sort -- direct and indirect sorting
   3  * Based on descriptions in Sedgewick "Algorithms in C"
   4  *
   5  * Copyright (C) 1999  Thomas Walter
   6  *
   7  * 18 February 2000: Modified for GSL by Brian Gough
   8  *
   9  * This is free software; you can redistribute it and/or modify it
  10  * under the terms of the GNU General Public License as published by the
  11  * Free Software Foundation; either version 3, or (at your option) any
  12  * later version.
  13  *
  14  * This source is distributed in the hope that it will be useful, but WITHOUT
  15  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  17  * for more details.
  18  */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_heapsort.h>
  23
  24 static inline void swap (void *base, size_t size, size_t i, size_t j);
  25 static inline void downheap (void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare);
  26
  27 /* Inline swap function for moving objects around */
  28
  29 static inline void
  30 swap (void *base, size_t size, size_t i, size_t j)
  31 {
  32   register char *a = size * i + (char *) base;
  33   register char *b = size * j + (char *) base;
  34   register size_t s = size;
  35
  36   if (i == j)
  37     return;
  38
  39   do
  40     {
  41       char tmp = *a;
  42       *a++ = *b;
  43       *b++ = tmp;
  44     }
  45   while (--s > 0);
  46 }
  47
  48 #define CMP(data,size,j,k) (compare((char *)(data) + (size) * (j), (char *)(data) + (size) * (k)))
  49
  50 static inline void
  51 downheap (void *data, const size_t size, const size_t N, size_t k, gsl_comparison_fn_t compare)
  52 {
  53   while (k <= N / 2)
  54     {
  55       size_t j = 2 * k;
  56
  57       if (j < N && CMP (data, size, j, j + 1) < 0)
  58         {
  59           j++;
  60         }
  61
  62       if (CMP (data, size, k, j) < 0)
  63         {
  64           swap (data, size, j, k);
  65         }
  66       else
  67         {
  68           break;
  69         }
  70
  71       k = j;
  72     }
  73 }
  74
  75 void
  76 gsl_heapsort (void *data, size_t count, size_t size, gsl_comparison_fn_t compare)
  77 {
  78   /* Sort the array in ascending order. This is a true inplace
  79      algorithm with N log N operations. Worst case (an already sorted
  80      array) is something like 20% slower */
  81
  82   size_t N;
  83   size_t k;
  84
  85   if (count == 0)
  86     {
  87       return;                   /* No data to sort */
  88     }
  89
  90   /* We have n_data elements, last element is at 'n_data-1', first at
  91      '0' Set N to the last element number. */
  92
  93   N = count - 1;
  94
  95   k = N / 2;
  96   k++;                          /* Compensate the first use of 'k--' */
  97   do
  98     {
  99       k--;
 100       downheap (data, size, N, k, compare);
 101     }
 102   while (k > 0);
 103
 104   while (N > 0)
 105     {
 106       /* first swap the elements */
 107       swap (data, size, 0, N);
 108
 109       /* then process the heap */
 110       N--;
 111
 112       downheap (data, size, N, 0, compare);
 113     }
 114 }