myqsort/qsort.h

/* Quick sort in C */
#ifndef MY_QUICKSORT_H
#define MY_QUICKSORT_H

/* Used in quicksort_rand3_sp under which we call regular quicksort_rand3 without singlepassing */
#ifndef RAND3_SP_VAL
#define RAND3_SP_VAL 64
#endif

#include <stdint.h>

/* Uncomment for debugging with if(someevent) raise(SIGINT); */
#include <signal.h>

/* Structure:
 *
 * - BASICS:	Basic quicksort
 * - RANDOMPIV:	Randomized pivoting and partitioning by given pivot-indexed element.
 * - MINMAX:	partitioning while doing min-max over the whole + related extras.
 * - PMINMAX:	Partitioning with min-max searches in generated partitions + related.
 * - VALPART:	Partitioning by value - not by random or given pivot index!
 */

/* BASICS */

/** Swap operation */
static inline void swapit(uint32_t *a, uint32_t *b) {
	uint32_t t = *a;
	*a = *b;
	*b = t;
}

/**
 * Partition the array and find the pivot element such that
 *
 * - Elements smaller than pivot are on left of pivot
 * - Elements greater than pivot are on right of pivot
 *
 * @param array The array to partition
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @returns The partition point.
 */
static inline int partition(uint32_t array[], int low, int high) {
	/* This is "Lomuto"s unidirectional partitioner - see algorithms book */
	
	/* select the rightmost element as pivot */
	uint32_t pivot = array[high];

	/* index until smaller or eq elements lay */
	int i = (low - 1);

	/* traverse each element of the array */
	/* compare them with the pivot */
	#pragma GCC unroll 4
	for (int j = low; j < high; ++j) {
		if (array[j] <= pivot) {
			/* if element smaller than pivot is found */
			/* swap it with the greater element pointed by i */
			++i;

			/* swap element at i with element at j */
			swapit(&array[i], &array[j]);
		}
	}

	/* swap the pivot element with the greater element at i */
	swapit(&array[i + 1], &array[high]);

	/* return the partition point */
	return (i + 1);
}

/** Simple in-place recursive quicksort on array for elements in [low, high) indices */
static inline void quicksort(uint32_t array[], int low, int high) {
	if (low < high) {
		int pi = partition(array, low, high);

		/* recursive call on the left of pivot */
		quicksort(array, low, pi - 1);

		/* recursive call on the right of pivot */
		quicksort(array, pi + 1, high);
	}
}

/* THREEWAY */

struct pret3 {
	int leftend;
	int rightend;
};
typedef struct pret3 pret3;

/**
 * Partition the array threeway by puutting ALL pivots to the middle and smaller/biggers left-right
 *
 * @param array The array to partition
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param pivotval This value is used to partition the array.
 * @returns The partition points, end of left and right (inclusive)
 */
static inline pret3 partition3(uint32_t array[], int low, int high, uint32_t pivotval) {
	/* index until smaller or eq elements lay */
	int i = (low - 1);
	uint32_t pc = 0;

	/* traverse each element of the array */
	/* compare them with the pivotval */
	#pragma GCC unroll 4
	for (int j = low; j <= high; ++j) {
		/* Branchless pivot-count */
		pc += (array[j] == pivotval);

		if(array[j] < pivotval) {
			/* if element smaller than pivotval is found */
			/* swap it with the greater element pointed by i */
			++i;

			/* swap element at i with element at j */
			swapit(&array[i], &array[j]);
		}
	}

	/* Can spare out the second loop in these cases */
	if(pc < 2) {
		pret3 ret;
		ret.leftend = i;
		ret.rightend = i + 1;
		return ret;
	}

	/* index until smaller or eq elements lay */
	int i2 = (high + 1);

	#pragma GCC unroll 4
	for (int j = high; j > i; --j) {
		if (array[j] > pivotval) {
			/* if element smaller than pivot is found */
			/* swap it with the greater element pointed by i */
			--i2;

			/* swap element at i with element at j */
			swapit(&array[i2], &array[j]);
		}
	}

	/* return the partition points */
	pret3 ret;
	ret.leftend = i;
	ret.rightend = i2;
	return ret;
}

