title: better algorithm for counting points in grids
body: <p>Suppose I have a point P in [0,1]*[0,1], and [0,1] is divided into m(say 200) grids. I use <code>A[m][m]</code> to indicate whether [a small square centred at P with length 2h] covers each grid or not. So for a point P, <code>A[i][j]</code> is either (increase by) 1 or 0.</p>
  
  <p>Suppose I have n such points(P1,...,Pn), I want to calculate A(for each point Pi, I redo the above procedure, adding 1 or not). How can I do this efficiently(with C++) rather than writing 3 layers of for loops to check for each grid and each point(So O(nm^2))?</p>
  
  <p>I tried the naive 3 for loops with C++. It takes longer time than using some of the vectorized operations(like <code>vector&lt;= number</code> for comparing n numbers together, <code>A[bool vector, bool vector]</code> for subsetting) in R. </p>
  
  <p>Since C++ is generally faster than R, is there any smart way to implement this process?</p>
  
  <pre><code>#include &lt;Rcpp.h&gt;
  #include &lt;cmath&gt;
  using namespace Rcpp;
  
  // [[Rcpp::export]]
  
  double myfun(NumericVector u, NumericVector v)
  {
      double n = u.size();
      double A[200][200] = {0};
      double pos[200];
      int i = 0, j = 0, k = 0;
  
      for (i = 0; i &lt; 200; i++)
      {
          pos[i] = (double)i / 201;
      }
  
      for (k = 0; k &lt; n; k++)
      {
          for (i = 0; i &lt; 200; i++)
          {
              for (j = 0; j &lt; 200; j++)
              {
                  if ( (fabs(u[k] - pos[i]) &lt;= h) &amp;&amp; (fabs(v[k] - pos[j]) &lt;=h ) )
                  {
                      A[i][j]++;
                  }
              }
          }
      }
  
      double s = 0, avg = 0;
  
      for (i = 0; i &lt;200; i++)
      {
          for (j = 0; j &lt; 200; j++)
          {
              s += A[i][j];
          }
      }
      avg = s / (200 * 200);
      return (avg);
  }
  </code></pre>