/* sorttemplates.c version 2.0, Mar 13, 2018.
 * Author: Brendan McKay; brendan.mckay@anu.edu.au
 *
 * This file contains templates for creating in-place sorting procedures
 * for different data types.  It cannot be compiled separately but
 * should be #included after defining a few preprocessor variables.
 *   SORT_OF_SORT, SORT_NAME and SORT_TYPE1 are always required, and
 *   SORT_TYPE2 is needed for SORT_OF_SORT = 3,
 *   SORT_COMPARE is needed for SORT_OF_SORT = 4.
 * 
 *   SORT_OF_SORT = 1: Creates a procedure
 *         static void SORT_NAME(SORT_TYPE1 *x, int n)
 *     which permutes x[0..n-1] so that x[0] <= ... <= x[n-1].
 *   SORT_OF_SORT = 2: Creates a procedure
 *         static void SORT_NAME(SORT_TYPE1 *x, SORT_TYPE2 *y, int n)
 *     which permutes x[0..n-1] so that x[0] <= ... <= x[n-1]
 *     and also permutes y[0..n-1] by the same permutation.
 *   SORT_OF_SORT = 3: Creates a procedure
 *         static void SORT_NAME(SORT_TYPE1 *x, SORT_TYPE2 *y, int n)
 *     which permutes x[0..n-1] so that y[x[0]] <= ... <= y[x[n-1]].
 *   SORT_OF_SORT = 4: Creates a procedure
 *         static void SORT_NAME(SORT_TYPE1 *x, int n)
 *     which permutes x[0..n-1] so that
 *     SORT_COMPARE(x[i],x[j]) <= 0 for i<j.
 *   SORT_COMPARE(x,y) is an expression <0,0,>0 for x<y,x=y,x>y.
 * 
 *   SORT_NAME = the name of the procedure to be created
 *
 *   SORT_TYPE1 = type of the first or only array (no default)
 *   	 This can be any numeric type for SORT_OF_SORT=1,2, but
 *   	 should be an integer type for SORT_OF_SORT=3.
 *       For SORT_OF_SORT=4, SORT_TYPE1 should be assignable;
 *         such as a simple type or a structure.
 *   SORT_TYPE2 = type of the second array if needed (no default)
 *       This can be any assignable type (including a structure) for
 *       SORT_OF_SORT=2, but must be a numeric type for SORT_OF_SORT=3.
 *
 *   Note that C doesn't like SORT_TYPE1 and SORT_TYPE2 to be defined
 *   as explicit pointer types like int*.  The simplest work-around
 *   is to define a type, e.g. typedef intptr int*  then use that.
 *
 *   SORT_MINPARTITION = least number of elements for using quicksort
 *           partitioning, otherwise insertion sort is used (default "11") 
 *   SORT_MINMEDIAN9 = least number of elements for using the median of 3
 *           medians of 3 for partitioning (default "320")
 *   SORT_FUNCTYPE = type of sort function (default "static void")
 *
 *   This file can be included any number of times provided the value
 *   of SORT_NAME is different each time.
 */
 
#define SORT_MEDIAN_OF_3(a,b,c) \
 ((a) <= (b) ? ((b) <= (c) ? (b) : (c) <= (a) ? (a) : (c)) \
             : ((a) <= (c) ? (a) : (c) <= (b) ? (b) : (c)))
#define F_SORT_MEDIAN_OF_3(f,a,b,c) \
 (f(a,b) <= 0 ? (f(b,c) <= 0 ? (b) : f(c,a) <= 0 ? (a) : (c)) \
             : (f(a,c) <= 0 ? (a) : f(c,b) <= 0 ? (b) : (c)))

#if !defined(SORT_OF_SORT) || !defined(SORT_NAME)
 #error Either SORT_OF_SORT or SORT_NAME is undefined
#endif

#if (SORT_OF_SORT < 1) || (SORT_OF_SORT > 4)
 #error Unknown value of SORT_OF_SORT
#endif

#ifndef SORT_TYPE1
 #error "SORT_TYPE1 must be defined before including sorttemplates.c"
#endif

#ifndef SORT_MINPARTITION
#define SORT_MINPARTITION 11
#endif

#ifndef SORT_MINMEDIAN9
#define SORT_MINMEDIAN9 320
#endif

#ifndef SORT_FUNCTYPE
#define SORT_FUNCTYPE static void
#endif

#define SORT_SWAP1(x,y) tmp1 = x; x = y; y = tmp1;
#define SORT_SWAP2(x,y) tmp2 = x; x = y; y = tmp2;

/*******************************************************************/

