CF-368.18.tar.gz

[apple/cf.git] / Base.subproj / CFSortFunctions.c

/*
 * Copyright (c) 2005 Apple Computer, Inc. All rights reserved.
 *
 * @APPLE_LICENSE_HEADER_START@
 * 
 * This file contains Original Code and/or Modifications of Original Code
 * as defined in and that are subject to the Apple Public Source License
 * Version 2.0 (the 'License'). You may not use this file except in
 * compliance with the License. Please obtain a copy of the License at
 * http://www.opensource.apple.com/apsl/ and read it before using this
 * file.
 * 
 * The Original Code and all software distributed under the License are
 * distributed on an 'AS IS' basis, WITHOUT WARRANTY OF ANY KIND, EITHER
 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT OR NON-INFRINGEMENT.
 * Please see the License for the specific language governing rights and
 * limitations under the License.
 * 
 * @APPLE_LICENSE_HEADER_END@
 */
/*      CFSortFunctions.c
        Copyright 1999-2002, Apple, Inc. All rights reserved.
        Responsibility: Christopher Kane
*/

/* This file contains code copied from the Libc project's files
   qsort.c, merge.c, and heapsort.c, and modified to suit the
   needs of CF, which needs the comparison callback to have an
   additional parameter. The code is collected into this one
   file so that the individual BSD sort routines can remain
   private and static.  We may not keep the bsd functions
   separate eventually.
*/

#include <CoreFoundation/CFBase.h>
#if defined(__MACH__) || defined(__LINUX__) || defined(__FREEBSD__)
#include <sys/types.h>
#else
#define EINVAL          22
#endif
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <stddef.h>
#include "CFInternal.h"

static int bsd_mergesort(void *base, int nmemb, int size, CFComparatorFunction cmp, void *context);

/* Comparator is passed the address of the values. */
void CFMergeSortArray(void *list, CFIndex count, CFIndex elementSize, CFComparatorFunction comparator, void *context) {
    bsd_mergesort(list, count, elementSize, comparator, context);
} 

#if 0
// non-recursing qsort with side index list which is the actual thing sorted
// XXX kept here for possible later reference
void cjk_big_isort(char *list, int size, int *idx_list, int (*comp)(const void *, const void *, void *), void *context, int lo, int hi) {
    int t, idx, idx2;
    char *v1;
    for (idx = lo + 1; idx < hi + 1; idx++) {
        t = idx_list[idx];
        v1 = list + t * size;
        for (idx2 = idx; lo < idx2 && (comp(v1, list + idx_list[idx2 - 1] * size, context) < 0); idx2--)
            idx_list[idx2] = idx_list[idx2 - 1];
        idx_list[idx2] = t;
    }
}

void cjk_big_qsort(char *list, int size, int *idx_list, int (*comp)(const void *, const void *, void *), void *context, int lo, int hi, int *recurse) {
    int p, idx, mid, g1, g2, g3, r1, r2;
    char *v1, *v2, *v3;
    g1 = lo;
    g2 = (hi + lo) / 2; //random() % (hi - lo + 1) + lo;
    g3 = hi;
    v1 = list + idx_list[g1] * size;
    v2 = list + idx_list[g2] * size;
    v3 = list + idx_list[g3] * size;
    if (comp(v1, v2, context) < 0) {
        p = comp(v2, v3, context) < 0 ? g2 : (comp(v1, v3, context) < 0 ? g3 : g1);
    } else {
        p = comp(v2, v3, context) > 0 ? g2 : (comp(v1, v3, context) < 0 ? g1 : g3);
    }
    if (lo != p) {int t = idx_list[lo]; idx_list[lo] = idx_list[p]; idx_list[p] = t;}
    r2 = 1;
    v2 = list + idx_list[lo] * size;
    for (mid = lo, idx = lo + 1; idx < hi + 1; idx++) {
        v1 = list + idx_list[idx] * size;
        r1 = comp(v1, v2, context);
        r2 = r2 && (0 == r1);
        if (r1 < 0) {
            mid++;
            if (idx != mid) {int t = idx_list[idx]; idx_list[idx] = idx_list[mid]; idx_list[mid] = t;}
        }
    }
    if (lo != mid) {int t = idx_list[lo]; idx_list[lo] = idx_list[mid]; idx_list[mid] = t;}
    *recurse = !r2 ? mid : -1;
}