/**
 * Partition the array threeway by puutting ALL pivots to the middle and smaller/biggers left-right.
 * Single-pass version.
 *
 * @param array The array to partition
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param pivotval This value is used to partition the array.
 * @returns The partition points, end of left and right (inclusive)
 */
static inline pret3 partition3_sp(uint32_t array[], int low, int high, uint32_t pivotval) {

	/* Invariant for left: index until smaller (than pivot) elements lay */
	int il = (low - 1);
	/* Invariant for right: index until (from top) bigger elements lay */
	int ir = (high + 1);
	/* Indices from where we swap left and right into "is" (and sometimes swap among here too) */
	int jl = low;
	int jr = high;

	while(jl <= jr) {
		/* Jump through any pivot elements */
		while((array[jl] == pivotval) && (jl < high)) ++jl;
		while((array[jr] == pivotval) && (jr > low)) --jr;
		if(jl > jr) break;

		if(array[jl] < pivotval) {
			/* if element smaller than pivotval is found */
			/* swap it with the element pointed by i */
			++il;
			if(il != jl) swapit(&array[il], &array[jl]);
		} else {
			/* Here its sure a[jl] and a[jr] are not pivots */
			/* Which means array[jl] > pivotval and on left */

			if(array[jr] < pivotval) {
				/* We can totally safely just swap them */
				/* Because both are in wrong places.... */
				swapit(&array[jl], &array[jr]);
				/* swap it with the element pointed by i */
				++il;
				if(il != jl) swapit(&array[il], &array[jl]);
			} else {
				/* Left->right move needed, but right is */
				/* at good place already so handle that! */
				/* ensured: if(array[jr] > pivotval) */

				/* Progress invariants */
				/* swap with the element pointed by i */
				--ir;
				if(ir != jr) swapit(&array[ir], &array[jr]);
				/* Because of swap we might swaped in pivot! */
				/* So continue our loop but do not step jl! */
				/* ++jl; // incomplete so no step! */
				--jr; /* STEP RIGHT because "continue"! */
				continue;
			}
		}
		/* Step left being processed */
		++jl;

		if(array[jr] > pivotval) {
			/* if element bigger than pivotval is found */
			/* swap it with the element pointed by i */
			--ir;
			if(ir != jr) swapit(&array[ir], &array[jr]);
		} else {
			/* On wrong side - so we just need to look */
			/* for left-wrong-sided to swap it with! */
			continue; /* I think no ++jl needed! */
		}
		/* Step right being processed */
		--jr;
	}

	/* return the partition points */
	pret3 ret;
	ret.leftend = il;
	ret.rightend = ir;
	return ret;
}

/**
 * Partition the array threeway by puutting ALL pivots to the middle and smaller/biggers left-right.
 * Single-pass version 2.0
 *
 * @param array The array to partition
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param pivotval This value is used to partition the array.
 * @returns The partition points, end of left and right (inclusive)
 */
static inline pret3 partition3_sp2(uint32_t array[], int low, int high, uint32_t pivotval) {

	/* Invariant for left: index until smaller (than pivot) elements lay */
	int il = (low - 1);
	/* Invariant for right: index until (from top) bigger elements lay */
	int ir = (high + 1);
	/* Indices from where we swap left and right into "is" (and sometimes swap among here too) */
	int jl = low;
	int jr = high;

	while(jl <= jr) {
		/* Handle left and find wrongly placed element */
		while((array[jl] <= pivotval) && (jl <= jr)) {
			int isNonPivot = (array[jl] != pivotval);
			int nonSameIndex = (il + 1 != jl);
			if(isNonPivot & nonSameIndex)
				swapit(&array[il + 1], &array[jl]);
			il += isNonPivot;
			++jl;
		}

		/* Handle right and find wrongly placed element */
		while((array[jr] >= pivotval) && (jl <= jr)) {
			int isNonPivot = (array[jr] != pivotval);
			int nonSameIndex = (ir - 1 != jr);
			if(isNonPivot & nonSameIndex)
					swapit(&array[ir - 1], &array[jr]);
			ir -= isNonPivot;
			--jr;
		}

		/* Swap the two found elements that are wrongly placed */
		if(jl < jr) swapit(&array[jl], &array[jr]);
	}

	/* return the partition points */
	pret3 ret;
	ret.leftend = il;
	ret.rightend = ir;
	return ret;
}

