 
   float Qii = oct_q11 + oct_q22 + oct_q33;
   float QijRiRj =
-    (oct_q11*dr.x*dr.x + oct_q22*dr.y*dr.y + oct_q33*dr.z*dr.z) +
+         (oct_q11*dr.x*dr.x + oct_q22*dr.y*dr.y + oct_q33*dr.z*dr.z) +
     2.0f*(oct_q12*dr.y*dr.x + oct_q13*dr.z*dr.x + oct_q23*dr.y*dr.z);
 
   acc.w        -= D0 + 0.5f*D1*Qii + 0.5f*D2*QijRiRj;


  const float rinv  = rsqrtf(r2);
  const float rinv2 = rinv *rinv;
  const float rinv3 = rinv *rinv2;
  const float rinv5 = rinv2*rinv3;
  const float rinv7 = rinv2*rinv5;

  float  D0  =  rinv*mass;
  float  D1  = -rinv3*mass;
  float  D2  =  rinv5*mass*( 3.0f);
  float  D3  = -rinv7*mass*(15.0f);

  float oct_q11 = Q0.x;
  float oct_q22 = Q0.y;

  acc.x         += C01a*dr.x + D2*(oct_q11*dr.x + oct_q12*dr.y + oct_q13*dr.z);
  acc.y         += C01a*dr.y + D2*(oct_q12*dr.x + oct_q22*dr.y + oct_q23*dr.z);
  acc.z         += C01a*dr.z + D2*(oct_q13*dr.x + oct_q23*dr.y + oct_q33*dr.z);
#else
  float  D0  =  mrinv;
  float  D1  = -mrinv3;
  float  D2  =  mrinv5*( 3.0f);
  float  D3  = -mrinv7*(15.0f); // 3

  const float q11 = Q0.x;
  const float q22 = Q0.y;
  const float q33 = Q0.z;
  const float q12 = Q1.x;
  const float q13 = Q1.y;
  const float q23 = Q1.z;

  const float  q  = q11 + q22 + q33;
  const float3 qR = make_float3(
      q11*dr.x + q12*dr.y + q13*dr.z,
      q12*dr.x + q22*dr.y + q23*dr.z,
      q13*dr.x + q23*dr.y + q33*dr.z);
  const float qRR = qR.x*dr.x + qR.y*dr.y + qR.z*dr.z;  // 22

  acc.w  -= D0 + 0.5f*(D1*q + D2*qRR);
  float C = D1 + 0.5f*(D2*q + D3*qRR);
  acc.x  += C*dr.x + D2*qR.x;
  acc.y  += C*dr.y + D2*qR.y;
  acc.z  += C*dr.z + D2*qR.z;               // 23
#endif  /* total: 16 + 3 + 22 + 23 = 64 flops */

  return acc;
}
