#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

struct Photon {
	float pos[3];
};

int nswap;
#define	swap(ph, a, b) { Photon *ph2=ph[a]; ph[a]=ph[b]; ph[b]=ph2; ++nswap; }

void median_split(Photon** p, int start, int end, int median, int axis) {
	int left = start;
	int right = end;

	while (right > left) {
		const float v = p[(left+right)/2]->pos[axis];
		int i = left-1;
		int j = right;
		for (;;) {
			while (p[++i]->pos[axis] < v);
			while (p[--j]->pos[axis] > v);
			if (i >= j) break;
			swap(p, i, j);
		}

		if (i >= median)	right = i-1;
		if (i <= median)	left = i+1;
	}
}

void median_split_using_qsort(Photon** p, int start, int end, int axis) {
	int i, j;
	float x;

	x = p[(start + end) / 2]->pos[axis];
	i = start; j = end;
	for (;;) {
		while (p[i]->pos[axis] < x) ++i;
		while (x < p[j]->pos[axis]) --j;
		if (i >= j) break;
		Photon* t = p[i]; p[i] = p[j]; p[j] = t;
		++i; --j;
	}
	if (start < i-1) median_split_using_qsort(p, start, i-1, axis);
	if (j+1 < end)   median_split_using_qsort(p, j+1, end, axis);
}

void balance_segment(Photon** pbal, Photon** porg, int index, int start, int end) {
	//--------------------
	// compute new median
	//--------------------

	int median = 1;
	while ((4*median) <= (end-start+1))
		median += median;

	if ((3*median) <= (end-start+1))
		median = median * 2 + start-1;
	else
		median = end-median+1;

	//--------------------------
	// find axis to split along
	//--------------------------

	int axis = (rand() >> 8) % 3;

	//------------------------------------------
	// partition photon block around the median
	//------------------------------------------

	median_split(porg, start, end, median, axis);
//	median_split_using_qsort(porg, start, end, axis);

	pbal[index] = porg[median];

	//----------------------------------------------
	// recursively balance the left and right block
	//----------------------------------------------

	if (median > start) {
		// balance left segment
		if (start < median-1) {
			balance_segment(pbal, porg, 2*index, start, median-1);
		} else {
			pbal[2*index] = porg[start];
		}
	}

	if (median < end) {
		// balance right segment
		if (median+1 < end) {
			balance_segment(pbal, porg, 2*index+1, median+1, end);
		} else {
			pbal[2*index+1] = porg[end];
		}
	}
}

float randf() {
	return rand() / ((float)RAND_MAX+1);
}

#ifdef _MSC_VER
const double log_2 = log(2.0);
double log2(double x) {
	return log(x) / log_2;
}
#endif

int main(int argc, char* argv[]) {
	srand(time(NULL));
	int N = 100;

	if (argc > 1) {
		N = atoi(argv[1]);
	}

	Photon* buf = new Photon[N+1];
	Photon** porg = new Photon*[N+1];
	Photon** pbal = new Photon*[N+1];
	for (int i=1; i<=N; ++i) {
		Photon* p = &buf[i];
		p->pos[0] = randf();
		p->pos[1] = randf();
		p->pos[2] = randf();

		porg[i] = p;
	}

	balance_segment(pbal, porg, 1, 1, N);

	double r = N * log2((double)N) / 2;
	printf("#photon:%d, #swap:%d, n.log2(n)/2:%d\n", N, nswap, (int)r);

	delete[] buf;
	delete[] porg;
	delete[] pbal;
}

#photon:100, #swap:181, n.log2(n)/2:332