快速平方根导数算法

简介

如果要求x的平方根倒数 , 一般我们都会写成这个形式1/sqrt(x)，里面包含一个除法操作和一个sqrt操作

针对于基本的加减乘除计算来说，而除法相对于计算机来说是最为复杂的运算

在某些场景下，大量的除法运算显然是不能满足需求的

快速平方根算法

float Q_rsqrt(float number)

{

 int32_t i; 

 float x2, y; 

 const float threehalfs = 1.5f;

 x2 = number * 0.5f;

 y = number;

 i = * (int32_t *) &y;

 i = 0x5f3759df - (i >> 1);

 y = * (float *) &i;

 y = y * (threehalfs - (x2 * y * y));

 //y = y * (threehalfs - (x2 * y * y));

 //y = y * (threehalfs - (x2 * y * y));

 return y;

}//来源于雷神之锤开源代码

浮点数在计算机内表示的方式（IEEE754）

S EEE EEEE E MMM MMMM MMMM MMMM MMMM MMMM

image.png

sign exp tails

二进制只有0和1，以科学记数法的形式表示二进制浮点数的时候，为了提升精度，首位的1会被省略掉

浮点数的二进制表示：

image.png

而浮点数实际的数值表示为如下公式

image.png image.png image.png image.png image.png

某种程度上来说，浮点数的二进制表示即为自身的对数 float_x ~= log(float_x)

c/c++ 数据类型，二进制移位

移位操作相比除法是很快的

110=6 11 = 110>>1 = 3 1100=110<<1=12

111=7 11 = 111>>1 = 3. 1110 = 111<<1 = 14

内存中存储的二进制数据是固定的，类型描述的是描述数据的元数据

如果普通的类型转换,会损失小数，也会打乱原来的二进制形式
i = (int32_t) y y(3.33) -> i(3)

如果我们想保留原二进制形式，就需要对描述该二进制的类型进行转换，而不是对应的值进行转换（有点像union）

i = * (int32_t *) &y; y(3.33) -> i(1079320248) 

某种程度上来说，浮点数的二进制表示即为自身的对数 float_x ~= log(float_x)

i = log(y)

如果要求1/sqrt(y) 就可以转换为 求 log(1/sqrt(y)) 就可以转换为 - 1/2 log(y) 就可以转换为求 -(i>>1)

0x5f3759df 这个近似项的计算：

假设w = 1/sqrt(y) =>

log(w) = -1/2 log(y) => 设

log(w_F) = -1/2 log(y_F) =>

w_F + u - 127 = -1/2(y_F + u - 127) =>

w_F = =>

w_F = 0x5f3759df - (i>>1)

牛顿迭代法

求f(x) = 0的解 x_n = x - dx

dx = (dy / dx) / (dy) = f'(x) / f(x)

x_n = x - f(x) / f'(x)

[图片上传失败...(image-65c956-1698394186269)]

这里的目标公式应该为 f(y)=1/y^2 - x

y_n =y - f(y) / f'(y)

=

=

y = y * (threehalfs - (x2 * y * y));

测试代码

#include<stdio.h>

#include<iostream>

#include<cmath>

#include<ctime>

using namespace std;

float Q_rsqrt(float number);

int main(){

 float a = 1.4444444444444444444444444444f;

 float ret, ret2;

 clock_t start,end;

 start = clock(); 

 for(int i=0; i <=100000000; i++) {

 ret = Q_rsqrt(a);

 }

 end = clock();

 cout<<"time_1: "<<(end - start)<<endl;

 start = clock(); 

 for(int i=0; i <=1000000000; i++)  {

 ret2 = 1/sqrt(a);

 }

 end = clock();

 cout<<"time_2: "<<(end - start)<<endl;

 std::cout<<"ret :" <<ret<<std::endl;

 std::cout<<"c++ 1/sqrt:" << ret2<<std::endl;

 return 0;

}

float Q_rsqrt(float number)

{

 int32_t i; //00000000 00000000 00000000

 //IEEE 754   //1       8(-127~128)   (23)

 float x2, y; //0(sign) 00000000(exp) 000000000000000(tails)

 const float threehalfs = 1.5f;

 x2 = number * 0.5f;

 y = number;

 i = * (int32_t *) &y;

 i = 0x5f3759df - (i >> 1);

 y = * (float *) &i;

 y = y * (threehalfs - (x2 * y * y));

 y = y * (threehalfs - (x2 * y * y));

 y = y * (threehalfs - (x2 * y * y));

 return y;

}

测试结果

编译器 gcc8.2

0次迭代

time_1: 599641

time_2: 5312553

ret :0.855104

c++ 1/sqrt:0.83205

1次迭代

time_1: 916537

time_2: 5310806

ret :0.831083

c++ 1/sqrt:0.83205

2次迭代

time_1: 1111124

time_2: 5208424

ret :0.832049

c++ 1/sqrt:0.83205

3次迭代

time_1: 1463755

time_2: 5279041

ret :0.83205

c++ 1/sqrt:0.83205

其他

eigen的rsqrt用的就是标准库的1/ std::sqrt(a)