Advice on optimization floating point calculations: solving a quadratic equation

I'm working on a granular dynamics problem. The computationally expensive part is the function below that solve a quadratic equation to detect the collision of two particles.

I was wondering if this can be easily optimized or if I'm doing something noticeably stupid? For example, is it a good idea to use these const double x1 = p1->x; constructions to give the compiler a hint? Looking at the assembler code, the compiler uses SSE instructions, but I have no idea if they are optimal in any way (probably not). According to the profiler, most of the time is spend in calculating the expressions a, b and c. What do you (in general) do when you are trying to optimize some kernel function like the one below?

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){
  开发者_StackOverflow                              if (time1>-dt){
                                        times = time1;
                                }else{
                                        times = time2;
                                }
                        }
                        resolve_collision(p1,p2,times);
                }
        }
}

Why have you expanded all the equations? This is not a good idea. Why not compute:

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

This is many fewer operations than what you're doing to compute a for example. The same type of thing can be done for b and c. You can also do similar with the positions, i.e. compute px, py,pz as the differences in positions and then square them.

On another note, get rid of all that const stuff, the compiler doesn't need those for local variables. In fact, you probably don't need to copy to local variables, just create locals for the things you need like the vx,vy,vz above. You're doing way too much here, and that's why the code is taking too much time :-)

Wrt to calculating a, b, and c, I don't see anything obviously wrong. Using the temporary variables x1 = p1->x and so forth should neither help nor hurt with a decent optimizing compiler, as in this case there's no aliasing issue (a common performance problem) due to p1 and p2 not being written to, and the compiler will decide which variables to keep in registers when doing the calculations of a,b, and c. If you want to be sure wrt aliasing, you can always mark the arguments as

struct particle * const restrict p1

and the same for p2.

I suppose one (minor?) improvement might be to move the multiplications by 2 after the summing, saving some multiplications. E.g.

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, one problem with your code is that the calculation of the roots of the quadratic equation with the obvious formula as you have done is prone to catastrophic cancellation. See e.g.

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

Besides what the others have mentioned... do you absolutely need to use doubles? Using floats would sure be one way to make this a lot faster. Or if you need that high of a precision, you could use compile it for a 64bit platform; as far as I'm aware, the 32bit compiler uses some "tricks" to do 64bit math, so just compiling for a 64bit processor should help you out there as there should be fewer instructions involved, or something like that.