/* Version 0.0
** From "Wayne's Little DSA Library" (DSA == Data Structures and
** Algorithms) Feel free to change, modify, or burn these sources, but if
** you modify them please don't distribute your changes without a clear
** indication that you've done so.  If you think your change is spiffy,
** send it to me and maybe I'll include it in the next release.
**
** Wayne Hayes, wayne@csri.utoronto.ca (preffered), or wayne@csri.toronto.edu
*/

#include "misc.h"
#include "sorts.h"
#include "heap.h"	/* for heapsort */
#include "sets.h"
#include <math.h>

#define TEST_SORTS 1
#if TEST_SORTS
int _compareCount;
#endif

int QuickSort(void *a, size_t n, size_t w, pfnCmpFcn compare)
{
    qsort(a, n, w, compare);
    return 0;
}

/* insertion sort. Good for lists that are already "nearly" sorted, ie no
** element is too far from it's sorted position.  Return the number of
** copies that had to be done (which on average would be O(n^2).)
*/
int InsertionSort(void *a, size_t n, size_t w, pfnCmpFcn compare)
{
    int i, numCopies = 0;
    char *A = (char*) a;
    void *x = alloca(w);     /* temporary storage */

    for(i=1; i < n; i++)
    {
	int j = i;  /* j scans left to find i's position */
	memcpy(x, (void*)(A + w*i), w);
	numCopies++;
	while(j >= 1 && compare((void*)(A + w*(j-1)), x) > 0)
	{
	    memcpy((void*)(A + w*j), (void*)(A + w*(j-1)), w);
	    j--;
	    numCopies++;
	}
	memcpy((void*)(A + w*j), x, w);
	numCopies++;
    }
    return numCopies;
}


/* Merge two sorted lists.
*/
static int merge(void *L1, size_t n1,    void *L2, size_t n2,
    size_t width, void *Dest, pfnCmpFcn CmpFcn)
{
    int numCopies = 0;
    int id = 0, i1 = 0, i2 = 0; /* indexes into the lists */
    char *l1 = L1, *l2 = L2, *dest = Dest;
    while(i1 < n1 || i2 < n2)
    {
	char *pd = dest + id * width,   *p1 = l1 + i1 * width,
					*p2 = l2 + i2 * width;
	if(i2 >= n2)    /* list 2 is finished, just add from list 1 */
	{
	    memcpy(pd, p1, width);
	    numCopies++;
	    i1++;
	}
	else if(i1 >= n1)   /* list 1 is finished, just add from list 2 */
	{
	    memcpy(pd, p2, width);
	    numCopies++;
	    i2++;
	}
	else if(CmpFcn(p1, p2) <= 0)
	{
	    memcpy(pd, p1, width);
	    numCopies++;
	    i1++;
	}
	else
	{
	    memcpy(pd, p2, width);
	    numCopies++;
	    i2++;
	}
	id++;
    }
    assert(i1 == n1 && i2 == n2);
    return numCopies;
}

int MergeSort(void *a, size_t n, size_t w, pfnCmpFcn compare)
{
    int numCopies = 0;
    char *A = (char*)a;
    void *tmp = alloca(n*w);
    assert(n >= 0);
    if(n < 2)
	return 0;
    numCopies += MergeSort(a, n/2, w, compare);
    numCopies += MergeSort((void*)(A + (n/2)*w), n-n/2, w, compare);
    numCopies += merge(a, n/2, (void*)(A+(n/2)*w), n-n/2, w, tmp, compare);
    memcpy(a, tmp, n*w);
    ++numCopies;
    return numCopies;
}


/* simple, in-place approximately nlogn sort; see BYTE Apr 1991.
 * I don't think it's been *proven* to be nlogn, although it
 * certainly tends to be, because until gap=1, the outer loop
 * has at most log_{COMB_SHRINK}(n) iterations.  However, often
 * more than one iteration is required at gap=1, and I don't know
 * that anybody has shown that the number of gap=1 iterations is
 * bounded by a constant.
 */
