快速平方根倒数算法C程序(包括单精度与双精度)

程序摘自：A Brief History of InvSqrt (By Matthew Robertson）

/**********************************************************************************************
Function: fast_invSqrt32
Description: 快速平方根倒数
Input: 单精度浮点数number
Output: 无
Input_Output: 无
Return: 快速平方根倒数的值
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
float fast_invSqrt32(float number)
{
	long i;
	float x2, y;
	const float threehalfs = 1.5F;

	x2 = number * 0.5F;
	y = number;
	i = *(long*) &y;
	i = 0x5f375a86 - (i >> 1);
	y = *(float*) &i;
	y = y * (threehalfs - (x2 * y * y));
	y = y * (threehalfs - (x2 * y * y));

	return y;
}


/**********************************************************************************************
Function: fast_sqrt32
Description: 快速平方根
Input: 单精度浮点数number
Output: 无
Input_Output: 无
Return: 快速平方根的值
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
float fast_sqrt32(float number)
{
	return 1.0F / fast_invSqrt32(number);
}


/**********************************************************************************************
Function: fast_invSqrt64
Description: 快速平方根倒数
Input: 双精度浮点数number
Output: 无
Input_Output: 无
Return: 快速平方根倒数的值
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
double fast_invSqrt64(double number)
{
	long long i;
	double x2, y;
	const double threehalfs = 1.5;

	x2 = number * 0.5;
	y = number;
	i = *(long long*)& y;
	i = 0x5fe6eb50c7b537a9 - (i >> 1);
	y = *(double*) &i;
	y = y * (threehalfs - (x2 * y * y));
	y = y * (threehalfs - (x2 * y * y));

	return y;
}


/**********************************************************************************************
Function: fast_sqrt64
Description: 快速平方根
Input: 双精度浮点数number
Output: 无
Input_Output: 无
Return: 快速平方根的值
Author: Marc Pony(marc_pony@163.com)
***********************************************************************************************/
double fast_sqrt64(double number)
{
	return 1.0 / fast_invSqrt64(number);
}