//////////////////////////////////////////////////////////////////////////
// Square Root Profiler
// Author: Shane Colliatie
// Date: 7/23/05
//////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <math.h>
using namespace std;

//#define UINT unsigned int
typedef unsigned int UINT;
#define ITERATIONS 10000000

//////////////////////////////////////////////////////////////////////////
// Assembly implementation of Square Root approximation
//////////////////////////////////////////////////////////////////////////
float FastSqrt (float f) {
	float RetVal;
	_asm{ mov eax, f
		sub eax, 0x3f800000
		sar eax, 1
		add eax, 0x3f800000
		mov RetVal, eax };
	return RetVal;
}

//////////////////////////////////////////////////////////////////////////
// System implementation of Square Root
//////////////////////////////////////////////////////////////////////////
float SlowSqrt (float f) {
	return sqrt(f);
}

//////////////////////////////////////////////////////////////////////////
// Quake 2 implementation of Square Root Approximation
//////////////////////////////////////////////////////////////////////////
float CarmSqrt(float x){
	union{
		int intPart;
		float floatPart;
	} convertor;
	union{
		int intPart;
		float floatPart;
	} convertor2;
	convertor.floatPart = x;
	convertor2.floatPart = x;
	convertor.intPart = 0x1FBCF800 + (convertor.intPart >> 1);
	convertor2.intPart = 0x5f3759df - (convertor2.intPart >> 1);
	return 0.5f*(convertor.floatPart + (x * convertor2.floatPart));
}

//////////////////////////////////////////////////////////////////////////
// Test function to calculate square roots
// Designed to ensure equality in all implementations
//
// Known issues:
//	Does not take initial cache misses into consideration
//////////////////////////////////////////////////////////////////////////
void testroot(float num_to_test, UINT NUM_ITERS)
{
	float answer = 0.0f;
	{
		__int64 FastStart, FastEnd, FastElapsed;
		// Get the time using assembly. VERY high precision
		_asm 
		{ 
			RDTSC
				mov DWORD PTR FastStart, eax
				mov DWORD PTR FastStart+4, edx 
		};

		UINT FastCalcs = 0;
		while( FastCalcs < NUM_ITERS)
		{
			FastSqrt(num_to_test);
			FastCalcs++;
		}
		answer = FastSqrt(num_to_test);
		_asm
		{ 
			RDTSC
				mov DWORD PTR FastEnd, eax
				mov DWORD PTR FastEnd+4, edx 
		};

		FastElapsed = FastEnd - FastStart;
		
		cout << "Fast:\t" << FastElapsed << "\t" << answer << '\n';
	}

	{
		__int64 SlowStart, SlowEnd, SlowElapsed;
		_asm 
		{ 
			RDTSC
				mov DWORD PTR SlowStart, eax
				mov DWORD PTR SlowStart+4, edx 
		};

		UINT SlowCalcs = 0;
		while( SlowCalcs < NUM_ITERS)
		{
			SlowSqrt(num_to_test);
			SlowCalcs++;
		}

		answer = SlowSqrt(num_to_test);

		_asm{ RDTSC
			mov DWORD PTR SlowEnd, eax
			mov DWORD PTR SlowEnd+4, edx 
		};

		SlowElapsed = SlowEnd - SlowStart;

		cout << "Slow:\t" << SlowElapsed << "\t" << answer << '\n';
	}

	{
		__int64 CarmStart, CarmEnd, CarmElapsed;

		_asm 
		{ 
			RDTSC
				mov DWORD PTR CarmStart, eax
				mov DWORD PTR CarmStart+4, edx 
		};

		UINT CarmCalcs = 0;
		while( CarmCalcs < NUM_ITERS)
		{
			CarmSqrt(num_to_test);
			CarmCalcs++;
		}

		answer = CarmSqrt(num_to_test);

		_asm
		{ 
			RDTSC
				mov DWORD PTR CarmEnd, eax
				mov DWORD PTR CarmEnd+4, edx 
		};

		CarmElapsed = CarmEnd - CarmStart;

		cout << "Carm:\t" << CarmElapsed << "\t" << answer << "\n\n";
	}
}

