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>

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;
}

#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];

  for (size_t i = 0; i < len; ++i) {
    if (rand() % 100 > 50) {
      x[i] = NA_REAL;
    } else {
      x[i] = 1;
    }
  }

  size_t num_na = 0;
  size_t num_na2 = 0;

  timing("R_IsNA2",
    for (size_t i = 0; i < len; ++i) {
      if (R_IsNA2(x[i])) {
        num_na2++;
      }
    }
  )

  timing("R_IsNA",
    for (size_t i = 0; i < len; ++i) {
      if (R_IsNA(x[i])) {
        num_na++;
      }
    }
  )
  setlocale(LC_NUMERIC, "");
  printf("num_na: %'zu\n", num_na);
  printf("num_na2: %'zu\n", num_na2);

  delete[] x;
  return 0;
}
/*
./timing 1000000000
R_IsNA2 secs: 0.792045
R_IsNA secs: 1.060737
num_na: 490,002,335
num_na2: 490,002,335
*/