vectorized is.na and is.nan

nan.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

typedef union {
  double value;
  unsigned int word[2];
} ieee_double;
static double R_ValueOfNA(void) {
  volatile ieee_double x;
  x.word[1] = 0x7ff00000;
  x.word[0] = 1954;
  return x.value;
}
#  define ISNAN(x)     (isnan(x)!=0)
int R_IsNA(double x) {
  if (isnan(x)) {
    ieee_double y;
    y.value = x;
    return (y.word[0] == 1954);
  }
  return 0;
}
int R_IsNaN(double x) {
  if (isnan(x)) {
    ieee_double y;
    y.value = x;
    return (y.word[0] != 1954);
  }
  return 0;
}
void R_vIsNA(double *x, int64_t nx, int *ans) {
  ieee_double y;
  for (int64_t i=0; i<nx; i++) {
    if (!isnan(x[i])) {
      ans[i] = 0;
    } else {
      y.value = x[i];
      ans[i] = y.word[0]==1954;
    }
  }
}
void R_vIsNaN(double *x, int64_t nx, int *ans) {
  ieee_double y;
  for (int64_t i=0; i<nx; i++) {
    if (!isnan(x[i])) {
      ans[i] = 0;
    } else {
      y.value = x[i];
      ans[i] = y.word[0]!=1954;
    }
  }
}

int main(int argc, char *argv[]) {
  if (argc < 3) {printf("usage: ./a.out ISNAN 6\n\n  runs ISNAN(x) for x of size 1e6\n\nsupported funs are ISNAN, R_IsNA, R_IsNaN and two new ones R_vIsNA, R_vIsNaN\n"); return 1;}

  int fun = 0;
  if (strcmp(argv[1], "ISNAN")==0) fun = 0;
  else if (strcmp(argv[1], "R_IsNA")==0) fun = 1;
  else if (strcmp(argv[1], "R_IsNaN")==0) fun = 2;
  else if (strcmp(argv[1], "R_vIsNA")==0) fun = 3;
  else if (strcmp(argv[1], "R_vIsNaN")==0) fun = 4;
  else {printf("unsupported function to test, see usage by calling ./a.out"); return 1;}
  
  int64_t nexp = atoi(argv[2]);
  int64_t n = pow(10, nexp);
  printf("running %s size 10^%ld: %ld\n", argv[1], nexp, n);
  
  struct timespec start, finish;
  double elapsed;
  
  // alloc memory for ans array storing 0/1
  clock_gettime(CLOCK_MONOTONIC, &start);
  int *ans = malloc(n*sizeof(int));
  clock_gettime(CLOCK_MONOTONIC, &finish);
  elapsed = (finish.tv_sec - start.tv_sec);
  elapsed += (finish.tv_nsec - start.tv_nsec) / 1000000000.0;
  double ans_alloc_elapsed = elapsed;
  if (!ans) {printf("failed to alloc ans: %.3f GB\n", ((double) (n*sizeof(int)))/1073741824); return 1;}
  else printf("ans allocated: %.3f GB (took %.3f)\n", ((double) (n*sizeof(int)))/1073741824, ans_alloc_elapsed);
  
  // genreate data
  double *x = malloc(n*sizeof(double));
  if (!x) {printf("failed to alloc x  : %.3f GB\n", ((double) (n*sizeof(double)))/1073741824); return 1;}
  else printf("x   allocated: %.3f GB\n", ((double) (n*sizeof(double)))/1073741824);
  double NA = R_ValueOfNA();
  for (int64_t i=0; i<n; i++) {
    if (i % 5) {
      x[i] = ((double) i)/2;
    } else if (i % 3) {
      x[i] = NA;
    } else {
      x[i] = NAN;
    }
  }

  if (fun==0) { // ISNAN
    clock_gettime(CLOCK_MONOTONIC, &start);
    for (int64_t i=0; i<n; i++) ans[i] = ISNAN(x[i]);
    clock_gettime(CLOCK_MONOTONIC, &finish);
  } else if (fun==1) { // R_IsNA
    clock_gettime(CLOCK_MONOTONIC, &start);
    for (int64_t i=0; i<n; i++) ans[i] = R_IsNA(x[i]);
    clock_gettime(CLOCK_MONOTONIC, &finish);
  } else if (fun==2) { // R_IsNaN
    clock_gettime(CLOCK_MONOTONIC, &start);
    for (int64_t i=0; i<n; i++) ans[i] = R_IsNaN(x[i]);
    clock_gettime(CLOCK_MONOTONIC, &finish);
  } else if (fun==3) { // R_vIsNA
    clock_gettime(CLOCK_MONOTONIC, &start);
    R_vIsNA(x, n, ans);
    clock_gettime(CLOCK_MONOTONIC, &finish);
  } else if (fun==4) { // R_vIsNaN
    clock_gettime(CLOCK_MONOTONIC, &start);
    R_vIsNaN(x, n, ans);
    clock_gettime(CLOCK_MONOTONIC, &finish);
  }
  elapsed = (finish.tv_sec - start.tv_sec);
  elapsed += (finish.tv_nsec - start.tv_nsec) / 1000000000.0;
  int w = 0;
  for (int64_t i=0; i<n; i++) w += ans[i];
  printf("%s test elapsed: %.3fs (count: %d)\n", argv[1], elapsed, w);
  return 0;
}

run.sh

gcc nan.c -lm

./a.out R_vIsNA 8
./a.out R_IsNA 8
./a.out R_vIsNaN 8
./a.out R_IsNaN 8
./a.out ISNAN 8