median

unsigned int select_double(double arr[], unsigned int len, unsigned int k);

void swap_double(double *a, double *b) {
	double tmp = *a;
	*a = *b;
	*b = tmp;
}

void median3_double(double *arr, unsigned int a, unsigned int b, unsigned int c) {
	if (arr[b] < arr[a]) {
		if (arr[b] < arr[c]) {
			if (arr[c] < arr[a]) {
				swap_double(&arr[b], &arr[c]);
			} else {
				swap_double(&arr[b], &arr[a]);
			}
		}
	} else if (arr[c] < arr[b]) {
		if (arr[c] < arr[a]) {
			swap_double(&arr[b], &arr[a]);
		} else {
			swap_double(&arr[b], &arr[c]);
		}
	}
}

// Hoare partition
unsigned int partition_double(double *arr, unsigned int len, unsigned int k) {
	swap_double(&arr[0], &arr[k]);
	unsigned int lo = 1;
	unsigned int hi = len - 1;

	for (;; ++lo, --hi) {
		for (;; ++lo) {
			if (lo > hi)
				goto loop_done;
			if (arr[lo] >= arr[0])
				break;
		}

		while (arr[0] < arr[hi])
			--hi;
		if (lo >= hi)
			break;

		swap_double(&arr[lo], &arr[hi]);
	}

loop_done:
	--lo;
	swap_double(&arr[0], &arr[lo]);
	return lo;
}

// repeated step
double pivot_double(double arr[], unsigned int len) {
	if (len < 9) {
		return partition_double(arr, len, len/2);
	}

	unsigned int j=0;
	for (unsigned int i=2; i<len; i+=3, j++) {
		median3_double(arr, i-2, i-1, i-0);
		swap_double(&arr[i-1], &arr[j]);
	}
	unsigned int m=0;
	for (unsigned int i=2; i<j; i+=3, m++) {
		median3_double(arr, i-2, i-1, i-0);
		swap_double(&arr[i-1], &arr[m]);
	}

	select_double(arr, m, m/2);
	return partition_double(arr, len, m/2);
}

// quickselect
unsigned int select_double(double arr[], unsigned int len, unsigned int k) {
	unsigned int p;
	unsigned int offset=0;
	while (1) {
		p = pivot_double(arr+offset, len);

		if (p == k) {
			return offset+p;
		}
		if (p > k) {
			len = p;
		} else {
			k = k - p - 1;
			offset = offset + (p + 1);
			len = len - (p + 1);
		}
	}
}

double median_double(double arr[], unsigned int len) {
	unsigned int m = select_double(arr, len, len/2);
	return arr[m];
}