//////////////////////////////////////////////////////////////////////////
// Test function to calculate random square roots
// Designed to ensure equality in all implementations
//
// Known issues:
//	int to float typecast conversion has a small impact on clock cycles
//////////////////////////////////////////////////////////////////////////
void testrootrand(UINT NUM_ITERS)
{
	float answer = 0.0f;
	{
		__int64 FastStart, FastEnd, FastElapsed;
		_asm 
		{ 
			RDTSC
				mov DWORD PTR FastStart, eax
				mov DWORD PTR FastStart+4, edx 
		};

		UINT FastCalcs = 0;
		while( FastCalcs < NUM_ITERS)
		{
			FastSqrt((float)rand());
			FastCalcs++;
		}
		_asm
		{ 
			RDTSC
				mov DWORD PTR FastEnd, eax
				mov DWORD PTR FastEnd+4, edx 
		};

		FastElapsed = FastEnd - FastStart;

		cout << "Fast:\t" << FastElapsed << '\n';
	}

	{
		__int64 SlowStart, SlowEnd, SlowElapsed;
		_asm 
		{ 
			RDTSC
				mov DWORD PTR SlowStart, eax
				mov DWORD PTR SlowStart+4, edx 
		};

		UINT SlowCalcs = 0;
		while( SlowCalcs < NUM_ITERS)
		{
			SlowSqrt((float)rand());
			SlowCalcs++;
		}

		_asm{ RDTSC
			mov DWORD PTR SlowEnd, eax
			mov DWORD PTR SlowEnd+4, edx 
		};

		SlowElapsed = SlowEnd - SlowStart;

		cout << "Slow:\t" << SlowElapsed << '\n';
	}

	{
		__int64 CarmStart, CarmEnd, CarmElapsed;

		_asm 
		{ 
			RDTSC
				mov DWORD PTR CarmStart, eax
				mov DWORD PTR CarmStart+4, edx 
		};

		UINT CarmCalcs = 0;
		while( CarmCalcs < NUM_ITERS)
		{
			CarmSqrt((float)rand());
			CarmCalcs++;
		}

		_asm
		{ 
			RDTSC
				mov DWORD PTR CarmEnd, eax
				mov DWORD PTR CarmEnd+4, edx 
		};

		CarmElapsed = CarmEnd - CarmStart;

		cout << "Carm:\t" << CarmElapsed << "\n\n";
	}
}

int main()
{
	__int64 nTime;
	_asm
	{ 
		RDTSC
		mov DWORD PTR nTime, eax
		mov DWORD PTR nTime+4, edx 
	};

	srand((UINT)nTime);
	printf("Square Root Profiler:\t%i iterations\n\n", ITERATIONS);
	printf(" -root of 1000.0f\n\tTime\t\tResult\n");
	testroot(1000.0f, ITERATIONS);
	printf(" -root of 100.0f\n\tTime\t\tResult\n");
	testroot(100.0f, ITERATIONS);
	printf(" -root of 10.0f\n\tTime\t\tResult\n");
	testroot(10.0f, ITERATIONS);
	printf(" -root of 4.0f\n\tTime\t\tResult\n");
	testroot(4.0f, ITERATIONS);
	printf(" -random square root\n\tTime\n");
	testrootrand(ITERATIONS);

	return 0;
}

inline float SseSqrt(float f){	_mm_store_ss(&f, _mm_sqrt_ss(_mm_load_ss(&f)));	return f;}

#include <xmmintrin.h>

float sum = 0.0f;while( FastCalcs < NUM_ITERS){    sum += FastSqrt(num_to_test);    FastCalcs++;}printf("%f\n", sum); // printf is an external function and can't be optimized away, hence sum must be calculated

float sum = 0.0f;while( FastCalcs < NUM_ITERS){    sum += num_to_test;    FastCalcs++;}printf("%f\n", sum);