#define COMB_SHRINK 1.3
int CombSort(void *a, size_t n, size_t w, pfnCmpFcn compare)
{
     /* Possible improvement: An initial gap=n/1.3 results in the first gap
      * being about 0.77n.  The for loop, then, only
      * considers elements from 0 through 0.23n for swapping with elements
      * in 0.77n through n.  This means that elements from 0.23n to 0.77n
      * aren't even considered during the first outer loop.  It makes
      * sense (to me) that it would be best if everybody was considered
      * for swapping from the start.  This would mean we want a gap of
      * n/2 in the first iteration, rather than gap/1.3.
      */
    int gap = MAX(n/COMB_SHRINK,1), gap1 = 0;
    char *A = (char*)a;
    void *tmp = alloca(w);
    Boolean done;
    int numCopies = 0;
    do
    {
	int i;
	done = true;
	/* Experiment has shown that the sequence of gaps starting at 11
	 * gives better results than the sequence starting at 9 or 10.
	 */
	if(gap == 9 || gap == 10)
	    gap = 11;
	if(gap == 1)
	    ++gap1;
	for(i=0; i < n-gap; i++)
	    if(compare((void*)(A + w*i), (void*)(A + w*(i+gap))) > 0)
	    {
		memcpy(tmp, (void*)(A + w*i), w);
		memcpy((void*)(A + w*i), (void*)(A + w*(i+gap)), w);
		memcpy((void*)(A + w*(i+gap)), tmp, w);
		numCopies += 3;
		done = false;
	    }
	gap = MAX(gap/COMB_SHRINK,1);
    } until(gap1 > 0 && done);
    /*printf("[gap1=%d]",gap1);*/
    return numCopies;
}

/*
** To make heads or tails of this function, take a look at heap.c.
*/
int HeapSort(void *a, size_t n, size_t w, pfnCmpFcn compare)
{
    int i, numCopies = 0, heapSize;
    char *A = a, *tmp = alloca(w), *top = alloca(w), *bubble = alloca(w);
    A -= w;		/* make A[1] the first element */

    /* The first element is a heap by itself, so start at 2. */
    heapSize = 1;
    for(i=2; i <= n; i++)
    {
	/* HeapInsert(A[i]); */
	int place = ++heapSize;	/* the new element starts at the bottom */
	/*
	** bubble it up.  Note the comparison is backwards; this is
	** necessary since we place things backwards at the next stage,
	** so we want to secretly create a *largest* at top heap, then
	** when we're done, the top of the heap goes to the back, etc.
	*/
	while(place > 1 && compare(A + (place/2)*w, A + place*w) <= 0)
	{
	    memcpy(tmp, A + place*w, w);
	    memcpy(A + place*w, A + (place/2)*w, w);
	    place /= 2;
	    memcpy(A + place*w, tmp, w);
	    numCopies += 3;
	}
    }

    /* Now we have a heap.  Build the sorted array. */
    for(i=0; i < n; i++)
    {
	int place, child1, child2;
	memcpy(top, A+w, w);	/* store the top element for later */
	memcpy(bubble, A + heapSize*w, w);
	memcpy(A+w, bubble, w);
	numCopies += 3;
	--heapSize;
	place = 1;
	/* bubble it down */
	while(child1 = 2*place, child2 = child1 + 1, child1 <= heapSize)
	{
	    int smallestChild = child1;
	    if(child2 <= heapSize && compare(A + child2*w, A + child1*w) >= 0)
		smallestChild = child2;
	    if(compare(bubble, A + smallestChild*w) > 0)	/* done */
		break;

	    /* exchange with smallest child */
	    memcpy(A + place*w, A + smallestChild*w, w);
	    place = smallestChild;
	    memcpy(A + place*w, bubble, w);
	    numCopies += 3;
	}
	/* Finally, put the top in the proper place */
	memcpy(A + (heapSize+1)*w, top, w);
	++numCopies;
    }

    return numCopies;
}


