title: Arithmetic in R faster on numerics as opposed to integers. What's going on?
body: <p>I was in the middle of converting some code that utilized mostly numeric data (i.e. doubles) to integers and did a quick benchmark to see how much efficiency I gained.</p>
  
  <p>To my surprise it was slower... by about 20%. I thought I had done something wrong, but the original code was only a few basic arithmetical operations on moderately sized vectors, so I knew it wasn't that. Maybe my environment was messed up? I restarted fresh, and the same result... integers were less efficient.</p>
  
  <p>This started a series of test and a dive into the rabbit hole. Here is my first test. We sum one million elements using base R's <code>sum</code>. Note that with R version <code>3.5.0</code> the timings are quite a bit different and with v 3.5.1, the timings are about the same (still not what one would expect):</p>
  
  <pre><code>set.seed(123)
  int1e6 &lt;- sample(1:10, 1e6, TRUE)
  dbl1e6 &lt;- runif(1e6, 1, 10)
  
  head(int1e6)
  # [1] 5 3 6 8 6 2
  class(int1e6)
  # [1] "integer"
  
  head(dbl1e6)
  # [1] 5.060628 2.291397 2.992889 5.299649 5.217105 9.769613
  class(dbl1e6)
  #[1] "numeric"
  
  mean(dbl1e6)
  # [1] 5.502034
  mean(int1e6)
  # [1] 5.505185
  
  ## R 3.5.0
  library(microbenchmark)
  microbenchmark(intSum = sum(int1e6), dblSum = sum(dbl1e6), times = 1000)
  Unit: microseconds
    expr      min       lq      mean   median       uq      max neval
  intSum 1033.677 1043.991 1147.9711 1111.438 1200.725 2723.834  1000
  dblSum  817.719  835.486  945.6553  890.529  998.946 2736.024  1000
  
  ## R 3.5.1
  Unit: microseconds
    expr     min       lq      mean   median        uq      max neval
  intSum 836.243 877.7655  966.4443 950.1525  997.9025 2077.257  1000
  dblSum 866.939 904.7945 1015.3445 986.4770 1046.4120 2541.828  1000
  
  class(sum(int1e6))
  # [1] "integer"
  class(sum(dbl1e6))
  #[1] "numeric"
  </code></pre>
  
  <p>From here on out both version 3.5.0 and 3.5.1 give nearly identical results.</p>
  
  <p>Here is our first dive into the <a href="https://github.com/wch/r-source/blob/trunk/src/library/methods/R/BasicFunsList.R#L53" rel="noreferrer">rabbit hole</a>. Along with the documentation for <code>sum</code> (see <code>?sum</code>), we see that <code>sum</code> is simply a generic function that is dispatched via <code>standardGeneric</code>. Digging deeper, we see it eventually calls <code>R_execMethod</code> <a href="https://github.com/wch/r-source/blob/trunk/src/library/methods/src/methods_list_dispatch.c" rel="noreferrer">here</a> on line 516. This is where I get lost. It looks to me, like <code>R_execClosure</code> is called next followed by many different possible branches. I think the standard path is to call <code>eval</code> next, but I'm not sure. My guess is that eventually, a function is called in <a href="https://github.com/wch/r-source/blob/trunk/src/main/arithmetic.c" rel="noreferrer">arithimetic.c</a> but I can't find anything that specifically sums a vector of numbers. Either way, based off of my limited knowledge of method dispatching and <code>C</code> in general, my naive assumption is that a function that looks like the following is called:</p>
  
  <pre><code>template &lt;typename T&gt;
  T sum(vector&lt;T&gt; x) {
  
      T mySum = 0;
  
      for (std::size_t i = 0; i &lt; x.size(); ++i)
          mySum += x[i];
  
      return mySum;
  }
  </code></pre>
  
  <p>I know there is no function overloading or vectors in <code>C</code>, but you get my point. My belief is that eventually, a bunch of the same type of elements are added to an element of the same type and eventually returned. In <code>Rcpp</code> we would have something like:</p>
  
  <pre><code>template &lt;typename typeReturn, typename typeRcpp&gt;
  typeReturn sumRcpp(typeRcpp x) {
      typeReturn mySum = 0;
      unsigned long int mySize = x.size();
  
      for (std::size_t i = 0; i &lt; mySize; ++i)
          mySum += x[i];
  
      return mySum;
  }
  
  // [[Rcpp::export]]
  SEXP mySumTest(SEXP Rx) {
      switch(TYPEOF(Rx)) {
          case INTSXP: {
              IntegerVector xInt = as&lt;IntegerVector&gt;(Rx);
              int resInt = sumRcpp&lt;int&gt;(xInt);
              return wrap(resInt);
          }
          case REALSXP: {
              NumericVector xNum = as&lt;NumericVector&gt;(Rx);
              double resDbl = sumRcpp&lt;double&gt;(xNum);
              return wrap(resDbl);
          }
          default: {
              Rcpp::stop("Only integers and numerics are supported");   
          }
      }
  }
  </code></pre>
  
  <p>And the benchmarks confirm my normal thinking about the inherit efficiency dominance of integers:</p>
  
  <pre><code>microbenchmark(mySumTest(int1e6), mySumTest(dbl1e6))
  Unit: microseconds
               expr      min       lq      mean    median        uq      max neval
  mySumTest(int1e6)  103.455  160.776  185.2529  180.2505  200.3245  326.950   100
  mySumTest(dbl1e6) 1160.501 1166.032 1278.1622 1233.1575 1347.1660 1644.494   100
  </code></pre>
  
  <h1>Binary Operators</h1>
  
  <p>This got me thinking further. Maybe it is just the complexity wrapped around <code>standardGeneric</code> that makes the different data types behave <em>strangely</em>. So, let's skip all that jazz and go straight to the binary operators (<code>+, -, *, /, %/%</code>)</p>
  
  <pre><code>set.seed(321)
  int1e6Two &lt;- sample(1:10, 1e6, TRUE)
  dbl1e6Two &lt;- runif(1e6, 1, 10)
  
  ## addition
  microbenchmark(intPlus = int1e6 + int1e6Two, 
                 dblPlus = dbl1e6 + dbl1e6Two, times = 1000)
  Unit: milliseconds
     expr      min       lq     mean   median       uq      max neval
  intPlus 2.531220 3.214673 3.970903 3.401631 3.668878 82.11871  1000
  dblPlus 1.299004 2.045720 3.074367 2.139489 2.275697 69.89538  1000
  
  ## subtraction
  microbenchmark(intSub = int1e6 - int1e6Two,
                 dblSub = dbl1e6 - dbl1e6Two, times = 1000)
  Unit: milliseconds
    expr      min       lq     mean   median       uq      max neval
  intSub 2.280881 2.985491 3.748759 3.166262 3.379755 79.03561  1000
  dblSub 1.302704 2.107817 3.252457 2.208293 2.382188 70.24451  1000
  
  ## multiplication
  microbenchmark(intMult = int1e6 * int1e6Two, 
                 dblMult = dbl1e6 * dbl1e6Two, times = 1000)
  Unit: milliseconds
     expr      min       lq     mean   median       uq      max neval
  intMult 2.913680 3.573557 4.380174 3.772987 4.077219 74.95485  1000
  dblMult 1.303688 2.020221 3.078500 2.119648 2.299145 10.86589  1000
  
  ## division
  microbenchmark(intDiv = int1e6 %/% int1e6Two,
                 dblDiv = dbl1e6 / dbl1e6Two, times = 1000)
  Unit: milliseconds
    expr      min       lq     mean   median       uq      max neval
  intDiv 2.892297 3.210666 3.720360 3.228242 3.373456 62.12020  1000
  dblDiv 1.228171 1.809902 2.558428 1.842272 1.990067 64.82425  1000
  </code></pre>
  
  <p>The classes are preserved as well:</p>
  
  <pre><code>unique(c(class(int1e6 + int1e6Two), class(int1e6 - int1e6Two),
           class(int1e6 * int1e6Two), class(int1e6 %/% int1e6Two)))
  # [1] "integer"
  
  unique(c(class(dbl1e6 + dbl1e6Two), class(dbl1e6 - dbl1e6Two),
           class(dbl1e6 * dbl1e6Two), class(dbl1e6 / dbl1e6Two)))
  # [1] "numeric"
  </code></pre>
  
  <p>With every case, we see that arithmetic is 40% - 70% faster on numeric data type. What is really strange is that we get an even larger discrepancy when the two vectors being operated on are identical:</p>
  
  <pre><code>microbenchmark(intPlus = int1e6 + int1e6, 
                 dblPlus = dbl1e6 + dbl1e6, times = 1000)
  Unit: microseconds
     expr      min       lq     mean   median       uq      max neval
  intPlus 2522.774 3148.464 3894.723 3304.189 3531.310 73354.97  1000
  dblPlus  977.892 1703.865 2710.602 1767.801 1886.648 77738.47  1000
  
  microbenchmark(intSub = int1e6 - int1e6,
                 dblSub = dbl1e6 - dbl1e6, times = 1000)
  Unit: microseconds
    expr      min       lq     mean   median       uq      max neval
  intSub 2236.225 2854.068 3467.062 2994.091 3214.953 11202.06  1000
  dblSub  893.819 1658.032 2789.087 1730.981 1873.899 74034.62  1000
  
  microbenchmark(intMult = int1e6 * int1e6, 
                 dblMult = dbl1e6 * dbl1e6, times = 1000)
  Unit: microseconds
     expr      min       lq     mean   median       uq      max neval
  intMult 2852.285 3476.700 4222.726 3658.599 3926.264 78026.18  1000
  dblMult  973.640 1679.887 2638.551 1754.488 1875.058 10866.52  1000
  
  microbenchmark(intDiv = int1e6 %/% int1e6,
                 dblDiv = dbl1e6 / dbl1e6, times = 1000)
  Unit: microseconds
    expr      min       lq     mean   median       uq      max neval
  intDiv 2879.608 3355.015 4052.564 3531.762 3797.715 11781.39  1000
  dblDiv  945.519 1627.203 2706.435 1701.512 1829.869 72215.51  1000
  
  unique(c(class(int1e6 + int1e6), class(int1e6 - int1e6),
           class(int1e6 * int1e6), class(int1e6 %/% int1e6)))
  # [1] "integer"
  
  unique(c(class(dbl1e6 + dbl1e6), class(dbl1e6 - dbl1e6),
           class(dbl1e6 * dbl1e6), class(dbl1e6 / dbl1e6)))
  # [1] "numeric"
  </code></pre>
  
  <p>That is nearly a 100% increase with every operator type!!!</p>
  
  <p>How about a regular for loop in base R?</p>
  
  <pre><code>funInt &lt;- function(v) {
      mySumInt &lt;- 0L
      for (element in v)
          mySumInt &lt;- mySumInt + element
      mySumInt
  }
  
  funDbl &lt;- function(v) {
      mySumDbl &lt;- 0
      for (element in v)
          mySumDbl &lt;- mySumDbl + element
      mySumDbl
  }
  
  microbenchmark(funInt(int1e6), funDbl(dbl1e6))
  Unit: milliseconds
            expr      min       lq     mean   median       uq      max neval
  funInt(int1e6) 25.44143 25.75075 26.81548 26.09486 27.60330 32.29436   100
  funDbl(dbl1e6) 24.48309 24.82219 25.68922 25.13742 26.49816 29.36190   100
  
  class(funInt(int1e6))
  # [1] "integer"
  class(funDbl(dbl1e6))
  # [1] "numeric"
  </code></pre>
  
  <p>The difference isn't amazing, but still one would expect the integer sum to outperform the double sum. I really don't know what to think about this.</p>
  
  <p>So my question is:</p>
  
  <blockquote>
    <p>Why exactly do numeric data types outperform integer data types on basic arithmetical operations in base R?</p>
  </blockquote>
  
  <p>Edit. Forgot to mention this:</p>
  
  <pre><code>sessionInfo()
  R version 3.5.1 (2018-07-02)
  Platform: x86_64-apple-darwin15.6.0 (64-bit)
  Running under: macOS High Sierra 10.13.6
  </code></pre>