/*
    Theseus - maximum likelihood superpositioning of macromolecular structures

    Copyright (C) 2004-2015 Douglas L. Theobald

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the:

    Free Software Foundation, Inc.,
    59 Temple Place, Suite 330,
    Boston, MA  02111-1307  USA

    -/_|:|_|_\-
*/

/* -/_|:|_|_\- */
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "quicksort.h"
/* 15 is pretty good -- any array segment below this gets insorted */
#ifndef CUTOFF
#  define CUTOFF 15
#endif

extern void
swapd(double *x, double *y)
{
    double          tempd;

    tempd = *x;
    *x = *y;
    *y = tempd;
}

extern void
swapc(char **x, char **y)
{
    char           *tempc;

    tempc = *x;
    *x = *y;
    *y = tempc;
}

/*
 |  void
 |  partial_quicksort (KEY_T  *array1, COARRAY_T  *array2, int lower, int upper)
 |
 |  Abstract:
 |  Sort array1[lower..upper] into a partial order
 |      leaving segments which are CUTOFF elements long
 |      unsorted internally.
 |
 |  Efficiency:
 |  I use a randomly generated pivot, which ensures maximum
 |  average efficiency. Also, you could use median of three
 |  method to choose a pivot, or just use the first item 
 |  in the array (very bad if array is almost presorted).
 |
 |  Method:
 |  Partial Quicksort with a sentinel (Robert Sedgewick)
 |
 |  !! NOTE:
 |  array1[upper+1] must hold the maximum possible key.
 */

void
partial_quicksort2 (KEY_T  *array1, COARRAY_T  *array2, int lower, int upper)
{
    int             i, j, random_index;
    KEY_T           pivot;

    if (upper - lower > CUTOFF)
    {
        swapd(&array1[lower], &array1[(upper+lower)/2]);
        swapc(&array2[lower], &array2[(upper+lower)/2]);
        i = lower;
        j = upper + 1;

        /* pivot = array1[lower]; */
        srand(time(NULL));
        random_index = (int) ( (double)(upper - lower) * (double) rand() / (double) RAND_MAX ) + lower;
        pivot = array1[random_index];

        while (1)
        {
            /*
             * ---------------------- !! NOTE ----------------------
             * ignoring NOTE above can lead to an infinite loop here
             * -----------------------------------------------------
             */
            do
                i++;
                while (LT(array1[i], pivot));
            do
                j--;
                while (GT(array1[j], pivot));

            if (j > i)
            {
                swapd(&array1[i], &array1[j]);
                swapc(&array2[i], &array2[j]);
            }
            else
                break;
        }

        swapd(&array1[lower], &array1[j]);
        swapc(&array2[lower], &array2[j]);

        partial_quicksort2 (array1, array2, lower, j - 1);
        partial_quicksort2 (array1, array2, i,     upper);
    }
}


void
partial_quicksort2d (KEY_T  *array1, KEY_T  *array2, int lower, int upper)
{
    int             i, j, random_index;
    KEY_T           pivot;

    if (upper - lower > CUTOFF)
    {
        swapd(&array1[lower], &array1[(upper+lower)/2]);
        swapd(&array2[lower], &array2[(upper+lower)/2]);
        i = lower;
        j = upper + 1;

        /* pivot = array1[lower]; */
        srand(time(NULL));
        random_index = (int) ( (double)(upper - lower) * (double) rand() / (double) RAND_MAX ) + lower;
        pivot = array1[random_index];

        while (1)
        {
            /*
             * ---------------------- !! NOTE ----------------------
             * ignoring NOTE above can lead to an infinite loop here
             * -----------------------------------------------------
             */
            do
                i++;
                while (LT(array1[i], pivot));
            do
                j--;
                while (GT(array1[j], pivot));

            if (j > i)
            {
                swapd(&array1[i], &array1[j]);
                swapd(&array2[i], &array2[j]);
            }
            else
                break;
        }

        swapd(&array1[lower], &array1[j]);
        swapd(&array2[lower], &array2[j]);

        partial_quicksort2d (array1, array2, lower, j - 1);
        partial_quicksort2d (array1, array2, i,     upper);
    }
}


