             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;

            for (int k = 0; k < NI; k++)
              acc_i[k] = add_acc(acc_i[k], pos_i[k], sh_mass[i], sh_pos[i], Q0, Q1, eps2);
          }
#endif /* _QUADRUPOLE_ */
          apprCount += WARP_SIZE*NI;
        }


  const float rinv  = rsqrtf(r2);
  const float rinv2 = rinv *rinv;
  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;
