Check all basic operators results in all possible cases

src/basic_op_div.cpp

@@ -26,6 +26,8 @@
 #include "libdivide.h" // for fast divison by integer constant
 // needed with gcc-3.3.2
 #include <cassert>
+static const std::complex<float> complex_float_nan(sqrt(-1), sqrt(-1));
+static const std::complex<double> complex_double_nan(sqrt(-1), sqrt(-1));
 
 // Div
 // division: left=left/right
@@ -227,7 +229,7 @@ Data_<SpDObj>* Data_<SpDObj>::DivInv(BaseGDL* r) {
   return this;
 }
 
-template<class Sp>  //no need to differentiate Sp Types, as the FP exception is produced only by s
+template<class Sp>
 Data_<Sp>* Data_<Sp>::DivS(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
   Data_* right = static_cast<Data_*> (r);
 
@@ -362,6 +364,24 @@ Data_<SpDULong64>* Data_<SpDULong64>::DivS(BaseGDL* r) { TRACE_ROUTINE(__FUNCTIO
  return this;
 }
 
+//floats & Complex
+template<>
+Data_<SpDFloat>* Data_<SpDFloat>::DivS(BaseGDL* r) {
+#include "snippets/basic_op_DivS.incpp"
+}
+
+template<>
+Data_<SpDDouble>* Data_<SpDDouble>::DivS(BaseGDL* r) {
+#include "snippets/basic_op_DivS.incpp"
+}
+template<>
+Data_<SpDComplex>* Data_<SpDComplex>::DivS(BaseGDL* r) {
+#include "snippets/basic_op_DivSCplx.incpp"
+}
+template<>
+Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::DivS(BaseGDL* r) {
+#include "snippets/basic_op_DivSCplxDbl.incpp"
+}
 // inverse division: left=right/left
 
 template<class Sp>
@@ -389,60 +409,21 @@ Data_<Sp>* Data_<Sp>::DivInvS(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,
 }
 
 template<>
-Data_<SpDFloat>* Data_<SpDFloat>::DivInvS(BaseGDL* r) { TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
-  Data_* right = static_cast<Data_*> (r);
-
-  ULong nEl = N_Elements();
-  assert(nEl); 
-  GDLStartRegisteringFPExceptions();
-
-  Ty s = (*right)[0];
-  SizeT i = 0;
-  if (nEl == 1) {
-	(*this)[0] = s / (*this)[0];
-	GDLStopRegisteringFPExceptions();
-	return this;
-  }
-  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
-	for (OMPInt ix = i; ix < nEl; ++ix) (*this)[ix] = s / (*this)[ix];
-  } else {
-	TRACEOMP(__FILE__, __LINE__)
-#pragma omp parallel for num_threads(GDL_NTHREADS)
-	  for (OMPInt ix = i; ix < nEl; ++ix) (*this)[ix] = s / (*this)[ix];
-  }
-
-  GDLStopRegisteringFPExceptions();
-
-  return this;
+Data_<SpDFloat>* Data_<SpDFloat>::DivInvS(BaseGDL* r) {
+#include "snippets/basic_op_DivInvSCplx.incpp"
 }
 
 template<>
 Data_<SpDDouble>* Data_<SpDDouble>::DivInvS(BaseGDL* r) {
-  TRACE_ROUTINE(__FUNCTION__, __FILE__, __LINE__)
-  Data_* right = static_cast<Data_*> (r);
-
-  ULong nEl = N_Elements();
-  assert(nEl);
-  GDLStartRegisteringFPExceptions();
-
-  Ty s = (*right)[0];
-  SizeT i = 0;
-  if (nEl == 1) {
-	(*this)[0] = s / (*this)[0];
-	GDLStopRegisteringFPExceptions();
-	return this;
-  }
-  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
-	for (OMPInt ix = i; ix < nEl; ++ix) (*this)[ix] = s / (*this)[ix];
-  } else {
-	TRACEOMP(__FILE__, __LINE__)
-#pragma omp parallel for num_threads(GDL_NTHREADS)
-	  for (OMPInt ix = i; ix < nEl; ++ix) (*this)[ix] = s / (*this)[ix];
-  }
-
-  GDLStopRegisteringFPExceptions();
-
-  return this;
+#include "snippets/basic_op_DivInvSCplx.incpp"
+}
+template<>
+Data_<SpDComplex>* Data_<SpDComplex>::DivInvS(BaseGDL* r) {
+#include "snippets/basic_op_DivInvSCplx.incpp"
+}
+template<>
+Data_<SpDComplexDbl>* Data_<SpDComplexDbl>::DivInvS(BaseGDL* r) {
+#include "snippets/basic_op_DivInvSCplx.incpp"
 }
 
 // invalid types

src/snippets/basic_op_AndOpCplx.incpp

@@ -0,0 +1,17 @@
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl);
+  if (nEl == 1) {
+    if ((*right)[0] == zero) (*this)[0] = zero;
+    return this;
+  }
+
+  if ((GDL_NTHREADS=parallelize( nEl))==1) {
+    for (OMPInt i = 0; i < nEl; ++i) if ((*right)[i] == zero) (*this)[i] = zero;
+  } else {
+    TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+      for (OMPInt i = 0; i < nEl; ++i) if ((*right)[i] == zero) (*this)[i] = zero;
+  }
+  return this;