#if SORT_OF_SORT == 1
SORT_FUNCTYPE
SORT_NAME(SORT_TYPE1 *x, int n)
{
    int i,j;
    int a,d,ba,dc,s,nn;
    SORT_TYPE1 tmp1,v,v1,v2,v3;
    SORT_TYPE1 *x0,*xa,*xb,*xc,*xd,*xh,*xl;
    struct {SORT_TYPE1 *addr; int len;} stack[40];
    int top;

    top = 0;
    if (n > 1)
    {
        stack[top].addr = x;
        stack[top].len = n;
        ++top;
    }

    while (top > 0)
    {
	--top;
	x0 = stack[top].addr;
	nn = stack[top].len;

        if (nn < SORT_MINPARTITION)
        {
            for (i = 1; i < nn; ++i)
            {
                tmp1 = x0[i];
                for (j = i; x0[j-1] > tmp1; )
                {
                    x0[j] = x0[j-1];
                    if (--j == 0) break;
                }
                x0[j] = tmp1;
            }
	    continue;
        }

        if (nn < SORT_MINMEDIAN9)
            v = SORT_MEDIAN_OF_3(x0[0],x0[nn/2],x0[nn-1]);
        else
        {
	    v1 = SORT_MEDIAN_OF_3(x0[0],x0[1],x0[2]);
	    v2 = SORT_MEDIAN_OF_3(x0[nn/2-1],x0[nn/2],x0[nn/2+1]);
	    v3 = SORT_MEDIAN_OF_3(x0[nn-3],x0[nn-2],x0[nn-1]);
	    v = SORT_MEDIAN_OF_3(v1,v2,v3);
        }

        xa = xb = x0;  xc = xd = x0+(nn-1);
        for (;;)
        {
	    while (xb <= xc && *xb <= v)
	    {
	        if (*xb == v)
	        {
		    *xb = *xa; *xa = v; ++xa;
	        }
	        ++xb;
	    }
	    while (xc >= xb && *xc >= v)
	    {
	        if (*xc == v)
	        {
		    *xc = *xd; *xd = v; --xd;
	        }
	        --xc;
	    }
	    if (xb > xc) break;
            SORT_SWAP1(*xb,*xc);
	    ++xb;
	    --xc;
        }
    
        a = xa - x0;
        ba = xb - xa;
        if (ba > a) s = a; else s = ba;
        for (xl = x0, xh = xb-s; s > 0; --s)
        {
	    *xl = *xh; *xh = v; ++xl; ++xh;
        }
        d = xd - x0;
        dc = xd - xc;
        if (dc > nn-1-d) s = nn-1-d; else s = dc;
        for (xl = xb, xh = x0 + (nn-s); s > 0; --s)
        {
	    *xh = *xl; *xl = v; ++xl; ++xh;
        }

	if (ba > dc)
	{
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	}
	else
	{
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	}
    }
}
#endif

#if SORT_OF_SORT == 2
#ifndef SORT_TYPE2
 #error "SORT_TYPE2 must be defined before including sorttemplates.c"
#endif