void cjk_big_sort(char *list, int n, int size, int (*comp)(const void *, const void *, void *), void *context) {
    int len = 0, los[40], his[40];
    // 40 is big enough for 2^40 elements, in theory; in practice, much
    // lower but still sufficient; should recurse if the 40 limit is hit
    int *idx_list = malloc(sizeof(int) * n);
    char *tmp_list = malloc(n * size);
    int idx;
    for (idx = 0; idx < n; idx++) {
        idx_list[idx] = idx;
    }
    los[len] = 0;
    his[len++] = n - 1;
    while (0 < len) {
        int lo, hi;
        len--;
        lo = los[len];
        hi = his[len];
        if (5 < (hi - lo)) {
            int mid;
            cjk_big_qsort(list, size, idx_list, comp, context, lo, hi, &mid);
            if (0 < mid) {
                los[len] = lo;
                his[len++] = mid - 1;
                los[len] = mid + 1;
                his[len++] = hi;
            }
        } else {
            cjk_big_isort(list, size, idx_list, comp, context, lo, hi);
        }
    }
    for (idx = 0; idx < n; idx++)
        memmove(tmp_list + idx * size, list + idx_list[idx] * size, size);
    memmove(list, tmp_list, n * size);
    free(tmp_list);
    free(idx_list);
}
#endif


/* stdlib.subproj/qsort.c ============================================== */

/*
 * Copyright (c) 1992, 1993
 *      The Regents of the University of California.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *      This product includes software developed by the University of
 *      California, Berkeley and its contributors.
 * 4. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#if !defined(min)
#define min(a, b)       (a) < (b) ? a : b
#endif

/*
 * Qsort routine from Bentley & McIlroy's "Engineering a Sort Function".
 */
#define swapcode(TYPE, parmi, parmj, n) {               \
        long i = (n) / sizeof (TYPE);                   \
        register TYPE *pi = (TYPE *) (parmi);           \
        register TYPE *pj = (TYPE *) (parmj);           \
        do {                                            \
                register TYPE   t = *pi;                \
                *pi++ = *pj;                            \
                *pj++ = t;                              \
        } while (--i > 0);                              \
}

#define SWAPINIT(a, es) swaptype = ((char *)a - (char *)0) % sizeof(long) || \
        es % sizeof(long) ? 2 : es == sizeof(long)? 0 : 1;

CF_INLINE void
swapfunc(char *a, char *b, int n, int swaptype) {
        if(swaptype <= 1) 
                swapcode(long, a, b, n)
        else
                swapcode(char, a, b, n)
}

#define swap(a, b)                                      \
        if (swaptype == 0) {                            \
                long t = *(long *)(a);                  \
                *(long *)(a) = *(long *)(b);            \
                *(long *)(b) = t;                       \
        } else                                          \
                swapfunc(a, b, es, swaptype)

#define vecswap(a, b, n)        if ((n) > 0) swapfunc(a, b, n, swaptype)

CF_INLINE char *
med3(char *a, char *b, char *c, int (*cmp)(const void *, const void *, void *), void *context) {
        return cmp(a, b, context) < 0 ?
               (cmp(b, c, context) < 0 ? b : (cmp(a, c, context) < 0 ? c : a ))
              :(cmp(b, c, context) > 0 ? b : (cmp(a, c, context) < 0 ? a : c ));
}

