/*
 * sort.c
 *      Comparing internal sort methods
 */
#ifdef WIN32
#include <windows.h>
typedef unsigned int uint;
#define srandom(x) srand((x))
#else
#include <unistd.h>
#include <sys/time.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define lengthof(array) \
	(sizeof (array) / sizeof ((array)[0]))

void            qsortPDQ(void *a, size_t n, size_t es, int (*cmp) (const void *, const void *));


void            qsortB(void *base, size_t nel, size_t width,
                                int (*compar) (const void *, const void *));

void qsb(void *base, size_t n, size_t size, int (*compar) (const void *, const void *));

typedef struct sort_t {
    char            name[32];
    void            (*fp_qsort) (void *base, size_t nel, size_t width,
                                int (*compar) (const void *, const void *));
}
                sort_t;

/* testing matrix */
static const sort_t d_method[] =
{
    {"glibc qsort", qsort},
    {"BSD qsort", qsortB},
    {"qsb", qsb},
    {"qsortPDQ", qsortPDQ}
};
static const size_t d_nelem[] =
{1000, 10000, 100000, 1000000, 5000000};
static const size_t d_range[] =
{2, 32, 1024, 0xFFFF0000L};
static const char *d_distr[] =
{"uniform", "gaussian", "95sorted", "95reversed"};
static const size_t d_ccost[] =
{2};

/* factor index */
static int      n,
                r,
                d,
                c;

static volatile int dummy;

static          uint
                uniform(uint min, uint max)
{
    return (min + (uint) (max * 1.0 * rand() / (RAND_MAX + 1.0)));
}

/* copied from C-FAQ 13.20 */
static double   gaussrand()
{
    static double   V1,
                    V2,
                    S;
    static int      phase = 0;
    double          X;

    if (phase == 0) {
        do {
            double          U1 = (double) rand() / RAND_MAX;
            double          U2 = (double) rand() / RAND_MAX;

            V1 = 2 * U1 - 1;
            V2 = 2 * U2 - 1;
            S = V1 * V1 + V2 * V2;
        }
        while (S >= 1 || S == 0);

        X = V1 * sqrt(-2 * log(S) / S);
    } else
        X = V2 * sqrt(-2 * log(S) / S);

    phase = 1 - phase;

    return X;
}

/* example of ruinning Gaussian ... */
static          uint
                gaussian(uint min, uint max)
{
    double          X = gaussrand();

    /* 99% X is in [-2.57, 2.57] */
    X = ((double) (min + max)) / 2 + X / 2.57 * (((double) (max - min)) / 2);
    if (X < min)
        return min;
    else if (X > max)
        return max;
    else
        return floor(X + .5);
}

static void     cpu_cost(size_t ccost)
{
    while (ccost--)
        dummy *= 37;
}

static int      key_comp(const void *i, const void *j)
{
    cpu_cost(d_ccost[c]);
    return (*(uint *) (*(uint **) i) - *(uint *) (*(uint **) j));
}

static int      _key_comp(const void *i, const void *j)
{
    return -key_comp(i, j);
}

#define swap(a, b) \
do{ \
        uint d = (a);        \
        (a) = (b);           \
        (b) = d;             \
}while(0)

static char    *
                generate_data()
{
    char           *array;
    uint          **parray;
    uint           *karray;

    size_t          nelem,
                    range;
    const char     *distr;
    int             i;

    nelem = d_nelem[n];
    range = d_range[r];
    distr = d_distr[d];

    /* int is used to store key, int * points to the key */
    array = (char *) malloc((sizeof(uint *) + sizeof(uint)) * nelem);
    if (array == NULL) {
        fprintf(stderr, "insufficient memory: nelem = %d", nelem);
        exit(-1);
    }
    parray = (uint **) array;
    karray = (uint *) (parray + nelem);

    for (i = 0; i < nelem; i++) {
        if (!strcmp(distr, "gaussian"))
            karray[i] = gaussian(1, range);
        else
            karray[i] = uniform(1, range);

        parray[i] = karray + i;
    }

    /* check if should do the 95% thing */
    if (!strcmp(distr, "95sorted")) {
        qsort(parray, nelem, sizeof(uint *), key_comp);

        /* 95% gaussian is in [-1.96, 1.96] */
        for (i = 0; i < nelem; i++) {
            double          r = gaussrand();
            if (r < -1.96 || r > 1.96)
                swap(karray[i], karray[nelem - i - 1]);
        }

    } else if (!strcmp(distr, "95reversed")) {
        qsort(parray, nelem, sizeof(uint *), _key_comp);
        for (i = 0; i < nelem; i++) {
            double          r = gaussrand();
            if (r < -1.96 || r > 1.96)
                swap(karray[i], karray[nelem - i - 1]);
        }
    }
    return array;
}

static void     check_sort(void *parray)
{
    int             i;
    for (i = 0; i < d_nelem[n] - 1; i++)
        if (key_comp((uint **) parray + i, (uint **) parray + i + 1) > 0) {
            fprintf(stderr, "Oooops ... \n");
            exit(-1);
        }
}

static void     do_sort(void *parray, int m)
{
    d_method[m].fp_qsort(parray, d_nelem[n], sizeof(uint *), key_comp);
    check_sort(parray);
}

int             main(int argc, char *argv[])
{
    int             rounds;
    struct timeval  start_t,
                    stop_t;

    rounds = 1;
    for (n = 0; n < lengthof(d_nelem); n++)
        for (r = 0; r < lengthof(d_range); r++)
            for (d = 0; d < lengthof(d_distr); d++)
                for (c = 0; c < lengthof(d_ccost); c++) {
                    int             m;
                    char           *array;

                    for (m = 0; m < lengthof(d_method); m++) {
                        /* use the same random seed */
                        srandom(getpid() + rounds);
                        array = generate_data();

                        gettimeofday(&start_t, NULL);
                        do_sort(array, m);
                        gettimeofday(&stop_t, NULL);

                        /* measure the time */
                        if (stop_t.tv_usec < start_t.tv_usec) {
                            stop_t.tv_sec--;
                            stop_t.tv_usec += 1000000;
                        }
                        fprintf(stdout, "[%d] [%s]: nelem(%d), range(%u) distr(%s) ccost(%d) : %ld.%03ld ms \n",
                                rounds,
                                d_method[m].name,
                             d_nelem[n], d_range[r], d_distr[d], d_ccost[c],
                           (long) ((stop_t.tv_sec - start_t.tv_sec) * 1000 +
                                 (stop_t.tv_usec - start_t.tv_usec) / 1000),
                          (long) (stop_t.tv_usec - start_t.tv_usec) % 1000);

                        free(array);
                    }

                    fprintf(stdout, "------\n");
                    rounds++;
                }

    fprintf(stdout, "done\n");
    exit(0);
}