关于优化浮点计算的建议：求解二次方程

Question

我正在研究一个细粒度的动力学问题。计算量大的部分是下面的函数，它求解一个二次方程来检测两个粒子的碰撞。

我想知道这是否可以轻松优化，或者我是否正在做一些明显愚蠢的事情？例如，使用这些const double x1 = p1->x; 结构来给编译器一个提示是个好主意吗？查看汇编代码，编译器使用 SSE 指令，但我不知道它们是否在任何方面都是最优的（可能不是）。根据分析器，大部分时间都花在计算表达式a、b 和c 上。当您尝试优化如下所示的某些内核功能时，您（通常）会做什么？

void detect_collision_of_pair(struct particle* p1, struct particle* p2){
        const double x1  = p1->x;
        const double y1  = p1->y;
        const double z1  = p1->z;
        const double x2  = p2->x;
        const double y2  = p2->y;
        const double z2  = p2->z;
        const double vx1 = p1->vx;
        const double vy1 = p1->vy;
        const double vz1 = p1->vz;
        const double vx2 = p2->vx;
        const double vy2 = p2->vy;
        const double vz2 = p2->vz;

        const double a = vx1*vx1 - 2.*vx1*vx2 + vx2*vx2 + vy1*vy1 - 2.*vy1*vy2 + vy2*vy2 + vz1*vz1 - 2.*vz1*vz2 + vz2*vz2;
        const double b = 2.*vx1*x1 - 2.*vx2*x1 - 2.*vx1*x2 + 2.*vx2*x2 + 2.*vy1*y1 - 2.*vy2*y1 - 2.*vy1*y2 + 2.*vy2*y2 + 2.*vz1*z1 - 2.*vz2*z1 - 2.*vz1*z2 + 2.*vz2*z2;
        const double c = -4.*particle_radius*particle_radius + x1*x1 - 2.*x1*x2 + x2*x2 + y1*y1 - 2.*y1*y2 + y2*y2 + z1*z1 - 2.*z1*z2 + z2*z2;

        double root = b*b-4.*a*c;
        if (root>=0.){
                root = sqrt(root);
                double time1 = (-b-root)/(2.*a);
                double time2 = (-b+root)/(2.*a);

                if ( (time1>-dt && time1<0.) || (time1<-dt && time2>0) ){
                        double times = -dt;
                        if (time1>-dt || time2<0){
                                if (time1>-dt){
                                        times = time1;
                                }else{
                                        times = time2;
                                }
                        }
                        resolve_collision(p1,p2,times);
                }
        }
}

Answer 1

为什么要扩展所有方程？这不是一个好主意。为什么不计算：

vx = v1x-v2x;
vy = v1y-v2y;
vz = v1z-v2z;
a = vx*vx + vy*vy + vz*vz;

这比例如计算 a 所执行的操作要少得多。 b 和 c 可以做同样的事情。您也可以对位置进行类似的操作，即计算 px、py、pz 作为位置的差异，然后将它们平方。

另一方面，摆脱所有 const 的东西，编译器不需要那些用于局部变量的东西。实际上，您可能不需要复制到局部变量，只需为您需要的东西创建局部变量，例如上面的 vx,vy,vz。您在这里做的太多了，这就是代码花费太多时间的原因:-)

Answer 2

对于计算 a、b 和 c，我看不出有什么明显的错误。使用临时变量 x1 = p1->x 等对一个体面的优化编译器既无帮助也无害，因为在这种情况下，由于未写入 p1 和 p2，因此不存在别名问题（常见的性能问题），并且编译器将决定在计算 a、b 和 c 时将哪些变量保存在寄存器中。如果你想确保 wrt 别名，你总是可以将参数标记为

struct particle * const restrict p1

p2 也是如此。

我想一个（小的？）改进可能是在求和之后将乘法移动 2，从而节省一些乘法。例如

const double b = 2 * (vx1*x1 - vx2*x1 - vx1*x2 + vx2*x2 + vy1*y1 - vy2*y1 - vy1*y2 + vy2*y2 + vz1*z1 - vz2*z1 - vz1*z2 + vz2*z2);

OTOH，您的代码的一个问题是，使用您所做的明显公式计算二次方程的根很容易发生灾难性的取消。参见例如

http://download.oracle.com/docs/cd/E19957-01/806-3568/ncg_goldberg.html

https://secure.wikimedia.org/wikipedia/en/wiki/Quadratic_equation#Floating_point_implementation

Answer 3

除了其他人提到的......你绝对需要使用双打吗？使用浮点数肯定是加快速度的一种方法。或者如果你需要这么高的精度，你可以使用 compile 来编译 64bit 平台；据我所知，32 位编译器使用一些“技巧”来进行 64 位数学运算，因此仅针对 64 位处理器进行编译应该会对您有所帮助，因为涉及的指令应该更少，或者类似的东西。