/* Comparator is passed the address of the values. */
void CFQSortArray(void *aVoidStar, CFIndex n, CFIndex es, CFComparatorFunction cmp, void *context)
{
        char *a = (char *)aVoidStar;
        char *pa, *pb, *pc, *pd, *pl, *pm, *pn;
        int d, r, swaptype, swap_cnt;

loop:   SWAPINIT(a, es);
        swap_cnt = 0;
        if (n < 7) {
                for (pm = a + es; pm < (char *) a + n * es; pm += es)
                        for (pl = pm; pl > (char *) a && cmp(pl - es, pl, context) > 0;
                             pl -= es)
                                swap(pl, pl - es);
                return;
        }
        pm = a + (n / 2) * es;
        if (n > 7) {
                pl = a;
                pn = a + (n - 1) * es;
                if (n > 40) {
                        d = (n / 8) * es;
                        pl = med3(pl, pl + d, pl + 2 * d, cmp, context);
                        pm = med3(pm - d, pm, pm + d, cmp, context);
                        pn = med3(pn - 2 * d, pn - d, pn, cmp, context);
                }
                pm = med3(pl, pm, pn, cmp, context);
        }
        swap(a, pm);
        pa = pb = a + es;

        pc = pd = a + (n - 1) * es;
        for (;;) {
                while (pb <= pc && (r = cmp(pb, a, context)) <= 0) {
                        if (r == 0) {
                                swap_cnt = 1;
                                swap(pa, pb);
                                pa += es;
                        }
                        pb += es;
                }
                while (pb <= pc && (r = cmp(pc, a, context)) >= 0) {
                        if (r == 0) {
                                swap_cnt = 1;
                                swap(pc, pd);
                                pd -= es;
                        }
                        pc -= es;
                }
                if (pb > pc)
                        break;
                swap(pb, pc);
                swap_cnt = 1;
                pb += es;
                pc -= es;
        }
        if (swap_cnt == 0) {  /* Switch to insertion sort */
                for (pm = a + es; pm < (char *) a + n * es; pm += es)
                        for (pl = pm; pl > (char *) a && cmp(pl - es, pl, context) > 0; 
                             pl -= es)
                                swap(pl, pl - es);
                return;
        }

        pn = a + n * es;
        r = min(pa - (char *)a, pb - pa);
        vecswap(a, pb - r, r);
        r = min(pd - pc, pn - pd - es);
        vecswap(pb, pn - r, r);
        if ((r = pb - pa) > es)
                CFQSortArray(a, r / es, es, cmp, context);
        if ((r = pd - pc) > es) { 
                /* Iterate rather than recurse to save stack space */
                a = pn - r;
                n = r / es;
                goto loop;
        }
/*              CFQSortArray(pn - r, r / es, es, cmp, context);*/
}

#undef min
#undef swapcode
#undef SWAPINIT
#undef swap
#undef vecswap

/* stdlib.subproj/mergesort.c ========================================== */

/*
 * Copyright (c) 1992, 1993
 *      The Regents of the University of California.  All rights reserved.
 *
 * This code is derived from software contributed to Berkeley by
 * Peter McIlroy.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *      This product includes software developed by the University of
 *      California, Berkeley and its contributors.
 * 4. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */


/*
 * Hybrid exponential search/linear search merge sort with hybrid
 * natural/pairwise first pass.  Requires about .3% more comparisons
 * for random data than LSMS with pairwise first pass alone.
 * It works for objects as small as two bytes.
 */

#define NATURAL
#define THRESHOLD 16    /* Best choice for natural merge cut-off. */

/* #define NATURAL to get hybrid natural merge.
 * (The default is pairwise merging.)
 */


#define ISIZE sizeof(int)
#define PSIZE sizeof(uint8_t *)
#define ICOPY_LIST(src, dst, last)                              \
        do                                                      \
        *(int*)dst = *(int*)src, src += ISIZE, dst += ISIZE;    \
        while(src < last)
#define ICOPY_ELT(src, dst, i)                                  \
        do                                                      \
        *(int*) dst = *(int*) src, src += ISIZE, dst += ISIZE;  \
        while (i -= ISIZE)

#define CCOPY_LIST(src, dst, last)              \
        do                                      \
                *dst++ = *src++;                \
        while (src < last)
#define CCOPY_ELT(src, dst, i)                  \
        do                                      \
                *dst++ = *src++;                \
        while (i -= 1)
                
/*
 * Find the next possible pointer head.  (Trickery for forcing an array
 * to do double duty as a linked list when objects do not align with word
 * boundaries.
 */
/* Assumption: PSIZE is a power of 2. */
#define EVAL(p) (uint8_t **)                                            \
        ((uint8_t *)0 +                                                 \
            (((uint8_t *)p + PSIZE - 1 - (uint8_t *) 0) & ~(PSIZE - 1)))


#define swap(a, b) {                                    \
                s = b;                                  \
                i = size;                               \
                do {                                    \
                        tmp = *a; *a++ = *s; *s++ = tmp; \
                } while (--i);                          \
                a -= size;                              \
        }
#define reverse(bot, top) {                             \
        s = top;                                        \
        do {                                            \
                i = size;                               \
                do {                                    \
                        tmp = *bot; *bot++ = *s; *s++ = tmp; \
                } while (--i);                          \
                s -= size2;                             \
        } while(bot < s);                               \
}

/*
 * This is to avoid out-of-bounds addresses in sorting the
 * last 4 elements.
 */