static pfnCmpFcn _PointerSortCompare;
static int _PointerSortCmp(const void *x, const void *y)
{
    const char *const*X = x, *const*Y = y;
    return _PointerSortCompare(*X, *Y);
}

#ifndef PointerSort
int PointerSort(SortFcn Sort, void *a, size_t n, size_t w, pfnCmpFcn compare)
{
    void *tmp = alloca(w);
    char *A = a, *p[n];
    int i, k, start, numCopies = 0;

    if(n < 2)
	return 0;

    for(i=0; i<n; i++)
	p[i] = A+i*w;
    _PointerSortCompare = compare;
    numCopies += Sort(p, n, sizeof(p[0]), _PointerSortCmp);

    /*
    ** This is the permutation array.  I[k] contains the j for which a[j]
    ** should be moved to a[k].  So, for example, if I[0] = 5, then the
    ** smallest element is currently in a[5], and should be moved to a[0],
    ** and of course *(p[0]) == a[5].
    */
#define I(k) ((p[k]-A)/w)
    /*printf("I = "); for(i=0; i<n; i++) printf("%d ", I(i)); printf("\n");*/

    /*
    ** Now do the actual permutation.
    ** start = the first element in a cycle.
    ** p[k] is set to NULL once it's a[k] is properly assigned.
    ** This is O(n).  Each element is touched no more than twice; once
    ** when start touches it, and once when it is touched by a cycle.
    ** We only need to loop to n-1 because the last one is either in
    ** it's proper place, or it was part of a previous cycle.
    */
    for(start = 0; start < n-1; start++)
    {
	if(p[start] == NULL || I(start) == start)
	    continue;

	/* temporarily store the beginning of the cycle, move it to the
	** end later. */
	memcpy(tmp, A + start*w, w); ++numCopies;
	k = start; /* now do the cycle */
	do
	{
	    int oldK = k;
	    memcpy(A + k*w, A + I(k)*w, w); ++numCopies;
	    k = I(k);
	    p[oldK] = NULL;
	} until(I(k) == start);
	memcpy(A + k*w, tmp, w); ++numCopies;
	p[k] = NULL;
    }
#undef I
    return numCopies;
}
#endif

#define FREDERICKS 1
#if FREDERICKS
/*
** This algorithm is called Frederickson's algorithm.  It was outlined
** in the following talk.  I didn't understand the author's extension
** very well, but Frederickson's algorithm is quite trivial.  In essence,
** it is this:
** 
**     Split the input into n/logn blocks of size logn.  Keep a priority
**     queue P of the smallest elements in each block (P has size n/logn).
**     After deleting an element in P, delete the same element from its
**     block, and use a brute-force search of the block for the next
**     smallest element.
** 
** Here is the abstract of the talk.
** 
**            Department of Computer Science, University of Toronto
**          (SF = Sandford Fleming Building, 10 King's College Road)
** 
**        -------------------------------------------------------------
** 
**                                   THEORY
**                 SF1101, at 11:00 a.m., Friday 12 June 1998
** 
**                                 Theis Rauhe
**                         BRICS, University of Aarhus
** 
**                 "Optimal Time-Space Trade-offs for Sorting"
** 
** This talk addresses the fundamental problem of sorting n elements from an
** ordered universe in a sequential model of computation and in particular
** consider the time-space trade-off (product of time and space) for this
** problem.  Beame has shown a lower bound of $Omega(n^2)$ for this product
** leaving a gap of a logarithmic factor up to the previously best known upper
** bound of $O(n^2log n)$ due to Frederickson.
** 
** In this talk a comparison based sorting algorithm which closes this gap is
** presented.  This algorithm obtain the optimal time-space product $O(n^2)$
** for the full range of space bounds between $log n$ and $n/log n$.
*/

static SET *fredDone;
static void *fredArray;
static int fredSize, fredWidth;
static pfnCmpFcn fredCmp;

