排序浮点数
我有一个法线向量的列表，我正在计算标量三乘积，并对它们进行排序。我比较了三种情况下的分类：
- 用Matlab排序求最大绝对三重乘积值
- 使用
std::sort函数在C++中得到std::vector的产品值。三重产品使用双倍。 - 在OpenCL C中使用基排序-将绝对浮点数转换为无符号整数并将它们转换回来。我正在使用
cl_float的三重产品。
它们都给出了不同的值以及不同的指标，从而导致了算法中的问题。在这种情况下有什么问题，我怎样才能保持它们的一致性？
回答 1
Stack Overflow用户
发布于 2014-04-07 19:51:47
眼前的问题是：
- 计算三维顶点的标量三乘积，即表示为
binary32 float的每个向量的每个分量. - 能够判断该计算的结果是否大于该计算的其他结果。
到目前为止还不错，但是如果我们直接将公式应用到向量上，可能会在运算中丢失一些位，我们将无法分辨出两个结果。正如@rcgldr所指出的，排序不是问题，而是精度问题。
浮点舍入问题的一个解决方案是增加位数，即使用double。您说您没有double，所以让我们自己来做:只要我们需要，就在一个unsigned char数组中执行整个计算。
好吧，好吧，实际考虑：
- 输入由归一化向量组成，因此长度不大于1，这意味着没有任何分量大于1。
binary32浮点的指数从-127 (零，非正态)到127 (或无穷大的128个)，但所有分量的指数从-127到0(否则它们会大于1)。- 输入的最大精度为24位。
- 标量三积涉及向量积和标量积。在向量乘积(首先发生)中有乘法结果的减法，在标量积中有乘法结果的和。
注意事项2和3告诉我们，整个输入组件家族都可以采用大小为127位的定点格式进行偏移，而尾数为24位，即19字节。让我们定在24点。
为了能够完全表示所有可能的和减，一个额外的位就足够了(在结转的出现)，但是要完全表示所有可能的乘法结果，我们需要双倍的位数，所以它决定，加倍的大小就足以表示向量乘法，三倍它将使它足够用于标量积中的下一个乘法。
下面是一个类的草稿，它将浮点数加载到非常大的不动点格式，将符号保留为bool标志(有一个帮助函数rollArrayRight()，我将单独发布，但希望它的名称解释它)：
const size_t initialSize=24;
const size_t sizeForMult1=initialSize+initialSize;
const size_t finalSize=sizeForMult1+initialSize;
class imSoHuge{
public:
bool isNegative;
uint8_t r[finalSize];
void load(float v){
isNegative=false;
for(size_t p=0;p<finalSize;p++)r[p]=0;
union{
uint8_t b[4];
uint32_t u;
float f;
} reunion;
reunion.f=v;
if((reunion.b[3]&0x80) != 0x00)isNegative=true;
uint32_t m, eu;
eu=reunion.u<<1; //get rid of the sign;
eu>>=24;
m=reunion.u&(0x007fffff);
if(eu==0){//zero or denormal
if(m==0)return; //zero
}else{
m|=(0x00800000); //implicit leading one if it's not denormal
}
int32_t e=(int32_t)eu-127; //exponent is now in [e]. Debiased (does this word exists?)
reunion.u=m;
r[finalSize-1]=reunion.b[3];
r[finalSize-2]=reunion.b[2];
r[finalSize-3]=reunion.b[1];
r[finalSize-4]=reunion.b[0];
rollArrayRight(r, finalSize, e-(sizeForMult1*8)); //correct position for fixed-point
}
explicit imSoHuge(float v){
load(v);
}
};当类是用数字1.0f构造的，例如，数组r有类似于00 00 00 00 80 00的东西，注意它被加载到它的下部，乘法将相应地将这个数字滚动到上字节，然后我们可以恢复我们的浮点数。
为了使它有用，我们需要实现与和乘法的等价性，只要我们记住，我们只能将相同次数的和数组相乘(就像在三重乘积中)，否则它们的大小就不匹配了。
举个例子，这样的类会产生不同的效果：
考虑以下三个向量：
float a[]={0.0097905760, 0.0223784577, 0.9997016787};
float b[]={0.8248013854, 0.4413521587, 0.3534274995};
float c[]={0.4152690768, 0.3959976136, 0.8189856410};下面是计算三重乘积的函数：(希望我做得对，哈哈)
float fTripleProduct(float*a, float*b, float*c){
float crossAB[3];
crossAB[0]=(a[1]*b[2])-(a[2]*b[1]);
crossAB[1]=(a[2]*b[0])-(a[0]*b[2]);
crossAB[2]=(a[0]*b[1])-(a[1]*b[0]);
float tripleP=(crossAB[0]*c[0])+(crossAB[1]*c[1])+(crossAB[2]*c[2]);
return tripleP;
}fTripleProduct(a,b,c);的结果是0.1336331
如果我们将a的fisrt组件的最后一个数字从0改为6，使其成为0.0097905766 (它具有不同的十六进制表示)并再次调用该函数，结果是相同的，但我们知道它应该更大。
现在考虑一下，我们已经为imSoHuge类实现了乘法、求和和减法，并有了一个函数来使用它计算三重乘积。
imSoHuge tripleProduct(float*a, float*b, float*c){
imSoHuge crossAB[3];
crossAB[0]=(imSoHuge(a[1])*imSoHuge(b[2]))-(imSoHuge(a[2])*imSoHuge(b[1]));
crossAB[1]=(imSoHuge(a[2])*imSoHuge(b[0]))-(imSoHuge(a[0])*imSoHuge(b[2]));
crossAB[2]=(imSoHuge(a[0])*imSoHuge(b[1]))-(imSoHuge(a[1])*imSoHuge(b[0]));
imSoHuge tripleP=(crossAB[0]*imSoHuge(c[0]))+(crossAB[1]*imSoHuge(c[1]))+(crossAB[2]*imSoHuge(c[2]));
return tripleP;
}如果我们为向量的上述两个版本调用该函数，则数组中的结果不同：
0 0 0 4 46 b9 4 69 39 3f 53 b8 19 e0 ...
0 0 0 4 46 b9 4 85 93 82 df ba 7d 80 ...实际上，它们在binary32 float的精度之后不同，这意味着如果我们将数组转换回浮点数，它将是相同的浮点数，但是如果我们比较数组，我们就可以知道哪个数组更大。
将这个推理应用到测试中，我已经做了一个完整的示例，您可以立即用GCC中的-O3 -Wall -std=c++11编译和运行，或者在另一个编译器上运行类似的代码，然后输出：
Using class: second result is greater
casting to float:
first reasult: 1.336331e-001
second result: 1.336331e-001
as floats, the results are the same: 1.336331e-001源代码在这里(在Ideone上运行得很好)：
C++11码的源代码
如果您还没有迁移到C++11，那么如果您自己定义了精确的宽度类型( uint8_t、uint16_t、uint32_t、int32_t )，那么代码将在C++98中编译和运行。
怎么用？
简单地使用输入调用函数tripleProduct，并使用提供的重载比较器操作符比较结果，还可以使用提供的重载转换操作符将类imSoHuge转换为浮动(在计算了三元乘积之后)。
您可以为任何排序算法提供该类的数组和比较器。
结论和考虑：
注意，浮点乘法现在被执行为两个70+字节长数组的乘法，这意味着需要数百倍的时钟周期，加上和、比较等，这将是慢几千倍的，但是它是准确的。
整个算法的设计是使用归一化向量(这里有一些空间，因为我不知道精度或您的规范化过程)，但是对于大多数“大于一”的向量，它都会溢出，毫无意义。
如果将所有数组保存在内存中太大，则可以轻松地将结果数组限制为任意多个字节。很少情况下会在12字节后产生发散结果。
我没有对所有的问题进行过压力测试，比如denormals，和角案例，在代码中有一些对关键点的注释。
当然还有：
你可以很容易地改进每件事，我只是愿意分享这个道理
源代码再一次
主要参考资料：
单精度浮点格式(维基百科)
https://stackoverflow.com/questions/22899563
复制相似问题
腾讯云开发者