/* RANDOMPIV */

/**
 * Partition the array and using the pivot index
 *
 * - Elements smaller than pivot are on left of pivot
 * - Elements greater than pivot are on right of pivot
 *
 * @param array The array to partition
 * @param pi The index of the pivot element to use. 0 or high is what OG quicksorts do.
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @returns The partition point.
 */
static inline int partition_with_pivot(uint32_t array[], int pi, int low, int high) {
	/*
	 * Rem.: This looks like overhead,
	 * but after seriously considering
	 * writing the whole out I can tell
	 * this is still fastests basically.
	 */

	/* swap pivot with rightmost */
	swapit(&array[high], &array[pi]);
	/* delegate to previous sol. */
	return partition(array, low, high);
}

// 32-bit LCG for fast random generations
static inline uint32_t lcg(uint32_t *state) {
	*state = *state * 1664525u + 1013904223u;
	return *state;
}

/** Get pivot index in [0, len-1] without modulus - see our fastrand.h */
static inline uint32_t pick_pivot(uint32_t *state, uint32_t len) {
	uint32_t rand = lcg(state);
	/* Multiply by len, take the upper 32 bits of the 64-bit result */
	return (uint32_t)(((uint64_t)rand * len) >> 32);
}

typedef uint32_t rpivotstate;

/** Randomized pivoting in-place recursive quicksort on array for elements in [low, high] indices */
static inline void quicksort_rand(uint32_t array[], int low, int high, rpivotstate *state) {
	if (low < high) {
		int pi = pick_pivot(state, (high + 1) - low) + low;
		pi = partition_with_pivot(array, pi, low, high);

		/* recursive call on the left of pivot */
		quicksort_rand(array, low, pi - 1, state);

		/* recursive call on the right of pivot */
		quicksort_rand(array, pi + 1, high, state);
	}
}

/* THREEWAY_RAND */

/** Randomized pivoting in-place recursive quicksort on array for elements in [low, high] indices */
static inline void quicksort_rand3(uint32_t array[], int low, int high, rpivotstate *state) {
	if (low < high) {
		/* partition threeway by random pivot */
		int pi = pick_pivot(state, (high + 1) - low) + low;
		pret3 res = partition3(array, low, high, array[pi]);

		/* recursive call on the left of pivot */
		quicksort_rand3(array, low, res.leftend, state);

		/* recursive call on the right of pivot */
		quicksort_rand3(array, res.rightend, high, state);
	}
}

/** Randomized pivoting in-place recursive quicksort on array for elements in [low, high] indices */
static inline void quicksort_rand3_sp(uint32_t array[], int low, int high, rpivotstate *state) {
	if (low < high) {
		/* partition threeway by random pivot */
		int pi = pick_pivot(state, (high + 1) - low) + low;
		pret3 res;
		if(high - low > RAND3_SP_VAL) res = partition3_sp(array, low, high, array[pi]);
		else res = partition3(array, low, high, array[pi]);

		/* recursive call on the left of pivot */
		quicksort_rand3_sp(array, low, res.leftend, state);

		/* recursive call on the right of pivot */
		quicksort_rand3_sp(array, res.rightend, high, state);
	}
}

/** Randomized pivoting in-place recursive quicksort on array for elements in [low, high] indices */
static inline void quicksort_rand3_sp2(uint32_t array[], int low, int high, rpivotstate *state) {
	if (low < high) {
		/* partition threeway by random pivot */
		int pi = pick_pivot(state, (high + 1) - low) + low;
		pret3 res;
		res = partition3_sp2(array, low, high, array[pi]);
		/* if(high - low > RAND3_SP_VAL) res = partition3_sp2(array, low, high, array[pi]);
		else res = partition3(array, low, high, array[pi]); */

		/* recursive call on the left of pivot */
		quicksort_rand3_sp2(array, low, res.leftend, state);

		/* recursive call on the right of pivot */
		quicksort_rand3_sp2(array, res.rightend, high, state);
	}
}

