Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
236 lines (192 sloc) 7.3 KB
/* Quicksort.
*
* Differs from stdlib's qsort() in that Easel provides a standardized
* and reentrant way to sort arbitrary data arrays/structures by
* sorting an array of unique integer identifiers.
*
* Contents:
* 1. Quicksort algorithm
* 2. Unit tests
* 3. Test driver
* 4. Example
*/
#include "esl_config.h"
#include "easel.h"
#include "esl_quicksort.h"
static void partition(const void *data, int (*comparison)(const void *data, int o1, int o2), int *ord, int lo, int hi);
/*****************************************************************
* 1. Quicksort
*****************************************************************/
/* Function: esl_quicksort()
* Synopsis: Quicksort algorithm.
* Incept: SRE, Wed Nov 11 16:34:05 2015
*
* Purpose: <data> points to an array or structure that consists of
* <n> elements that we can identify with unique integer
* identifiers 0..<n-1>. Sort these elements, using a
* <*comparison()> function that returns -1, 0, or 1 when
* element <o1> should be sorted before, equal to, or after
* element <o2>. Caller provides an allocated array
* <sorted_at> to hold the result, allocated for at least
* <n> integers. The result in <sorted_at> consists of the
* ordered array of identifiers. For example,
* <sorted_at[0]> is the index of the best (first) element
* in the original <data> array, and <sorted_at[n-1]> is
* the id of the worst (last) element; <data[sorted_at[0]]>
* and <data[sorted_at[n-1]]> are the best and worst
* elements themselves. The original <data> array is
* unaltered.
*
* Compared to standard <qsort()>, the main advantage of
* <esl_qsort()> is that it is reentrant: it passes an
* arbitrary <data> pointer to the comparison function, so
* it can sort the elements of arbitrary structures or
* arrays without having to globally or statically expose
* those data. (A reentrant <qsort_r()> function is
* available in the libc of some platforms, but is not
* standard POSIX ANSI C, so we cannot rely on it.) The
* main disadvantage is that <esl_quicksort()> takes
* additional memory (the <sorted_at> array), whereas
* <qsort()> sorts in place.
*
* Args: data - generic pointer to the data to be sorted.
* n - <data> consists of <n> elements numbered 0..n-1
* comparison() - returns -1,0,1 if o1 is better, equal, worse than o2
* sorted_at - sorted array of the data[] idx's 0..n-1.
*
* Returns: <eslOK> on success.
*/
int
esl_quicksort(const void *data, int n, int (*comparison)(const void *data, int o1, int o2), int *sorted_at)
{
int i;
for (i = 0; i < n; i++) sorted_at[i] = i;
partition(data, comparison, sorted_at, 0, n-1);
return eslOK;
}
/*
* We're sorting a subrange of <ord>, from <ord[lo]..ord[hi]>
* Values in <ord> are unique identifiers for <data>.
* We're sorting <ord> by data[ord[i]] better than data[ord[i+1]]
*/
static void
partition(const void *data, int (*comparison)(const void *data, int o1, int o2), int *ord, int lo, int hi)
{
int i = lo;
int j = hi+1;
int swap, pivot;
/* Select pivot position from lo, hi, lo + (hi-lo)/2 by median-of-three */
if (comparison(data, ord[hi], ord[lo]) < 0) { swap = ord[hi]; ord[hi] = ord[lo]; ord[hi] = swap; }
pivot = lo + (hi-lo)/2;
if (comparison(data, ord[pivot], ord[lo]) < 0) pivot = lo;
else if (comparison(data, ord[pivot], ord[hi]) > 0) pivot = hi;
swap = ord[pivot]; ord[pivot] = ord[lo]; ord[lo] = swap;
/* Partition */
while (1)
{
do { i++; } while (i <= hi && comparison(data, ord[i], ord[lo]) < 0);
do { j--; } while ( comparison(data, ord[j], ord[lo]) > 0);
if (j > i) { swap = ord[j]; ord[j] = ord[i]; ord[i] = swap; }
else break;
}
swap = ord[lo]; ord[lo] = ord[j]; ord[j] = swap;
/* Recurse, doing the smaller partition first */
if (j-lo < hi-j) {
if (j-lo > 1) partition(data, comparison, ord, lo, j-1);
if (hi-j > 1) partition(data, comparison, ord, j+1, hi);
} else {
if (hi-j > 1) partition(data, comparison, ord, j+1, hi);
if (j-lo > 1) partition(data, comparison, ord, lo, j-1);
}
}
/*****************************************************************
* 2. Unit tests
*****************************************************************/
#ifdef eslQUICKSORT_TESTDRIVE
#include "esl_random.h"
int
sort_floats_ascending(const void *data, int o1, int o2)
{
float *x = (float *) data;
if (x[o1] < x[o2]) return -1;
if (x[o1] > x[o2]) return 1;
return 0;
}
static void
utest_floatsort(ESL_RANDOMNESS *rng, int N, int K)
{
char msg[] = "esl_quicksort: float sort test failed";
float *x = malloc(sizeof(float) * N);
int *sorted_at = malloc(sizeof(int) * N);
int i,r;
for (i = 0; i < N; i++)
x[i] = esl_random(rng) * K;
esl_quicksort(x, N, sort_floats_ascending, sorted_at);
for (r = 1; r < N; r++)
if (x[sorted_at[r]] < x[sorted_at[r-1]]) esl_fatal(msg);
free(x);
free(sorted_at);
}
#endif /*eslQUICKSORT_TESTDRIVE*/
/*****************************************************************
* 3. Test driver
*****************************************************************/
#ifdef eslQUICKSORT_TESTDRIVE
#include "easel.h"
#include "esl_getopts.h"
#include "esl_random.h"
#include "esl_quicksort.h"
#include <stdio.h>
static ESL_OPTIONS options[] = {
/* name type default env range togs reqs incomp help docgrp */
{"-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show help and usage", 0},
{"-s", eslARG_INT, "0", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>", 0},
{ 0,0,0,0,0,0,0,0,0,0},
};
static char usage[] = "[-options]";
static char banner[] = "test driver for quicksort module";
int
main(int argc, char **argv)
{
ESL_GETOPTS *go = esl_getopts_CreateDefaultApp(options, 0, argc, argv, banner, usage);
ESL_RANDOMNESS *rng = esl_randomness_Create(esl_opt_GetInteger(go, "-s"));
int N = 100;
int K = 10;
utest_floatsort(rng, N, K);
esl_randomness_Destroy(rng);
esl_getopts_Destroy(go);
return 0;
}
#endif /*eslQUICKSORT_TESTDRIVE*/
/*****************************************************************
* 4. Example
*****************************************************************/
#ifdef eslQUICKSORT_EXAMPLE
#include "easel.h"
#include "esl_random.h"
#include "esl_quicksort.h"
int
sort_floats_descending(const void *data, int o1, int o2)
{
float *x = (float *) data;
if (x[o1] > x[o2]) return -1;
if (x[o1] < x[o2]) return 1;
return 0;
}
int
main(int argc, char **argv)
{
ESL_RANDOMNESS *rng = esl_randomness_Create(0);
int N = 100;
float K = 10.0;
float *x = malloc(sizeof(float) * N);
int *sorted_at = malloc(sizeof(int) * N);
int i,r;
for (i = 0; i < N; i++)
x[i] = esl_random(rng) * K;
esl_quicksort(x, N, sort_floats_descending, sorted_at);
for (r = 0; r < N; r++)
printf("%.4f\n", x[sorted_at[r]]);
return 0;
}
#endif /*eslQUICKSORT_EXAMPLE*/
You can’t perform that action at this time.