static void
insertionsort(char *a, int n, int size, int (*cmp)(const void *, const void *, void *), void *context) {
        char *ai, *s, *t, *u, tmp;
        int i;

        for (ai = a+size; --n >= 1; ai += size)
                for (t = ai; t > a; t -= size) {
                        u = t - size;
                        if (cmp(u, t, context) <= 0)
                                break;
                        swap(u, t);
                }
}

/*
 * Optional hybrid natural/pairwise first pass.  Eats up list1 in runs of
 * increasing order, list2 in a corresponding linked list.  Checks for runs
 * when THRESHOLD/2 pairs compare with same sense.  (Only used when NATURAL
 * is defined.  Otherwise simple pairwise merging is used.)
 */
static void
setup(char *list1, char *list2, int n, int size, int (*cmp)(const void *, const void *, void *), void *context) {
        int i, length, size2, tmp, sense;
        char *f1, *f2, *s, *l2, *last, *p2;

        size2 = size*2;
        if (n <= 5) {
                insertionsort(list1, n, size, cmp, context);
                *EVAL(list2) = (uint8_t*) list2 + n*size;
                return;
        }
        /*
         * Avoid running pointers out of bounds; limit n to evens
         * for simplicity.
         */
        i = 4 + (n & 1);
        insertionsort(list1 + (n - i) * size, i, size, cmp, context);
        last = list1 + size * (n - i);
        *EVAL(list2 + (last - list1)) = list2 + n * size;

#ifdef NATURAL
        p2 = list2;
        f1 = list1;
        sense = (cmp(f1, f1 + size, context) > 0);
        for (; f1 < last; sense = !sense) {
                length = 2;
                                        /* Find pairs with same sense. */
                for (f2 = f1 + size2; f2 < last; f2 += size2) {
                        if ((cmp(f2, f2+ size, context) > 0) != sense)
                                break;
                        length += 2;
                }
                if (length < THRESHOLD) {               /* Pairwise merge */
                        do {
                                p2 = *EVAL(p2) = f1 + size2 - list1 + list2;
                                if (sense > 0)
                                        swap (f1, f1 + size);
                        } while ((f1 += size2) < f2);
                } else {                                /* Natural merge */
                        l2 = f2;
                        for (f2 = f1 + size2; f2 < l2; f2 += size2) {
                                if ((cmp(f2-size, f2, context) > 0) != sense) {
                                        p2 = *EVAL(p2) = f2 - list1 + list2;
                                        if (sense > 0)
                                                reverse(f1, f2-size);
                                        f1 = f2;
                                }
                        }
                        if (sense > 0)
                                reverse (f1, f2-size);
                        f1 = f2;
                        if (f2 < last || cmp(f2 - size, f2, context) > 0)
                                p2 = *EVAL(p2) = f2 - list1 + list2;
                        else
                                p2 = *EVAL(p2) = list2 + n*size;
                }
        }
#else           /* pairwise merge only. */
        for (f1 = list1, p2 = list2; f1 < last; f1 += size2) {
                p2 = *EVAL(p2) = p2 + size2;
                if (cmp (f1, f1 + size, context) > 0)
                        swap(f1, f1 + size);
        }
#endif /* NATURAL */
}

/*
 * Arguments are as for qsort.
 */