/*
** Compare Fred things using the array index rather than the objects
** themselves.  Use fredArray and fredCmp to call the "real" comparison
** function.
*/
static int FredCmpIndex(const foint fi, const foint fj)
{
    int i = fi.i, j = fj.i;
    int cmp = fredCmp((char*)fredArray + i*fredWidth,
	(char*)fredArray + j*fredWidth);
    assert(0 <= i && i < fredSize);
    assert(0 <= j && j < fredSize);
    if(cmp == 0)
	return i - j;
    else
	return cmp;
}

/*
** Use brute force to find the found element in a block, starting at
** "start" and going for n elements.  Return index of biggest, relative
** to the entire array, ie, start <= biggest < start + n.
*/
static int BruteSearch(int start, int end)
{
    int i, found = -1;

    for(i=start; i<end; i++)
    {
	if(SetIn(fredDone, i))
	    continue;

	if(found == -1)
	{
	    found = i;
	    continue;
	}

	if(FredCmpIndex((foint)i, (foint)found) < 0)
	    found = i;
    }
    return found;
}

int FredSort(void *a, size_t n, size_t w, pfnCmpFcn compare)
{
    int i, block, numCopies=0;
    int blockSize = Log2(n);	/* floor(log_2(n)) */
    int numBlocks = (n+blockSize-1)/blockSize;
    HEAP *P = HeapAlloc(numBlocks, FredCmpIndex);
    int outputNext = 0;
    char output[n*w];

	int next, nextsBlock, found;
    
    fredArray = a;
    fredSize = n;
    fredWidth = w;
    fredCmp = compare;
    fredDone = SetAlloc(n);

    /* Initialize the Heap */
    for(block=0; block < numBlocks; block++)
    {
	found = BruteSearch(block*blockSize, MIN((block+1)*blockSize, n));
	assert(block*blockSize <= found && found < MIN((block+1)*blockSize, n));
	HeapInsert(P, (foint)found);
    }
    assert(HeapSize(P) == numBlocks);
    
    for(i=0; i<n; i++)
    {
	next = HeapNext(P).i;
	assert(0 <= next && next < n);
	nextsBlock = next / blockSize;
	assert(0 <= nextsBlock && nextsBlock < numBlocks);
	memcpy(output + w*outputNext++, (char*)a + w*next, w);
	++numCopies;

	assert(!SetIn(fredDone, next));
	SetAdd(fredDone, next);
	found=BruteSearch(nextsBlock*blockSize,
	    MIN((nextsBlock+1)*blockSize, n));
	if(found != -1)
	    HeapInsert(P, (foint)found);
    }

    /* sanity checks */
    for(i=0; i<n; i++)
	assert(SetIn(fredDone,i));
    assert(HeapSize(P) == 0);

    /*
    ** Copy the output back into the correct places.
    */
    memcpy(a, output, n*w);
    numCopies += n;
    return numCopies;
}
#endif	/* FREDERICKS */

/***************************   END OF LIBRARY SOURCE *******************/

#if TEST_SORTS
#include <stdio.h>

#define TOTBYTES (1*1024*1024)
#define ElementSize (1)
#define NUM (TOTBYTES/ElementSize)

#define PreSortCmp ElementCmpj
#define SortCmp ElementCmpi

typedef struct _foo {unsigned char i,j; unsigned char f[ElementSize]; } ELEMENT;

int ElementCmp(const void *X, const void *Y)
{
    const ELEMENT *x = X, *y = Y;
    int i;
    _compareCount++;
    if(x->i != y->i) return x->i - y->i;
    if(x->j != y->j) return x->j - y->j;

    for(i=0; i<ElementSize; i++)
	if(x->f[i] != y->f[i]) return x->f[i] - y->f[i];

    return 0;
}

int ElementCmpi(const void *X, const void *Y)
{
    const ELEMENT *x = X, *y = Y;
    _compareCount++;
    return x->i - y->i;
}

int ElementCmpj(const void *X, const void *Y)
{
    const ELEMENT *x = X, *y = Y;
    _compareCount++;
    return x->j - y->j;
}

