timing.cc

/* clang++ -std=c++14 timing.cc -o timing -g -Wall -O3 */

#include <time.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <locale.h>
#include <algorithm>

typedef union
{
  double value;
  unsigned int word[2];
} ieee_double;

double R_ValueOfNA(void)
{
  volatile ieee_double x;
  x.word[1] = 0x7ff00000;
  x.word[0] = 1954;
  return x.value;
}

static double NA_REAL = R_ValueOfNA();

int R_IsNA(double x)
{
  if (isnan(x)) {
    ieee_double y;
    y.value = x;
    return (y.word[0] == 1954);
  }
  return 0;
}

uint64_t na = *reinterpret_cast<uint64_t*>(&NA_REAL) ;

// the mask to nuke the 13th bit
uint64_t mask = ~( (uint64_t(1) << 51 ) );

int R_IsNA2(double x)
{
  return ( *reinterpret_cast<uint64_t*>(&x) & mask ) == na;
}

int R_IsNA3(uint64_t x)
{
  return ( x & mask ) == na;
}

#define timing(name, a) start=clock(); a; cpu_time_used = ((double) (clock() - start)) / CLOCKS_PER_SEC; printf("%s secs: %f\n", name, cpu_time_used);

int main(int argc, char* argv[]) {
  
  time_t t;
  /* Initializes random number generator */
  srand((unsigned) time(&t));
  
  clock_t start;
  double cpu_time_used;
  
  size_t len = atol(argv[1]);
  double* x = new double[len];
  
  timing("generation",
    for (size_t i = 0; i < len; ++i) {
      if (rand() % 100 > 50) {
        x[i] = NA_REAL;
      } else {
        x[i] = (double)rand() / RAND_MAX ;
      }
    }
  ) ;
  
  size_t num_na  = 0;
  size_t num_na2 = 0;
  size_t num_na3 = 0;
  
  // first run is not timed
  size_t dummy = std::count_if( x, x+len, R_IsNA ) ;
  
  
  timing("R_IsNA",
    num_na = std::count_if( x, x+len, R_IsNA )
  )
    
  timing("R_IsNA2",
    num_na2 = std::count_if( x, x+len, R_IsNA2 )
  )
    
  timing("R_IsNA3", 
    num_na3 = std::count_if( 
      reinterpret_cast<uint64_t*>(x), 
      reinterpret_cast<uint64_t*>(x+len), 
      R_IsNA3 
    )
  )  
    
  setlocale(LC_NUMERIC, "");
  printf("num_na: %'zu\n", num_na);
  printf("num_na2: %'zu\n", num_na2);
  printf("num_na2: %'zu\n", num_na3);
  printf("num_na2: %'zu\n", dummy); // just so that it is not optimized away
  
  delete[] x;
  return 0;
}