src/snippets/basic_op_AndOpInv.incpp

@@ -0,0 +1,18 @@
+ TRACE_ROUTINE(__FUNCTION__, __FILE__, __LINE__)
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl);
+  if (nEl == 1) {
+	if ((*this)[0] != zero) (*this)[0] = (*right)[0];
+	return this;
+  }
+
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] != zero) (*this)[i] = (*right)[i];
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] != zero) (*this)[i] = (*right)[i];
+  }
+  return this;

src/snippets/basic_op_AndOpInvS.incpp

@@ -0,0 +1,27 @@
+TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+
+  assert(nEl);
+  Ty s = (*right)[0];
+  if (s == zero) {
+    {
+      for (SizeT i = 0; i < nEl; ++i)
+        (*this)[i] = zero;
+    }
+  } else {
+    if (nEl == 1) {
+      if ((*this)[0] != zero) (*this)[0] = s;
+      return this;
+    }
+
+    if ((GDL_NTHREADS=parallelize( nEl))==1) {
+      for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] != zero) (*this)[i] = s;
+    } else {
+      TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+        for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] != zero) (*this)[i] = s;
+    }
+  }
+  return this;

src/snippets/basic_op_DivInvSCplx.incpp

@@ -0,0 +1,24 @@
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  GDLStartRegisteringFPExceptions();
+
+  Ty s = (*right)[0];
+  SizeT i = 0;
+  if (nEl == 1) {
+	(*this)[0] = s / (*this)[0];
+	GDLStopRegisteringFPExceptions();
+	return this;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (OMPInt ix = i; ix < nEl; ++ix) (*this)[ix] = s / (*this)[ix];
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = i; ix < nEl; ++ix) (*this)[ix] = s / (*this)[ix];
+  }
+
+  GDLStopRegisteringFPExceptions();
+
+  return this;

src/snippets/basic_op_DivInvSNew.incpp

@@ -0,0 +1,26 @@
+ TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+ Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  Data_* res = NewResult();
+  GDLStartRegisteringFPExceptions();
+
+  Ty s = (*right)[0];
+  SizeT i = 0;
+  if (nEl == 1) {
+	(*res)[0] = s / (*this)[0];
+	GDLStopRegisteringFPExceptions();
+	return res;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (OMPInt ix = i; ix < nEl; ++ix) (*res)[ix] = s / (*this)[ix];
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = i; ix < nEl; ++ix) (*res)[ix] = s / (*this)[ix];
+  }
+
+  GDLStopRegisteringFPExceptions();
+
+  return res;

src/snippets/basic_op_DivS.incpp

@@ -0,0 +1,23 @@
+TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  Ty s = (*right)[0];
+
+  GDLStartRegisteringFPExceptions();
+  if (nEl == 1) {
+	(*this)[0]/= s;
+	GDLStopRegisteringFPExceptions();
+	return this;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (SizeT ix = 0; ix < nEl; ++ix) (*this)[ix] /= s;
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = 0; ix < nEl; ++ix) (*this)[ix] /= s;
+  }
+
+  GDLStopRegisteringFPExceptions();
+  return this;

src/snippets/basic_op_DivSCplx.incpp

@@ -0,0 +1,30 @@
+TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  Ty s = (*right)[0];
+  // the general complex division code below will cost time and produce a slightly different result
+  // from IDL when s=complex(0,0). Best to treat that here directly
+  if (s == this->zero) {
+    for (SizeT ix = 0; ix < nEl; ++ix) (*this)[ix]=complex_float_nan;
+	  GDLRegisterADivByZeroException();
+	  return this; 
+  }
+
+  GDLStartRegisteringFPExceptions();
+  if (nEl == 1) {
+	(*this)[0]/= s;
+	GDLStopRegisteringFPExceptions();
+	return this;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (SizeT ix = 0; ix < nEl; ++ix) (*this)[ix] /= s;
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = 0; ix < nEl; ++ix) (*this)[ix] /= s;
+  }
+
+  GDLStopRegisteringFPExceptions();
+  return this;

src/snippets/basic_op_DivSCplxDbl.incpp