void
partial_quicksort (KEY_T  *array, int lower, int upper)
{
    int             i, j, random_index;
    KEY_T           pivot;

    if (upper - lower > CUTOFF)
    {
        swapd(&array[lower], &array[(upper+lower)/2]);
        i = lower;
        j = upper + 1;

        /* pivot = array1[lower]; */
        srand(time(NULL));
        random_index = (int) ( (double)(upper - lower) * (double) rand() / (double) RAND_MAX ) + lower;
        pivot = array[random_index];

        while (1)
        {
            /*
             * ---------------------- !! NOTE ----------------------
             * ignoring NOTE above can lead to an infinite loop here
             * -----------------------------------------------------
             */
            do
                i++;
                while (LT(array[i], pivot));
            do
                j--;
                while (GT(array[j], pivot));

            if (j > i)
            {
                swapd(&array[i], &array[j]);
            }
            else
                break;
        }

        swapd(&array[lower], &array[j]);

        partial_quicksort (array, lower, j - 1);
        partial_quicksort (array, i,     upper);
    }
}

/*
 |  void  insort (KEY_T array1[], int len)
 |
 |  Abstract:   Sort array1[0..len-1] into increasing order.
 |
 |  Method: Optimized insertion-sort (ala Jon Bentley)
 */

void
insort2 (KEY_T *array1, COARRAY_T  *array2, int len)
{
    int             i, j;
    KEY_T           temp1;
    COARRAY_T       temp2;

    for (i = 1; i < len; ++i)
    {
        j = i;
        temp1 = array1[j];
        temp2 = array2[j];

        while ((j > 0) && GT(array1[j-1], temp1))
        {
            array1[j] = array1[j-1];
            array2[j] = array2[j-1];
            j--;
        }

        array1[j] = temp1;
        array2[j] = temp2;
    }
}

void
insort2d (KEY_T *array1, KEY_T  *array2, int len)
{
    int             i, j;
    KEY_T           temp1;
    KEY_T           temp2;

    for (i = 1; i < len; ++i)
    {
        j = i;
        temp1 = array1[j];
        temp2 = array2[j];

        while ((j > 0) && GT(array1[j-1], temp1))
        {
            array1[j] = array1[j-1];
            array2[j] = array2[j-1];
            j--;
        }

        array1[j] = temp1;
        array2[j] = temp2;
    }
}

void
insort (KEY_T *array, int len)
{
    int             i, j;
    KEY_T           temp;

    for (i = 1; i < len; ++i)
    {
        j = i;
        temp = array[j];

        while ((j > 0) && GT(array[j-1], temp))
        {
            array[j] = array[j-1];
            j--;
        }

        array[j] = temp;
    }
}

/*
 |  void  quicksort (KEY_T  array1[], int len)
 |
 |  Abstract:
 |  Sort array1[0..len-1] into increasing order.
 |
 |  Method:
 |  Use partial_quicksort() with a sentinel (ala Sedgewick)
 |  to reach a partial order, leave the unsorted segments of
 |  length <= CUTOFF to low-overhead straight insertion sort.
 |
 |  !! NOTE:
 |  array1[len] must be a "sentinel" -- the largest
 |  possible value (e.g, array1[len] = HUGE_VAL;)
 */
void
quicksort2 (KEY_T  *array1, COARRAY_T  *array2, int len)
{
    partial_quicksort2 (array1, array2, 0, len - 1);
    insort2 (array1, array2, len);
}

void
quicksort2d (KEY_T  *array1, KEY_T  *array2, int len)
{
    partial_quicksort2d (array1, array2, 0, len - 1);
    insort2d (array1, array2, len);
}

void
quicksort (KEY_T  *array, int len)
{
    partial_quicksort (array, 0, len - 1);
    insort (array, len);
}