/* MINMAX */

/**
 * Partition the array while doing a min-max search and find the pivot element such that
 *
 * - Elements smaller than pivot are on left of pivot
 * - Elements greater than pivot are on right of pivot
 *
 * @param array The array to partition
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param minout OUT: Will be filled with the minimum key
 * @param maxout OUT: Will be filled with the maximum key
 * @returns The partition point.
 */
static inline int partition_and_minmax(uint32_t array[], int low, int high, uint32_t *minout, uint32_t *maxout) {
	/* This is "Lomuto"s unidirectional partitioner - see algorithms book */
	
	/* select the rightmost element as pivot */
	uint32_t pivot = array[high];
	*minout = pivot;
	*maxout = pivot;
	
	/* index until smaller or eq elements lay */
	int i = (low - 1);

	/* traverse each element of the array */
	/* compare them with the pivot */
	#pragma GCC unroll 4
	for (int j = low; j < high; ++j) {
		/* Branchless min-max */
		*minout = array[j] < *minout ? array[j] : *minout;
		*maxout = array[j] > *maxout ? array[j] : *maxout;

		/* Lomuto partitioning */
		if (array[j] <= pivot) {
			/* if element smaller than pivot is found */
			/* swap it with the greater element pointed by i */
			++i;

			/* swap element at i with element at j */
			swapit(&array[i], &array[j]);
		}
	}

	/* swap the pivot element with the greater element at i */
	swapit(&array[i + 1], &array[high]);

	/* return the partition point */
	return (i + 1);
}

/**
 * Partition the array and min-max and using the pivot index
 *
 * - Elements smaller than pivot are on left of pivot
 * - Elements greater than pivot are on right of pivot
 *
 * @param array The array to partition
 * @param pi The index of the pivot element to use. 0 or high is what OG quicksorts do.
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param minout OUT: Will be filled with the minimum key
 * @param maxout OUT: Will be filled with the maximum key
 * @returns The partition point.
 */
static inline int partition_and_minmax_with_pivot(uint32_t array[], int pi, int low, int high, uint32_t *minout, uint32_t *maxout) {
	/*
	 * Rem.: This looks like overhead,
	 * but after seriously considering
	 * writing the whole out I can tell
	 * this is still fastests basically.
	 */

	/* swap pivot with rightmost */
	swapit(&array[high], &array[pi]);
	/* delegate to previous sol. */
	return partition_and_minmax(array, low, high, minout, maxout);
}

/* PMINMAX */

/**
 * Partition the array with partition-based min-max search (4 values: 2 per partition) and find the pivot element such that
 *
 * - Elements smaller than pivot are on left of pivot
 * - Elements greater than pivot are on right of pivot
 *
 * (***): The left min-max outputs can return not just the left partition, but left partition ++ pivotpoint region min-max!
 *        This does not only happen if partition is empty, but also when pivot was already highest and last earlier.
 *        Would be hard to handle this edge-case and min-max out is generally used as output hints only so I prefer speed..
 *
 * @param array The array to partition
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param minout_left OUT: Will be filled with the minimum key of left partition or pivot value when partition empty. Also (***)
 * @param maxout_left OUT: Will be filled with the maximum key of left partition or pivot value when partition empty. Also (***)
 * @param minout_right OUT: Will be filled with the minimum key of right partition or pivot value when partition empty
 * @param maxout_right OUT: Will be filled with the maximum key of right partition or pivot value when partition empty
 * @returns The partition point.
 */
