简介
如果要求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)
网友评论