//Pigeon sort
//see Nov 1988 Australian Personal Computer (Vol 9 No 11) p233
//Tue 22 Feb 1999
//Why don't they just use radix sort on the no. of bits in a double?
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

enum {
	p_size = 10000,
	nmax = 10024+1,
	nocommaontheendofthisone
};

void pigeon_sort(double *p, int n, int *ip) {
	static int ipp[nmax];
	static int no_in_pig[nmax]; //should start from 1, I'll ignore 0
	static int no_of_slot[nmax];
	static int no_in_slot[nmax]; //these two have to start from zero
	static int npig[nmax+1];
	double pmin, pmax, slot;

	int i, istart, istart1, istart2, istart3,
		iend, iptemp, ip0, ip1, ip2, ip3, j;
	int jtemp;
	
	for (i=0; i<=nmax; i++) {
		ipp[i] = 0;
		no_in_pig[i] = 0;
		no_of_slot[i] = 0;
		no_in_slot[i] = 0;
		npig[i] = 0;
	}

	istart = 1;
	no_in_pig[1] = n;

	while (istart < n) {
		while (no_in_pig[istart] > 1) {
			if (no_in_pig[istart] > 2) {
				if (no_in_pig[istart] > 3) {
					if (no_in_pig[istart] > 4) {
						iend = istart + no_in_pig[istart] - 1;
						istart1 = istart + 1;
						pmax = p[ip[istart]];
						pmin = pmax;
						for (i=istart1; i<=iend; i++) {
							if (p[ip[i]] > pmax) {
								pmax = p[ip[i]];
							} else if (p[ip[i]] < pmin) {
								pmin = p[ip[i]];
							}
						}
						if (pmax == pmin) goto areequal;
						slot = ((double) no_in_pig[istart])/(pmax-pmin);
						for (i=istart; i<=iend; i++) {
							no_of_slot[i] = (int) rint((p[ip[i]] - pmin)*slot);
							no_in_slot[no_of_slot[i]]++;
						}
						npig[0] = istart;
						jtemp = no_in_pig[istart];
						for (j=0; j<=jtemp; j++) {
							npig[j+1] = npig[j] + no_in_slot[j];
							if (no_in_slot[j]) {
								no_in_pig[npig[j]] = no_in_slot[j];
							}
							no_in_slot[j] = 0;
						}
						for (i=istart; i<=iend; i++) {
							ipp[npig[no_of_slot[i]]] = ip[i];
							npig[no_of_slot[i]]++;
						}
						for (i=istart; i<=iend; i++) ip[i] = ipp[i];
					} else {
						// sort 4
						istart1 = istart + 1;
						istart2 = istart + 2;
						istart3 = istart + 3;
						if (p[ip[istart]] > p[ip[istart2]]) {
							iptemp = ip[istart2];
							ip[istart2] = ip[istart];
							ip[istart] = iptemp;
						}
						if (p[ip[istart1]] > p[ip[istart3]]) {
							iptemp = ip[istart3];
							ip[istart3] = ip[istart1];
							ip[istart1] = iptemp;
						}
						if (p[ip[istart]] > p[ip[istart1]]) {
							iptemp = ip[istart1];
							ip[istart1] = ip[istart];
							ip[istart] = iptemp;
						}
						if (p[ip[istart2]] > p[ip[istart3]]) {
							iptemp = ip[istart3];
							ip[istart3] = ip[istart2];
							ip[istart2] = iptemp;
						}
						if (p[ip[istart1]] > p[ip[istart2]]) {
							iptemp = ip[istart2];
							ip[istart2] = ip[istart1];
							ip[istart1] = iptemp;
						}
						istart = istart3;
						no_in_pig[istart] = 1;
					}
				} else {
					//sort 3
					istart1 = istart + 1;
					istart2 = istart + 2;
					ip0 = ip[istart];
					ip1 = ip[istart1];
					ip2 = ip[istart2];
					if (p[ip0] > p[ip1]) {
						if (p[ip2] > p[ip1]) {
							if (p[ip0] > p[ip2]) {
								ip[istart1] = ip2;
								ip[istart2] = ip0;
								ip[istart] = ip1;
							} else {
								ip[istart1] = ip0;
								ip[istart] = ip1;
							}
						} else {
							ip[istart2] = ip0;
							ip[istart] = ip2;
						}
					} else {
						if (p[ip1] > p[ip2]) {
							if (p[ip0] > p[ip2]) {
								ip[istart1] = ip0;
								ip[istart] = ip2;
								ip[istart2] = ip1;
							} else {
								ip[istart1] = ip2;
								ip[istart2] = ip1;
							}
						}
					}
					istart = istart2;
					no_in_pig[istart] = 1;
				}
			} else {
				//sort 2
				istart1 = istart + 1;
				if (p[ip[istart]] > p[ip[istart1]]) {
					iptemp = ip[istart];
					ip[istart] = ip[istart1];
					ip[istart1] = iptemp;
				}
				istart = istart1;
				no_in_pig[istart1] = 1;
			}
		}
areequal:
		istart=istart+no_in_pig[istart];
	}
}

int main(void) {
	double p[p_size+1];
	int ip[p_size+1];
	int i;

	for (i=1; i<=p_size; i++) {
		p[i] = ((double) (rand() % 1000000))/1000;
		ip[i] = i;
	}

	for (i=1; i<=p_size; i++) {
		printf("%lf\n", p[i]);
	}

	fprintf(stderr, "Sorting...\n");
	pigeon_sort(p, p_size, ip);
	fprintf(stderr, "...done\n");

	for (i=1; i<=p_size; i++) {
		printf("%lf\n", p[ip[i]]);
	}
}