/* [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include <math.h>
# include <ieee754.h>
# include "genlib.h"

void mexFunction(int nlhs,       mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{

# define MEXNAME nansmex

 /* out = nansmex(in); */

# define OUT pout[0]

# define IN  prhs[0] /* double, optional */

 CHECK_NARGIN_LE(1);

 if (nrhs>0) CHECK_DOUBLE(IN);

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];

 int k;
 double *pr;
 union ieee754_double u;
 bool *pl;
 if (nrhs == 0) {
  OUT = mxCreateDoubleMatrix(1, 1, mxREAL);
  mxGetPr(OUT)[0] = __builtin_nans("") /* DBL_SNAN */ ;
 } else {
  OUT = mxCreateLogicalArray(mxGetNumberOfDimensions(prhs[0]),mxGetDimensions(prhs[0]));
  pl = mxGetLogicals(OUT);
  pr = mxGetPr(IN);
  for (k = mxGetNumberOfElements(IN); k--; ) {
   /* pl[k] = fpclassify(pr[k]) == FP_SNAN; */
   u.d = pr[k];
   pl[k] = (double)( fpclassify(pr[k]) == FP_NAN && u.ieee_nan.quiet_nan == 0);
  }
 }

 ASSIGN_PLHS;
}

/*
format hex
x = [NaN nans];
disp(x)
format short
disp(x)
disp(isnan(x))
disp(isnans(x))
*/