SORT_FUNCTYPE
SORT_NAME(SORT_TYPE1 *x, SORT_TYPE2 *y, int n)
{
    int i,j;
    int a,d,ba,dc,s,nn;
    SORT_TYPE2 tmp2,*y0,*ya,*yb,*yc,*yd,*yl,*yh;
    SORT_TYPE1 tmp1,v,v1,v2,v3;
    SORT_TYPE1 *x0,*xa,*xb,*xc,*xd,*xh,*xl;
    struct {SORT_TYPE1 *addr; int len;} stack[40];
    int top;

    top = 0;
    if (n > 1)
    {
        stack[top].addr = x;
        stack[top].len = n;
        ++top;
    }

    while (top > 0)
    {
	--top;
	x0 = stack[top].addr;
	y0 = y + (x0-x);
	nn = stack[top].len;

        if (nn < SORT_MINPARTITION)
        {
            for (i = 1; i < nn; ++i)
            {
                tmp1 = x0[i];
                tmp2 = y0[i];
                for (j = i; x0[j-1] > tmp1; )
                {
                    x0[j] = x0[j-1];
                    y0[j] = y0[j-1];
                    if (--j == 0) break;
                }
                x0[j] = tmp1;
                y0[j] = tmp2;
            }
	    continue;
        }

        if (nn < SORT_MINMEDIAN9)
            v = SORT_MEDIAN_OF_3(x0[0],x0[nn/2],x0[nn-1]);
        else
        {
	    v1 = SORT_MEDIAN_OF_3(x0[0],x0[1],x0[2]);
	    v2 = SORT_MEDIAN_OF_3(x0[nn/2-1],x0[nn/2],x0[nn/2+1]);
	    v3 = SORT_MEDIAN_OF_3(x0[nn-3],x0[nn-2],x0[nn-1]);
	    v = SORT_MEDIAN_OF_3(v1,v2,v3);
        }

        xa = xb = x0;  xc = xd = x0+(nn-1);
        ya = yb = y0;  yc = yd = y0+(nn-1);
        for (;;)
        {
	    while (xb <= xc && *xb <= v)
	    {
	        if (*xb == v)
	        {
		    *xb = *xa; *xa = v; ++xa;
		    SORT_SWAP2(*ya,*yb); ++ya;
	        }
	        ++xb; ++yb;
	    }
	    while (xc >= xb && *xc >= v)
	    {
	        if (*xc == v)
	        {
		    *xc = *xd; *xd = v; --xd;
		    SORT_SWAP2(*yc,*yd); --yd;
	        }
	        --xc; --yc;
	    }
	    if (xb > xc) break;
            SORT_SWAP1(*xb,*xc);
            SORT_SWAP2(*yb,*yc);
	    ++xb; ++yb;
	    --xc; --yc;
        }
    
        a = xa - x0;
        ba = xb - xa;
        if (ba > a) s = a; else s = ba;
        for (xl = x0, xh = xb-s, yl = y0, yh = yb-s; s > 0; --s)
        {
	    *xl = *xh; *xh = v; ++xl; ++xh;
	    SORT_SWAP2(*yl,*yh); ++yl; ++yh;
        }
        d = xd - x0;
        dc = xd - xc;
        if (dc > nn-1-d) s = nn-1-d; else s = dc;
        for (xl = xb, xh = x0+(nn-s), yl = yb, yh = y0+(nn-s); s > 0; --s)
        {
	    *xh = *xl; *xl = v; ++xl; ++xh;
	    SORT_SWAP2(*yl,*yh); ++yl; ++yh;
        }

	if (ba > dc)
	{
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	}
	else
	{
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	}
    }
}
#endif

#if SORT_OF_SORT == 3
#ifndef SORT_TYPE2
 #error "SORT_TYPE2 must be defined before including sorttemplates.c"
#endif

SORT_FUNCTYPE
SORT_NAME(SORT_TYPE1 *x, SORT_TYPE2 *y, int n)
{
    int i,j;
    int a,d,ba,dc,s,nn;
    SORT_TYPE2 tmp2,v,v1,v2,v3;
    SORT_TYPE1 tmp1,*x0,*xa,*xb,*xc,*xd,*xh,*xl;
    struct {SORT_TYPE1 *addr; int len;} stack[40];
    int top;

    top = 0;
    if (n > 1)
    {
        stack[top].addr = x;
        stack[top].len = n;
        ++top;
    }

    while (top > 0)
    {
	--top;
	x0 = stack[top].addr;
	nn = stack[top].len;

        if (nn < SORT_MINPARTITION)
        {
            for (i = 1; i < nn; ++i)
            {
                tmp1 = x0[i];
		tmp2 = y[tmp1];
                for (j = i; y[x0[j-1]] > tmp2; )
                {
                    x0[j] = x0[j-1];
                    if (--j == 0) break;
                }
                x0[j] = tmp1;
            }
	    continue;
        }

        if (nn < SORT_MINMEDIAN9)
            v = SORT_MEDIAN_OF_3(y[x0[0]],y[x0[nn/2]],y[x0[nn-1]]);
        else
        {
	    v1 = SORT_MEDIAN_OF_3(y[x0[0]],y[x0[1]],y[x0[2]]);
	    v2 = SORT_MEDIAN_OF_3(y[x0[nn/2-1]],y[x0[nn/2]],y[x0[nn/2+1]]);
	    v3 = SORT_MEDIAN_OF_3(y[x0[nn-3]],y[x0[nn-2]],y[x0[nn-1]]);
	    v = SORT_MEDIAN_OF_3(v1,v2,v3);
        }

        xa = xb = x0;  xc = xd = x0+(nn-1);
        for (;;)
        {
	    while (xb <= xc && y[*xb] <= v)
	    {
	        if (y[*xb] == v)
	        {
		    SORT_SWAP1(*xa,*xb); ++xa;
	        }
	        ++xb;
	    }
	    while (xc >= xb && y[*xc] >= v)
	    {
	        if (y[*xc] == v)
	        {
		    SORT_SWAP1(*xc,*xd); --xd;
	        }
	        --xc;
	    }
	    if (xb > xc) break;
            SORT_SWAP1(*xb,*xc);
	    ++xb;
	    --xc;
        }
    
        a = xa - x0;
        ba = xb - xa;
        if (ba > a) s = a; else s = ba;
        for (xl = x0, xh = xb-s; s > 0; --s)
        {
	    SORT_SWAP1(*xl,*xh); ++xl; ++xh;
        }
        d = xd - x0;
        dc = xd - xc;
        if (dc > nn-1-d) s = nn-1-d; else s = dc;
        for (xl = xb, xh = x0 + (nn-s); s > 0; --s)
        {
	    SORT_SWAP1(*xl,*xh); ++xl; ++xh;
        }

	if (ba > dc)
	{
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	}
	else
	{
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	}
    }
}
#endif

