c - 关于优化浮点计算的建议 : solving a quadratic equation

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。你在这里做的太多了，这就是代码花费太多时间的原因:-)