@@ -0,0 +1,30 @@
+TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  Ty s = (*right)[0];
+  // the general complex division code below will cost time and produce a slightly different result
+  // from IDL when s=complex(0,0). Best to treat that here directly
+  if (s == this->zero) {
+    for (SizeT ix = 0; ix < nEl; ++ix) (*this)[ix]=complex_double_nan;
+	  GDLRegisterADivByZeroException();
+	  return this; 
+  }
+
+  GDLStartRegisteringFPExceptions();
+  if (nEl == 1) {
+	(*this)[0]/= s;
+	GDLStopRegisteringFPExceptions();
+	return this;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (SizeT ix = 0; ix < nEl; ++ix) (*this)[ix] /= s;
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = 0; ix < nEl; ++ix) (*this)[ix] /= s;
+  }
+
+  GDLStopRegisteringFPExceptions();
+  return this;

src/snippets/basic_op_DivSNew.incpp

@@ -0,0 +1,24 @@
+TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  Ty s = (*right)[0];
+  Data_* res=NewResult(); 
+
+  GDLStartRegisteringFPExceptions();
+  if (nEl == 1) {
+	(*res)[0]=(*this)[0] / s;
+	GDLStopRegisteringFPExceptions();
+	return res;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (SizeT ix = 0; ix < nEl; ++ix) (*res)[ix]=(*this)[ix] / s;
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = 0; ix < nEl; ++ix) (*res)[ix]=(*this)[ix] / s;
+  }
+
+  GDLStopRegisteringFPExceptions();
+  return res;

src/snippets/basic_op_DivSNewCplx.incpp

@@ -0,0 +1,32 @@
+  TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  Ty s = (*right)[0];
+  // the general complex division code below will cost time and produce a slightly different result
+  // from IDL when s=complex(0,0). Best to treat that here directly
+  Data_* res; 
+  if (s == this->zero) {
+	res = this->Dup(); //faster 
+    for (SizeT ix = 0; ix < nEl; ++ix) (*res)[ix]=complex_float_nan;
+	GDLRegisterADivByZeroException();
+	return res; 
+  }
+  res = NewResult(); 
+  GDLStartRegisteringFPExceptions();
+  if (nEl == 1) {
+	(*res)[0]=(*this)[0] / s;
+	GDLStopRegisteringFPExceptions();
+	return res;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (SizeT ix = 0; ix < nEl; ++ix) (*res)[ix]=(*this)[ix] / s;
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = 0; ix < nEl; ++ix) (*res)[ix]=(*this)[ix] / s;
+  }
+
+  GDLStopRegisteringFPExceptions();
+  return res;

src/snippets/basic_op_DivSNewCplxDbl.incpp

@@ -0,0 +1,32 @@
+  TRACE_ROUTINE(__FUNCTION__,__FILE__,__LINE__)
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  assert(nEl); 
+  Ty s = (*right)[0];
+  // the general complex division code below will cost time and produce a slightly different result
+  // from IDL when s=complex(0,0). Best to treat that here directly
+  Data_* res; 
+  if (s == this->zero) {
+	res = this->Dup(); //faster 
+    for (SizeT ix = 0; ix < nEl; ++ix) (*res)[ix]=complex_double_nan;
+	GDLRegisterADivByZeroException();
+	return res; 
+  }
+  res = NewResult(); 
+  GDLStartRegisteringFPExceptions();
+  if (nEl == 1) {
+	(*res)[0]=(*this)[0] / s;
+	GDLStopRegisteringFPExceptions();
+	return res;
+  }
+  if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
+	for (SizeT ix = 0; ix < nEl; ++ix) (*res)[ix]=(*this)[ix] / s;
+  } else {
+	TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+	  for (OMPInt ix = 0; ix < nEl; ++ix) (*res)[ix]=(*this)[ix] / s;
+  }
+
+  GDLStopRegisteringFPExceptions();
+  return res;

src/snippets/basic_op_GtMarkNewCplx.incpp

@@ -0,0 +1,22 @@
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  Data_* res = NewResult();
+
+  assert(nEl);
+
+  if (nEl == 1) {
+    if (abs((*this)[0]) < abs((*right)[0])) (*res)[0] = (*right)[0];
+    else (*res)[0] = (*this)[0];
+    return res;
+  }
+  if ((GDL_NTHREADS=parallelize( nEl))==1) {
+    for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) < abs((*right)[i])) (*res)[i] = (*right)[i];
+      else (*res)[i] = (*this)[i];
+  } else {
+    TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+      for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) < abs((*right)[i])) (*res)[i] = (*right)[i];
+      else (*res)[i] = (*this)[i];
+  }
+  return res;

src/snippets/basic_op_GtMarkSNewCplx.incpp

