/*

Tragically, I have bollixed it up badly. It's such a simple

algorithm, but somehow I've gone off, because it does not always

work. Any ideas where I went off?

*/

#include <stdlib.h>

#include <string.h>

#include <errno.h>

/*

We will swap sizeof(long) bytes at a time if possible (if the data is

long aligned),

else a char at a time.

*/

static void swapfunc(char *a, char *b, const size_t e_size, const

int swaptype)

{

if (swaptype <= 1) {

size_t i = (e_size) / sizeof(long);

long *pi = (long *) (void *) a;

long *pj = (long *) (void *) b;

do {

const long t = *pi;

*pi++ = *pj;

*pj++ = t;

} while (--i 0);

} else {

size_t i = (e_size);

char *pi = (char *) (void *) a;

char *pj = (char *) (void *) b;

do {

const char t = *pi;

*pi++ = *pj;

*pj++ = t;

} while (--i 0);

}

}

/*

Citation:

http://www.ics.uci.edu/~eppstein/161/960125.html
This beautiful algorithm:

heapselect(L,k)

{

heap H = heapify(L)

for (i = 1; i < k; i++) remove min(H)

return min(H)

}

has slowly transmongrified itself into this awful slop:

*/

int heapselect(void *voidbase, const size_t rank, const

size_t nmemb, const size_t e_size, int (*cmp) (const void *, const

void *))

{

char *base = voidbase; /* Base of the array */

char *tmp; /* will hold temporary element after

allocation */

size_t i, /* array position indicator */

i_stride, /* array position indicator, adjusted

for element stride */

j, /* array position indicator */

j_stride, /* array position indicator, adjusted

for element stride */

left, /* array position indicator */

mid, /* array position indicator */

mid_stride, /* array position indicator, adjusted

for element stride */

right; /* array position indicator */

/* swaptype will tell us about the alignment of the data type, for

opimization of exchange */

const unsigned swaptype = ((char *) (base) - (char *) 0) %

sizeof(long) || (e_size) %sizeof(long) ? 2 : (e_size) ==

sizeof(long) ? 0 : 1;

if (e_size == 0 || rank nmemb) {

return -EDOM;

}

if (nmemb <= 1) {

return 0;

}

left = (rank + 1) / 2;

right = rank - 1;

if (!(tmp = malloc(e_size))) {

return -ENOMEM;

}

while (0 < left) {

left--;

i = left;

i_stride = i * e_size;

memcpy(tmp, base + i_stride, e_size);

while ((j = 2 * i + 1) <= right) {

j_stride = j * e_size;

char *current = base + j_stride;

if (j < right && (*cmp) (current, current + e_size) < 0) {

j++;

j_stride += e_size;

}

if ((*cmp) (base + j_stride, tmp) < 0) {

break;

}

memcpy(base + i_stride, base + j_stride, e_size);

i = j;

i_stride = j_stride;

}

if (i != left) {

memcpy(base + i_stride, tmp, e_size);

}

}

for (mid = right + 1; mid < nmemb; mid++) {

mid_stride = mid * e_size;

if ((*cmp) (base + mid_stride, base) < 0) {

if (swaptype == 0) {

char *p = base + mid_stride;

const long t = *(long *) (void *) p;

*(long *) (void *) p = *(long *) (void *) base;

*(long *) (void *) base = t;

} else

swapfunc(base + mid_stride, base, e_size, swaptype);

i = 0;

i_stride = 0;

memcpy(tmp, base + i_stride, e_size);

while ((j = 2 * i + 1) <= right) {

j_stride = j * e_size;

char *current = base + j_stride;

if (j < right && (*cmp) (current, current + e_size) <

0) {

j++;

j_stride += e_size;

}

if ((*cmp) (base + j_stride, tmp) < 0) {

break;

}

memcpy(base + i_stride, base + j_stride, e_size);

i = j;

i_stride = j_stride;

}

if (i != 0) {

memcpy(base + i_stride, tmp, e_size);

}

}

}

free(tmp);

return 0;

}

#ifdef UNIT_TEST

#include <time.h>

#include <stdio.h>

int arr[100000];

int compare(const void *aa, const void *bb)

{

int a = *(int *) aa;

int b = *(int *) bb;

if (a b) return 1;

if (a < b) return -1;

return 0;

}

int main(void)

{

size_t index;

size_t rank;

size_t arrlen = sizeof arr / sizeof arr[0];

int position;

srand((unsigned) time(NULL));

for (index = 0; index < arrlen; index++) {

arr[index] = rand();

}

rank = rand() / (RAND_MAX / arrlen + 1);

position = heapselect(arr, rank, arrlen, sizeof arr[0], compare);

if (position < 0)

printf("Error in heapselect = %d\n", position);

else

printf("item (%u) is %d\n", (unsigned) rank, arr[position]);

qsort(arr, arrlen, sizeof arr[0], compare);

printf("sorted item (%u) is %d\n", (unsigned) rank, arr[rank]);

return 0;

}

#endif