#if SORT_OF_SORT == 4
#ifndef SORT_COMPARE
 #error "SORT_COMPARE must be defined before including sorttemplates.c"
#endif

SORT_FUNCTYPE
SORT_NAME(SORT_TYPE1 *x, int n)
{
    int i,j;
    int a,d,ba,dc,s,nn;
    SORT_TYPE1 tmp2,v,v1,v2,v3;
    SORT_TYPE1 tmp1,*x0,*xa,*xb,*xc,*xd,*xh,*xl;
    struct {SORT_TYPE1 *addr; int len;} stack[40];
    int top;

    top = 0;
    if (n > 1)
    {
        stack[top].addr = x;
        stack[top].len = n;
        ++top;
    }

    while (top > 0)
    {
	--top;
	x0 = stack[top].addr;
	nn = stack[top].len;

        if (nn < SORT_MINPARTITION)
        {
            for (i = 1; i < nn; ++i)
            {
                tmp1 = x0[i];
                for (j = i; SORT_COMPARE(tmp1,x0[j-1]) < 0; )
                {
                    x0[j] = x0[j-1];
                    if (--j == 0) break;
                }
                x0[j] = tmp1;
            }
	    continue;
        }

        if (nn < SORT_MINMEDIAN9)
            v = F_SORT_MEDIAN_OF_3(SORT_COMPARE,x0[0],x0[nn/2],x0[nn-1]);
        else
        {
	    v1 = F_SORT_MEDIAN_OF_3(SORT_COMPARE,x0[0],x0[1],x0[2]);
	    v2 = F_SORT_MEDIAN_OF_3(SORT_COMPARE,x0[nn/2-1],x0[nn/2],x0[nn/2+1]);
	    v3 = F_SORT_MEDIAN_OF_3(SORT_COMPARE,x0[nn-3],x0[nn-2],x0[nn-1]);
	    v = F_SORT_MEDIAN_OF_3(SORT_COMPARE,v1,v2,v3);
        }

        xa = xb = x0;  xc = xd = x0+(nn-1);
        for (;;)
        {
	    while (xb <= xc && SORT_COMPARE(*xb,v) <= 0)
	    {
	        if (SORT_COMPARE(*xb,v) == 0)
	        {
		    SORT_SWAP1(*xa,*xb); ++xa;
	        }
	        ++xb;
	    }
	    while (xc >= xb && SORT_COMPARE(*xc,v) >= 0)
	    {
	        if (SORT_COMPARE(*xc,v) == 0)
	        {
		    SORT_SWAP1(*xc,*xd); --xd;
	        }
	        --xc;
	    }
	    if (xb > xc) break;
            SORT_SWAP1(*xb,*xc);
	    ++xb;
	    --xc;
        }
    
        a = xa - x0;
        ba = xb - xa;
        if (ba > a) s = a; else s = ba;
        for (xl = x0, xh = xb-s; s > 0; --s)
        {
	    SORT_SWAP1(*xl,*xh); ++xl; ++xh;
        }
        d = xd - x0;
        dc = xd - xc;
        if (dc > nn-1-d) s = nn-1-d; else s = dc;
        for (xl = xb, xh = x0 + (nn-s); s > 0; --s)
        {
	    SORT_SWAP1(*xl,*xh); ++xl; ++xh;
        }

	if (ba > dc)
	{
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	}
	else
	{
	    if (dc > 1)
	    {
	        stack[top].addr = x0+(nn-dc); stack[top].len = dc; ++top;
	    }
	    if (ba > 1)
	    {
	        stack[top].addr = x0; stack[top].len = ba; ++top;
	    }
	}
    }
}
#endif

#undef SORT_NAME
#undef SORT_OF_SORT
#undef SORT_TYPE1
#undef SORT_TYPE2
#undef SORT_COMPARE