title: Advice on optimization floating point calculations: solving a quadratic equation
body: <p>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. </p>
  
  <p>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 <code>const double x1  = p1-&gt;x;</code> 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 <code>a</code>, <code>b</code> and <code>c</code>. What do you (in general) do when you are trying to optimize some kernel function like the one below? </p>
  
  <pre><code>void detect_collision_of_pair(struct particle* p1, struct particle* p2){
          const double x1  = p1-&gt;x;
          const double y1  = p1-&gt;y;
          const double z1  = p1-&gt;z;
          const double x2  = p2-&gt;x;
          const double y2  = p2-&gt;y;
          const double z2  = p2-&gt;z;
          const double vx1 = p1-&gt;vx;
          const double vy1 = p1-&gt;vy;
          const double vz1 = p1-&gt;vz;
          const double vx2 = p2-&gt;vx;
          const double vy2 = p2-&gt;vy;
          const double vz2 = p2-&gt;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&gt;=0.){
                  root = sqrt(root);
                  double time1 = (-b-root)/(2.*a);
                  double time2 = (-b+root)/(2.*a);
  
                  if ( (time1&gt;-dt &amp;&amp; time1&lt;0.) || (time1&lt;-dt &amp;&amp; time2&gt;0) ){
                          double times = -dt;
                          if (time1&gt;-dt || time2&lt;0){
                                  if (time1&gt;-dt){
                                          times = time1;
                                  }else{
                                          times = time2;
                                  }
                          }
                          resolve_collision(p1,p2,times);
                  }
          }
  }
  </code></pre>