static inline int partition_and_pminmax(
		uint32_t array[],
		int low,
		int high,
		uint32_t *minout_left,
		uint32_t *maxout_left,
		uint32_t *minout_right,
		uint32_t *maxout_right) {
	/* This is "Lomuto"s unidirectional partitioner - see algorithms book */

	/* select the rightmost element as pivot */
	uint32_t pivot = array[high];
	*minout_left = pivot;
	*maxout_left = pivot;
	*minout_right = pivot;
	*maxout_right = pivot;

	/* index until smaller or eq elements lay */
	int i = (low - 1);

	/* traverse each element of the array */
	/* compare them with the pivot */
	#pragma GCC unroll 4
	for (int j = low; j < high; ++j) {
		/* Lomuto partitioning */
		if (array[j] <= pivot) {
			/* Branchless min-max */
			*minout_left = array[j] < *minout_left ? array[j] : *minout_left;
			*maxout_left = array[j] > *maxout_left ? array[j] : *maxout_left;

			/* if element smaller than pivot is found */
			/* swap it with the greater element pointed by i */
			++i;

			/* swap element at i with element at j */
			swapit(&array[i], &array[j]);
		} else {
			/* Branchless min-max */
			*minout_right = array[j] < *minout_right ? array[j] : *minout_right;
			*maxout_right = array[j] > *maxout_right ? array[j] : *maxout_right;
		}
	}

	/* swap the pivot element with the greater element at i */
	swapit(&array[i + 1], &array[high]);

	/* return the partition point */
	return (i + 1);
}

/**
 * Partition the array with partition-based min-max search (4 values: 2 per partition) and using the pivot index
 *
 * - Elements smaller than pivot are on left of pivot
 * - Elements greater than pivot are on right of pivot
 *
 * @param array The array to partition
 * @param pi The index of the pivot element to use. 0 or high is what OG quicksorts do.
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param minout OUT: Will be filled with the minimum key
 * @param maxout OUT: Will be filled with the maximum key
 * @returns The partition point.
 */
static inline int partition_and_pminmax_with_pivot(
		uint32_t array[],
		int pi,
		int low,
		int high,
		uint32_t *minout_left,
		uint32_t *maxout_left,
		uint32_t *minout_right,
		uint32_t *maxout_right) {
	/*
	 * Rem.: This looks like overhead,
	 * but after seriously considering
	 * writing the whole out I can tell
	 * this is still fastests basically.
	 */

	/* swap pivot with rightmost */
	swapit(&array[high], &array[pi]);
	/* delegate to previous sol. */
	return partition_and_pminmax(array, low, high, minout_left, maxout_left, minout_right, maxout_right);
}

/* VALPART */

/**
 * Partition the array using pivot value - and find pivot closest to that value (and place them at proper pivot index)
 *
 * - Elements smaller-eq than pivotval are on left of pivot
 * - Elements greater than pivotval are on right of pivot
 * - The "pivot" element we find is the biggest among the ones smaller-eq to the pivot value
 *   or if there is no such (all is greater) we return first greater-than value index (right[0])
 *
 * @param array The array to partition
 * @param low From when. (inclusive)
 * @param high Until when. (inclusive too!)
 * @param pivotval This value is used to partition the array.
 * @returns The partition point.
 */
static inline int partition_with_pivotval(uint32_t array[], int low, int high, uint32_t pivotval) {
	/* This is "Lomuto"s unidirectional partitioner - see algorithms book */

	/* Select the rightmost element as pivot just because */
	/* Need some start-value for min(abs(pv - p)) search! */
	int64_t leftmax = -1;
	/* Index of currently found pivot value */
	uint32_t pivoti = low;

	/* index until smaller or eq elements lay */
	int i = (low - 1);

	/* traverse each element of the array */
	/* compare them with the pivotval */
	/* The "<=" is needed for our trickz here too */
	#pragma GCC unroll 4
	for (int j = low; j <= high; ++j) {
		/* This "<=" ensures pivoti must be only searched among "left" values! */
		if (array[j] <= pivotval) {
			/* if element smaller than pivot is found */
			/* swap it with the greater element pointed by i */
			++i;

			/* swap element at i with element at j */
			swapit(&array[i], &array[j]);

			/* After this, array[i] can never change - so we can save it as a found pivot-index */
			/* Max-search on elements by telling which is closest to pivotval by abs difference! */
			if(array[i] > leftmax) {
				pivoti = i;
				leftmax = array[i];
			}
		}
	}

	/* Must check if above loop found elem at all - because its guessing */
	if(i != (low - 1)) {
		/* swap the pivot element into its place */
		swapit(&array[i], &array[pivoti]);
	}

	/* return the partition point: index of pivot element */
	return i;
}

#endif /* MY_QUICKSORT_H */