Skip to content

Instantly share code, notes, and snippets.

@jangorecki jangorecki/nan.c
Last active Sep 29, 2019

Embed
What would you like to do?
vectorized is.na and is.nan
#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;
}
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.