 
   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);
+  const float mrinv  =  mass*rinv;
+  const float mrinv3 = rinv2*mrinv;
+  const float mrinv5 = rinv2*mrinv3; 
+  const float mrinv7 = rinv2*mrinv5;   // 16
+
+#if 0
+  float  D0  =  mrinv;
+  float  D1  = -mrinv3;
+  float  D2  =  mrinv5*( 3.0f);
+  float  D3  = -mrinv7*(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);

  return acc;
}

            for (int k = 0; k < NI; k++)
              acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i], eps2);
#else
#if 1    /*  a bit faster */
          const float4  Q0  = tex1Dfetch(texMultipole, address + 1);
          const float4  Q1  = tex1Dfetch(texMultipole, address + 2);
          for (int i = 0; i < WARP_SIZE; i++)
          {
            const float4 jQ0 = make_float4(__shfl(Q0.x, i), __shfl(Q0.y, i), __shfl(Q0.z, i), 0.0f);
            const float4 jQ1 = make_float4(__shfl(Q1.x, i), __shfl(Q1.y, i), __shfl(Q1.z, i), 0.0f);
            for (int k = 0; k < NI; k++)
              acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i], jQ0, jQ1, eps2);
          }
#else
          for (int i = 0; i < WARP_SIZE; i++)
          {
            const int address = approxM[n_approx + i] * 3;
