/*** HELP START ***//*
 
## >>> `qsortInCbyProcProto()` proto function: <<< <a name="qsortincbyprocproto-proto-function"></a> #######################  

The **qsortInCbyProcProto()** is external *C* function, 
this is the implementation of the *Quick Sort* algorithm. 

The function is used **internally** by 
functions in the *BasePlus* package.

Asumptions:
- smaller subarray is sorted first, 
- subarrays of *size < 11* are sorted by *insertion sort*, 
- pivot is selected as median of low index value, 
  high index value, and (low+high)/2 index value.

`!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!`<br>
`!CAUTION! Sorted array CANNOT contains SAS missing values !`<br>
`!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!`<br>

### SYNTAX: ###################################################################

The basic syntax is the following:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~sas
qsortInCbyProcProto(arr, low, high)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

**Arguments description**:

1. `arr`  - An array of double type to be sorted.

2. `low`  - An integer low index of starting position (from which the sorting is done).

3. `high` - An integer high index of ending position (up to which the sorting is done).

*//*** HELP END ***/

void qsortInCbyProcProto(double *arr/iotype=u, int low, int high)
  label = "This is the implementation of the Quick sort algorithm";


externc qsortInCbyProcProto;

/* A utility function to swap two elements */
void swapX(double *x,double *y)
{
    double temp;
    temp = *x;
    *x = *y;
    *y = temp;
}

/* A utility function to find middle of three elements */
int median3(double *arr, int l, int h)
{
  int k;
  k = ((l + h) / 2);

  if (&arr[k] > &arr[h]) {
      if (&arr[h] > &arr[l]) return h;
      else { if (&arr[k] > &arr[l]) return l; } 
    } 
  else {
      if (&arr[l] > &arr[h]) return h;  
      else { if (&arr[l] > &arr[k]) return l; }
    }
  return k;
}

/* A utility function to execute partition step */
int partition(double *arr, int l, int h) 
{ 
    int i,j,k;
    double x;

    /* find the pivot */ 
    k = median3(arr, l, h);
    if (h != k) { swapX(&arr[h],&arr[k]); }
    
    x = arr[h]; 
    i = (l - 1); 
  
    for (j = l; j <= h - 1; j++) { 
        if (arr[j] <= x) { 
            i++; 
            swapX(&arr[i], &arr[j]); 
        } 
    } 
    swapX(&arr[i + 1], &arr[h]); 
    return (i + 1); 
} 

/*** HELP START ***//*
 
### REFERENCES: ####################################################

*Reference 1.*

Insertion sort for arrays smaller then 11 elements: 

Based on the code from the following WikiBooks page [2020.08.14]:

[https://pl.wikibooks.org/wiki/Kody_%C5%BAr%C3%B3d%C5%82owe/Sortowanie_przez_wstawianie](https://pl.wikibooks.org/wiki/Kody_%C5%BAr%C3%B3d%C5%82owe/Sortowanie_przez_wstawianie)
  
*//*** HELP END ***/

void insertSort(double *arr, int l, int h)
{
    int i, j;
    double value;

    for (i = l + 1; i <= h; ++i) 
    {
        value = arr[i];
        for (j = i - 1; j >= 0 && arr[j] > value; --j) 
        {
            arr[j + 1] = arr[j];
        }
        arr[j + 1] = value;
    }
}

 
/* Arr[] --> Array to be sorted,  
      l  --> Starting index,  
      h  --> Ending index 
*/
/*** HELP START ***//*
 
*Reference 2.*

Iterative Quick Sort:  

Based on the code from the following pages [2020.08.14]:

[https://www.geeksforgeeks.org/iterative-quick-sort/](https://www.geeksforgeeks.org/iterative-quick-sort/)

[https://www.geeksforgeeks.org/c-program-for-iterative-quick-sort/](https://www.geeksforgeeks.org/c-program-for-iterative-quick-sort/)
  
---

*//*** HELP END ***/

void qsortInCbyProcProto(double *arr, int l, int h) 
{ 
    int p, ls, rs; 
      
    /* Create an auxiliary stack */
    /* since smaller subarray is sorted first 
       the stack should not exceed log2(#Arr) 
       for 1e9 element array it should be up 
       to 2*30 = 60. With 1000 it should cover 
       array of size up to 2*500 (i.e. ~3.273391e+150).
    */
    int stack[10000]; /*increased from 1000 to 10000 */
  
    /* initialize top of stack */
    int top = -1; 
  
    /* push initial values of l and h to stack */
    stack[++top] = l; 
    stack[++top] = h; 
  
    /* Keep popping from stack while is not empty */
    while (top >= 0) { 
        /* Pop h and l */
        h = stack[top--]; 
        l = stack[top--]; 
  
        /* Set pivot element at its correct position 
           in sorted array */
        p = partition(arr, l, h); 
        
        ls = ((p - 1) - l); /* calculate the size of the left side */
        rs = (h - (p + 1)); /* calculate the size of the right size  */

        if (ls > rs) /* to sort smaller part first */
          {
            /* If there are elements on left side of pivot, 
               then push left side to stack */
            if (ls > 10) { 
                stack[++top] = l; 
                stack[++top] = p - 1; 
            }
            else {if (ls > 0) { insertSort(arr, l, p - 1); } } /* sort "by hand" */
            /*##*/

            /* If there are elements on right side of pivot, 
               then push right side to stack */
            if (rs > 10) { 
                stack[++top] = p + 1; 
                stack[++top] = h; 
            }
            else {if (rs > 0) { insertSort(arr, p + 1, h); } } /* sort "by hand" */ 
          }
        else 
          {
            /* If there are elements on right side of pivot, 
               then push right side to stack */
            if (rs > 10) { 
                stack[++top] = p + 1; 
                stack[++top] = h; 
            }
            else {if (rs > 0) { insertSort(arr, p + 1, h); } } /* sort "by hand" */ 
            /*##*/
            /* If there are elements on left side of pivot, 
               then push left side to stack */
            if (ls > 10) { 
                stack[++top] = l; 
                stack[++top] = p - 1; 
            }
            else {if (ls > 0) { insertSort(arr, l, p - 1);  } } /* sort "by hand" */
          } 
    }
    
    top = -1;
} 
externcend;
