c++ - 排序 float

Question

我有一个法向 vector 列表，正在计算标量三乘积并对它们进行排序。我在三种情况下比较了排序:

使用Matlab排序查找最大的绝对三元乘积值

在C++中使用std::sort函数获取std::vector的乘积值。对三联产品使用双打。

在OpenCL C中使用Radix排序-将绝对浮点值转换为无符号整数并将其转换回。我正在将cl_float用于三重产品。

所有这些都给出了不同的值以及不同的索引，这导致了我的算法出现问题。在这种情况下有什么问题，如何使它们保持一致？

Answer 1

手头的问题:

计算3个3维顶点的标量三乘积，即每个 vector 的每个分量，表示为binary32 float。

能够判断该计算的结果是否大于该计算的其他结果。

到目前为止，一切都很好，但是如果直接将公式应用于 vector ，操作中可能会丢失一些位，因此我们将无法分辨出两个结果。正如@rcgldr指出的那样，排序不是问题，而是精度。

解决浮点舍入问题的一种方法是增加位数，即使用 double。您说您没有 double，所以我们自己做:只要需要，就可以在 unsigned char数组中执行整个计算。

好的，好的，实际考虑因素:

输入由归一化 vector 组成，因此长度不大于1，这意味着没有任何组件大于1

binary32浮点型的指数范围为-127(零，非正规数)至127(或无穷大为128)，但是所有分量的指数都将从-127变为0(否则它们将大于1)。

输入的最大精度为24位。

标量三元积涉及 vector 积和标量积。在 vector 乘积(将首先出现)中，相乘的结果相减，而在标量乘积中，相乘的结果相加。

注意事项2和3告诉我们，整个输入组件系列都可以采用定点格式，其中偏移量为127位，尾数为24位，即19字节。设为24。

为了能够完全表示所有可能的和与减，一个额外的位就足够了(随着进位的到来)，但是要完全表示所有可能的乘法结果，我们需要将位数增加一倍，因此坚决地说加倍大小就足够了表示 vector 乘法，将其乘以三将足以满足标量积的下一个乘法。

这是一个类的草稿，它将浮点数加载为非常大的定点格式，并将符号保留为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的内容，请注意，该数组已加载到其下部，乘法将〜roll〜该数字对应于高字节，我们可以然后恢复我们的 float 。

为了使它有用，我们需要实现求和与乘法的等效，非常简单，只要我们记得我们只能对已被乘以相同次数(如三乘积)或它们的大小的数组求和不匹配。

这样的类会有所作为的一个示例:

考虑以下三个 vector :

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;
}

如果我们针对 vector 的以上两个版本调用该函数，则数组中的结果将不同:

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的精度之后有所不同，这意味着，如果将该数组强制转换为float，则将是相同的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上可以正常工作):

Source Code on IDEONE C++11 code

如果您尚未迁移到C++ 11，并且您自己定义了精确宽度类型 uint8_t， uint16_t， uint32_t和 int32_t，则代码将编译并在C++ 98中运行。

如何使用它？

使用输入简单地调用函数 tripleProduct并使用提供的过载比较器运算符比较结果，还可以使用提供的重载强制转换运算符将 imSoHuge类强制转换为浮点数(在计算三乘积之后)。

您可以为任何排序算法提供该类的数组和比较器。

结论和注意事项:

请注意，浮点乘法现在是两个70+字节长的数组的乘法，这意味着要多数百倍的时钟周期，再加上求和，比较等，这将慢数千倍，但是，这是正确的。

该算法的整个设计是使用归一化 vector (这里有一定空间，因为我不知道精度或归一化过程)，但是对于大多数“大于一”的方法，它都会溢出并且毫无意义。 vector 。

如果将所有结果数组存储在内存中过多，则可以轻松地将结果数组的大小限制为所需的字节数。极少数情况下，产生的结果在〜12bytes后会发散

我还没有对所有事情进行压力测试，例如异常，异常情况，代码中有一些对关键点的注释。

而且当然:

您可以轻松地改善一切，我只是愿意分享推理=)

Source code again

主要引用:

Single-precision floating-point format (Wikipedia)