static int
bsd_mergesort(void *base, int nmemb, int size, CFComparatorFunction cmp, void *context) {
        register int i, sense;
        int big, iflag;
        register uint8_t *f1, *f2, *t, *b, *tp2, *q, *l1, *l2;
        uint8_t *list2, *list1, *p2, *p, *last, **p1;

        if (size < (int)PSIZE / 2) {            /* Pointers must fit into 2 * size. */
                errno = EINVAL;
                return (-1);
        }

        /*
         * XXX
         * Stupid subtraction for the Cray.
         */
        iflag = 0;
        if (!(size % ISIZE) && !(((char *)base - (char *)0) % ISIZE))
                iflag = 1;

        if ((list2 = malloc(nmemb * size + PSIZE)) == NULL)
                return (-1);

        list1 = base;
        setup(list1, list2, nmemb, size, cmp, context);
        last = list2 + nmemb * size;
        i = big = 0;
        while (*EVAL(list2) != last) {
            l2 = list1;
            p1 = EVAL(list1);
            for (tp2 = p2 = list2; p2 != last; p1 = EVAL(l2)) {
                p2 = *EVAL(p2);
                f1 = l2;
                f2 = l1 = list1 + (p2 - list2);
                if (p2 != last)
                        p2 = *EVAL(p2);
                l2 = list1 + (p2 - list2);
                while (f1 < l1 && f2 < l2) {
                        if ((*cmp)(f1, f2, context) <= 0) {
                                q = f2;
                                b = f1, t = l1;
                                sense = -1;
                        } else {
                                q = f1;
                                b = f2, t = l2;
                                sense = 0;
                        }
                        if (!big) {     /* here i = 0 */
/* LINEAR: */                           while ((b += size) < t && cmp(q, b, context) >sense)
                                        if (++i == 6) {
                                                big = 1;
                                                goto EXPONENTIAL;
                                        }
                        } else {
EXPONENTIAL:                    for (i = size; ; i <<= 1)
                                        if ((p = (b + i)) >= t) {
                                                if ((p = t - size) > b &&
                                                    (*cmp)(q, p, context) <= sense)
                                                        t = p;
                                                else
                                                        b = p;
                                                break;
                                        } else if ((*cmp)(q, p, context) <= sense) {
                                                t = p;
                                                if (i == size)
                                                        big = 0; 
                                                goto FASTCASE;
                                        } else
                                                b = p;
/* SLOWCASE: */                 while (t > b+size) {
                                        i = (((t - b) / size) >> 1) * size;
                                        if ((*cmp)(q, p = b + i, context) <= sense)
                                                t = p;
                                        else
                                                b = p;
                                }
                                goto COPY;
FASTCASE:                       while (i > size)
                                        if ((*cmp)(q,
                                                p = b + (i >>= 1), context) <= sense)
                                                t = p;
                                        else
                                                b = p;
COPY:                           b = t;
                        }
                        i = size;
                        if (q == f1) {
                                if (iflag) {
                                        ICOPY_LIST(f2, tp2, b);
                                        ICOPY_ELT(f1, tp2, i);
                                } else {
                                        CCOPY_LIST(f2, tp2, b);
                                        CCOPY_ELT(f1, tp2, i);
                                }
                        } else {
                                if (iflag) {
                                        ICOPY_LIST(f1, tp2, b);
                                        ICOPY_ELT(f2, tp2, i);
                                } else {
                                        CCOPY_LIST(f1, tp2, b);
                                        CCOPY_ELT(f2, tp2, i);
                                }
                        }
                }
                if (f2 < l2) {
                        if (iflag)
                                ICOPY_LIST(f2, tp2, l2);
                        else
                                CCOPY_LIST(f2, tp2, l2);
                } else if (f1 < l1) {
                        if (iflag)
                                ICOPY_LIST(f1, tp2, l1);
                        else
                                CCOPY_LIST(f1, tp2, l1);
                }
                *p1 = l2;
            }
            tp2 = list1;        /* swap list1, list2 */
            list1 = list2;
            list2 = tp2;
            last = list2 + nmemb*size;
        }
        if (base == list2) {
                memmove(list2, list1, nmemb*size);
                list2 = list1;
        }
        free(list2);
        return (0);
}

#undef NATURAL
#undef THRESHOLD
#undef ISIZE
#undef PSIZE
#undef ICOPY_LIST
#undef ICOPY_ELT
#undef CCOPY_LIST
#undef CCOPY_ELT
#undef EVAL
#undef swap
#undef reverse

#if 0
/* stdlib.subproj/heapsort.c =========================================== */

/*
 * Copyright (c) 1991, 1993
 *      The Regents of the University of California.  All rights reserved.
 *
 * This code is derived from software contributed to Berkeley by
 * Ronnie Kon at Mindcraft Inc., Kevin Lew and Elmer Yglesias.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *      This product includes software developed by the University of
 *      California, Berkeley and its contributors.
 * 4. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

/*
 * Swap two areas of size number of bytes.  Although qsort(3) permits random
 * blocks of memory to be sorted, sorting pointers is almost certainly the
 * common case (and, were it not, could easily be made so).  Regardless, it
 * isn't worth optimizing; the SWAP's get sped up by the cache, and pointer
 * arithmetic gets lost in the time required for comparison function calls.
 */
#define SWAP(a, b, count, size, tmp) { \
        count = size; \
        do { \
                tmp = *a; \
                *a++ = *b; \
                *b++ = tmp; \
        } while (--count); \
}