void ElementDmp(const ELEMENT *x)
{
    int i;
    printf("%d.%d[", x->i, x->j);
    for(i=0; i<ElementSize; i++)
	printf("%d,", x->f[i]);
    printf("]");
}

int main(void)
{
    /*assert(NUM*ElementSize <= 16384*1024);*/
    srand48(time(0) + getpid());
    printf("NUM %d, ElementSize %d\n", NUM, sizeof(ELEMENT));
do
{
    ELEMENT a[NUM], b[NUM], c[NUM], d[NUM], e[NUM], f[NUM], g[NUM];
    int i, cpCnt;
    double t;
    Boolean stable;

    printf("generating random list..."); fflush(stdout);

    for(i=0; i<NUM; i++)
    {
	int j;
	a[i].i = lrand48();
	a[i].j = lrand48();
	for(j=0; j<ElementSize; j += 4)
	{
	    unsigned l = lrand48();
	    a[i].f[j] = l;
	    a[i].f[j+1] = (l <<= 8);
	    a[i].f[j+2] = (l <<= 8);
	    a[i].f[j+3] = (l <<= 8);
	}
	/*ElementDmp(a+i); printf("\n");*/
    }

    memcpy(b, a, sizeof(ELEMENT)*NUM);
    memcpy(c, a, sizeof(ELEMENT)*NUM);
    memcpy(d, a, sizeof(ELEMENT)*NUM);
    memcpy(e, a, sizeof(ELEMENT)*NUM);
    memcpy(f, a, sizeof(ELEMENT)*NUM);
    memcpy(g, a, sizeof(ELEMENT)*NUM);

    printf("Starting sorts\n");

    printf("qsort: "); fflush(stdout);
    _compareCount = 0;
    t = uTime();
    cpCnt = PointerSort(QuickSort,a, NUM, sizeof(ELEMENT), PreSortCmp);
    cpCnt += PointerSort(QuickSort,a, NUM, sizeof(ELEMENT), SortCmp);
    stable = true;
    for(i=0; i<NUM-1; i++)
    {
	assert(SortCmp(a+i, a+i+1) <= 0);
	if(SortCmp(a+i, a+i+1) == 0 && PreSortCmp(a+i, a+i+1) > 0)
	{
	    stable = false;
	    break;
	}
    }
    printf("cm: %9d cp: %9d (%gs) %sstable\n", _compareCount, cpCnt, uTime()-t,
	stable?"":"un");

    printf("csort: "); fflush(stdout);
    _compareCount = 0;
    t = uTime();
    cpCnt = PointerSort(CombSort,b, NUM, sizeof(ELEMENT), PreSortCmp);
    cpCnt += PointerSort(CombSort,b, NUM, sizeof(ELEMENT), SortCmp);
    stable = true;
    for(i=0; i<NUM-1; i++)
    {
	assert(SortCmp(b+i, b+i+1) <= 0);
        /*ElementDmp(b+i); printf("\n");*/
	assert(SortCmp(a+i,b+i) == 0);
	if(SortCmp(b+i, b+i+1) == 0 && PreSortCmp(b+i, b+i+1) > 0)
	{
	    stable = false;
	    break;
	}
    }
    printf("cm: %9d cp: %9d (%gs) %sstable\n", _compareCount, cpCnt, uTime()-t,
	stable?"":"un");

    printf("msort: "); fflush(stdout);
    _compareCount = 0;
    t = uTime();
    cpCnt = PointerSort(MergeSort,c, NUM, sizeof(ELEMENT), PreSortCmp);
    cpCnt += PointerSort(MergeSort,c, NUM, sizeof(ELEMENT), SortCmp);
    stable = true;
    for(i=0; i<NUM-1; i++)
    {
	assert(SortCmp(c+i, c+i+1) <= 0);
	/*ElementDmp(c+i); printf("\n");*/
	assert(SortCmp(a+i,c+i) == 0);
	if(SortCmp(c+i, c+i+1) == 0 && PreSortCmp(c+i, c+i+1) > 0)
	{
	    stable = false;
	    break;
	}
    }
    printf("cm: %9d cp: %9d (%gs) %sstable\n", _compareCount, cpCnt, uTime()-t,
	stable?"":"un");

    printf("hsort: "); fflush(stdout);
    _compareCount = 0;
    t = uTime();
    cpCnt = PointerSort(HeapSort,d, NUM, sizeof(ELEMENT), PreSortCmp);
    cpCnt += PointerSort(HeapSort,d, NUM, sizeof(ELEMENT), SortCmp);
    stable = true;
    for(i=0; i<NUM-1; i++)
    {
	assert(SortCmp(d+i, d+i+1) <= 0);
	/*ElementDmp(d+i); printf("\n");*/
	assert(SortCmp(a+i,d+i) == 0);
	if(SortCmp(d+i, d+i+1) == 0 && PreSortCmp(d+i, d+i+1) > 0)
	{
	    stable = false;
	    break;
	}
    }
    printf("cm: %9d cp: %9d (%gs) %sstable\n", _compareCount, cpCnt, uTime()-t,
	stable?"":"un");

#if 0
    printf("isort: "); fflush(stdout);
    _compareCount = 0;
    t = uTime();
    cpCnt = PointerSort(InsertionSort, e, NUM, sizeof(ELEMENT), PreSortCmp);
    cpCnt += PointerSort(InsertionSort, e, NUM, sizeof(ELEMENT), SortCmp);
    stable = true;
    for(i=0; i<NUM-1; i++)
    {
	assert(SortCmp(e+i, e+i+1) <= 0);
	/*ElementDmp(e+i); printf("\n");*/
	assert(SortCmp(a+i,e+i) == 0);
	if(SortCmp(e+i, e+i+1) == 0 && PreSortCmp(e+i, e+i+1) > 0)
	{
	    stable = false;
	    break;
	}
    }
    printf("cm: %9d cp: %9d (%gs) %sstable\n", _compareCount, cpCnt, uTime()-t,
	stable?"":"un");
#endif

    printf("fsort: "); fflush(stdout);
    _compareCount = 0;
    t = uTime();
    cpCnt = FredSort(f, NUM, sizeof(ELEMENT), PreSortCmp);
    cpCnt += FredSort(f, NUM, sizeof(ELEMENT), SortCmp);
    stable = true;
    for(i=0; i<NUM-1; i++)
    {
	assert(SortCmp(f+i, f+i+1) <= 0);
        /*ElementDmp(b+i); printf("\n");*/
	assert(SortCmp(a+i,f+i) == 0);
	if(SortCmp(f+i, f+i+1) == 0 && PreSortCmp(f+i, f+i+1) > 0)
	{
	    stable = false;
	    break;
	}
    }
    printf("cm: %9d cp: %9d (%gs) %sstable\n", _compareCount, cpCnt, uTime()-t,
	stable?"":"un");

#if 0	/* David Neto's FredSort */
    printf("dfsort: "); fflush(stdout);
    _compareCount = 0;
    t = uTime();
    cpCnt = FredSort(g, NUM, sizeof(ELEMENT), PreSortCmp);
    cpCnt += FredSort(g, NUM, sizeof(ELEMENT), SortCmp);
    stable = true;
    for(i=0; i<NUM-1; i++)
    {
	assert(SortCmp(g+i, g+i+1) <= 0);
        /*ElementDmp(b+i); printf("\n");*/
	assert(SortCmp(a+i,g+i) == 0);
	if(SortCmp(g+i, g+i+1) == 0 && PreSortCmp(g+i, g+i+1) > 0)
	{
	    stable = false;
	    break;
	}
    }
    printf("cm: %9d cp: %9d (%gs) %sstable\n", _compareCount, cpCnt, uTime()-t,
	stable?"":"un");
#endif
} while(1);
    return 0;
}
#endif	/* TEST_SORTS */
