Untitled

Advice on optimization floating point calculations: solving a quadratic equation

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

                }

        }

}

       

vx = v1x-v2x;

vy = v1y-v2y;

vz = v1z-v2z;

a = vx*vx + vy*vy + vz*vz;

       

struct particle * const restrict p1

       

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