/* Copy one block of size size to another. */
#define COPY(a, b, count, size, tmp1, tmp2) { \
        count = size; \
        tmp1 = a; \
        tmp2 = b; \
        do { \
                *tmp1++ = *tmp2++; \
        } while (--count); \
}

/*
 * Build the list into a heap, where a heap is defined such that for
 * the records K1 ... KN, Kj/2 >= Kj for 1 <= j/2 <= j <= N.
 *
 * There two cases.  If j == nmemb, select largest of Ki and Kj.  If
 * j < nmemb, select largest of Ki, Kj and Kj+1.
 */
#define CREATE(initval, nmemb, par_i, child_i, par, child, size, count, tmp) { \
        for (par_i = initval; (child_i = par_i * 2) <= nmemb; \
            par_i = child_i) { \
                child = base + child_i * size; \
                if (child_i < nmemb && compar(child, child + size, context) < 0) { \
                        child += size; \
                        ++child_i; \
                } \
                par = base + par_i * size; \
                if (compar(child, par, context) <= 0) \
                        break; \
                SWAP(par, child, count, size, tmp); \
        } \
}

/*
 * Select the top of the heap and 'heapify'.  Since by far the most expensive
 * action is the call to the compar function, a considerable optimization
 * in the average case can be achieved due to the fact that k, the displaced
 * elememt, is ususally quite small, so it would be preferable to first
 * heapify, always maintaining the invariant that the larger child is copied
 * over its parent's record.
 *
 * Then, starting from the *bottom* of the heap, finding k's correct place,
 * again maintianing the invariant.  As a result of the invariant no element
 * is 'lost' when k is assigned its correct place in the heap.
 *
 * The time savings from this optimization are on the order of 15-20% for the
 * average case. See Knuth, Vol. 3, page 158, problem 18.
 *
 * XXX Don't break the #define SELECT line, below.  Reiser cpp gets upset.
 */
#define SELECT(par_i, child_i, nmemb, par, child, size, k, count, tmp1, tmp2) { \
        for (par_i = 1; (child_i = par_i * 2) <= nmemb; par_i = child_i) { \
                child = base + child_i * size; \
                if (child_i < nmemb && compar(child, child + size, context) < 0) { \
                        child += size; \
                        ++child_i; \
                } \
                par = base + par_i * size; \
                COPY(par, child, count, size, tmp1, tmp2); \
        } \
        for (;;) { \
                child_i = par_i; \
                par_i = child_i / 2; \
                child = base + child_i * size; \
                par = base + par_i * size; \
                if (child_i == 1 || compar(k, par, context) < 0) { \
                        COPY(child, k, count, size, tmp1, tmp2); \
                        break; \
                } \
                COPY(child, par, count, size, tmp1, tmp2); \
        } \
}

/*
 * Heapsort -- Knuth, Vol. 3, page 145.  Runs in O (N lg N), both average
 * and worst.  While heapsort is faster than the worst case of quicksort,
 * the BSD quicksort does median selection so that the chance of finding
 * a data set that will trigger the worst case is nonexistent.  Heapsort's
 * only advantage over quicksort is that it requires little additional memory.
 */
static int
bsd_heapsort(void *vbase, int nmemb, int size, int (*compar)(const void *, const void *, void *), void *context) {
        register int cnt, i, j, l;
        register char tmp, *tmp1, *tmp2;
        char *base, *k, *p, *t;

        if (nmemb <= 1)
                return (0);

        if (!size) {
                errno = EINVAL;
                return (-1);
        }

        if ((k = malloc(size)) == NULL)
                return (-1);

        /*
         * Items are numbered from 1 to nmemb, so offset from size bytes
         * below the starting address.
         */
        base = (char *)vbase - size;

        for (l = nmemb / 2 + 1; --l;)
                CREATE(l, nmemb, i, j, t, p, size, cnt, tmp);

        /*
         * For each element of the heap, save the largest element into its
         * final slot, save the displaced element (k), then recreate the
         * heap.
         */
        while (nmemb > 1) {
                COPY(k, base + nmemb * size, cnt, size, tmp1, tmp2);
                COPY(base + nmemb * size, base + size, cnt, size, tmp1, tmp2);
                --nmemb;
                SELECT(i, j, nmemb, t, p, size, k, cnt, tmp1, tmp2);
        }
        free(k);
        return (0);
}

#undef SWAP
#undef COPY
#undef CREATE
#undef SELECT

#endif

/* ===================================================================== */

#undef EINVAL