@@ -0,0 +1,22 @@
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  Data_* res = NewResult();
+  assert(nEl);
+  if (nEl == 1) {
+    if (abs((*this)[0]) < abs((*right)[0])) (*res)[0] = (*right)[0];
+    else (*res)[0] = (*this)[0];
+    return res;
+  }
+  Ty s=(*right)[0];
+  DDouble v= abs(s); // note: we compare to doubles.
+   if ((GDL_NTHREADS=parallelize( nEl))==1) {
+    for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) < v) (*res)[i] = s;
+      else (*res)[i] = (*this)[i];
+  } else {
+    TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+      for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) < v) (*res)[i] = s;
+      else (*res)[i] = (*this)[i];
+  };
+  return res;

src/snippets/basic_op_LtMarkNewCplx.incpp

@@ -0,0 +1,22 @@
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  Data_* res = NewResult();
+
+  assert(nEl);
+
+  if (nEl == 1) {
+    if (abs((*this)[0]) > abs((*right)[0])) (*res)[0] = (*right)[0];
+    else (*res)[0] = (*this)[0];
+    return res;
+  }
+  if ((GDL_NTHREADS=parallelize( nEl))==1) {
+    for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) > abs((*right)[i])) (*res)[i] = (*right)[i];
+      else (*res)[i] = (*this)[i];
+  } else {
+    TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+      for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) > abs((*right)[i])) (*res)[i] = (*right)[i];
+      else (*res)[i] = (*this)[i];
+  }
+  return res;

src/snippets/basic_op_LtMarkSNewCplx.incpp

@@ -0,0 +1,22 @@
+  Data_* right = static_cast<Data_*> (r);
+
+  ULong nEl = N_Elements();
+  Data_* res = NewResult();
+  assert(nEl);
+  if (nEl == 1) {
+    if (abs((*this)[0]) > abs((*right)[0])) (*res)[0] = (*right)[0];
+    else (*res)[0] = (*this)[0];
+    return res;
+  }
+  Ty s=(*right)[0];
+  DDouble v= abs(s); // note: we compare to doubles.
+   if ((GDL_NTHREADS=parallelize( nEl))==1) {
+    for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) > v) (*res)[i] = s;
+      else (*res)[i] = (*this)[i];
+  } else {
+    TRACEOMP(__FILE__, __LINE__)
+#pragma omp parallel for num_threads(GDL_NTHREADS)
+      for (OMPInt i = 0; i < nEl; ++i) if (abs((*this)[i]) > v) (*res)[i] = s;
+      else (*res)[i] = (*this)[i];
+  };
+  return res;

src/snippets/basic_op_OrOpSCplx.incpp

@@ -5,16 +5,16 @@ assert(nEl);
 Ty s = (*right)[0];
 if (s != zero) {
   if (nEl == 1) {
-    if ((*this)[0] == zero) (*this)[0] = s;
+    (*this)[0] = s;
     return this;
   }
 
   if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
-    for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] == zero) (*this)[i] = s;
+    for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = s;
   } else {
     TRACEOMP(__FILE__, __LINE__)
 #pragma omp parallel for num_threads(GDL_NTHREADS)
-      for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] == zero) (*this)[i] = s;
+      for (OMPInt i = 0; i < nEl; ++i) (*this)[i] = s;
   }
 }
 return this;

src/snippets/basic_op_OrOpSNewCplx.incpp

@@ -4,23 +4,18 @@ ULong nEl = N_Elements();
 Data_* res = NewResult();
 assert(nEl);
 Ty s = (*right)[0];
-// right->Scalar(s);
 if (s != zero) {
   if (nEl == 1) {
-    if ((*this)[0] == zero) (*res)[0] = s;
-    else (*res)[0] = (*this)[0];
+    (*res)[0] = s;
     return res;
   }
   if ((GDL_NTHREADS = parallelize(nEl)) == 1) {
-    for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] == zero) (*res)[i] = s;
-      else (*res)[i] = (*this)[i];
+    for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = s;
   } else {
     TRACEOMP(__FILE__, __LINE__)
 #pragma omp parallel for num_threads(GDL_NTHREADS)
-      for (OMPInt i = 0; i < nEl; ++i) if ((*this)[i] == zero) (*res)[i] = s;
-      else (*res)[i] = (*this)[i];
+      for (OMPInt i = 0; i < nEl; ++i) (*res)[i] = s;
   }
   return res;
 }
-// s == zero
 return this->Dup();

testsuite/LIST

@@ -1,4 +1,5 @@
 test_all_basic_functions.pro
+test_indepth_basic_functions.pro
 test_zero_divide.pro
 test_angles.pro
 test_antlr_issues.pro