Big overhaul and simplification of the NumericTraits system.

Add a utility header for random numbers (might be merged into NumericTraits)

tvmet-1.7.1/include/tvmet/NumericTraits.h

@@ -27,193 +27,79 @@
 #if defined(EIGEN_USE_COMPLEX)
 #  include <complex>
 #endif
+
 #include <cmath>
 #include <limits>
 
 #include <tvmet/CompileTimeError.h>
-
+#include <tvmet/util/Random.h>
 
 namespace tvmet {
 
-
 /**
  * \class NumericTraits NumericTraits.h "tvmet/NumericTraits.h"
  * \brief Traits for integral types for operations.
  *
  * For each type we have to specialize this traits.
- *
- * \note Keep in mind that the long types long long and long double doesn't
- *       have traits. This is due to the sum_type. We can't give a guarantee
- *       that there is a type of holding the sum. Therefore using this traits
- *       is only safe if you have long long resp. long double types by
- *       working on long ints and doubles. Otherwise you will get not expected
- *       result for some circumstances. Anyway, you can use big integer/float
- *       libraries and specialize the traits by your own.
- *
- * \todo The abs function of complex<non_float_type> can have an
- *       overrun due to numeric computation. Solve it (someone
- *       using value_type=long here?)
  */
 template<class T>
 struct NumericTraits {
-  typedef T					base_type;
+  typedef T real_type;
-  typedef T					value_type;
+  typedef T value_type;
-  typedef value_type				sum_type;
+  typedef T float_type;
-  typedef value_type				diff_type;
+  typedef const T & argument_type;
-  typedef value_type				float_type;
-  typedef value_type				signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef const value_type&			argument_type;
 
   static inline
-  base_type real(argument_type x);
+  real_type real(argument_type x);
 
   static inline
-  base_type imag(argument_type x);
+  real_type imag(argument_type x);
 
   static inline
   value_type conj(argument_type x);
 
   static inline
-  base_type abs(argument_type x);
+  real_type abs(argument_type x);
 
   static inline
   value_type sqrt(argument_type x);
 
-  static inline
+  enum{ is_complex = false };
-  base_type norm_1(argument_type x) {
-    return NumericTraits<base_type>::abs(traits_type::real(x))
-         + NumericTraits<base_type>::abs(traits_type::imag(x));
-  }
-
-  static inline
-  base_type norm_2(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_inf(argument_type x) {
-    return std::max(NumericTraits<base_type>::abs(traits_type::real(x)),
-		    NumericTraits<base_type>::abs(traits_type::imag(x)));
-   }
-
-  static inline
-  bool equals(argument_type lhs, argument_type rhs) {
-    static base_type sqrt_epsilon(
-      NumericTraits<base_type>::sqrt(
-        std::numeric_limits<base_type>::epsilon()));
-
-    return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon *
-      std::max(std::max(traits_type::norm_inf(lhs),
-			traits_type::norm_inf(rhs)),
-	       std::numeric_limits<base_type>::min());
-  }
 };
 
-
 /*
- * numeric traits for standard types
+ * numeric traits for built-in types
  */
 
-
-/**
- * \class NumericTraits<char> NumericTraits.h "tvmet/NumericTraits.h"
- * \brief Traits specialized for char.
- */
-template<>
-struct NumericTraits<char> {
-  typedef char					value_type;
-  typedef value_type				base_type;
-  typedef long					sum_type;
-  typedef int					diff_type;
-  typedef float					float_type;
-  typedef char					signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef value_type				argument_type;
-
-  static inline
-  base_type real(argument_type x) { return x; }
-
-  static inline
-  base_type imag(argument_type x) { TVMET_UNUSED(x); return 0; }
-
-  static inline
-  value_type conj(argument_type x) { return x; }
-
-  static inline
-  base_type abs(argument_type x) { return std::abs(x); }
-
-  static inline
-  value_type sqrt(argument_type x) {
-    return static_cast<value_type>(std::sqrt(static_cast<float_type>(x)));
-  }
-
-  static inline
-  base_type norm_1(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_2(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_inf(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; }
-
-  enum { is_complex = false };
-
-  /** Complexity on operations. */
-  enum {
-    ops_plus = 1,	/**< Complexity on plus/minus ops. */
-    ops_muls = 1	/**< Complexity on multiplications. */
-  };
-};
-
-
 /**
  * \class NumericTraits<int> NumericTraits.h "tvmet/NumericTraits.h"
  * \brief Traits specialized for int.
  */
 template<>
 struct NumericTraits<int> {
-  typedef int 					value_type;
+  typedef int value_type;
-  typedef value_type				base_type;
+  typedef value_type real_type;
-  typedef long 					sum_type;
+  typedef double float_type;
-  typedef int			 		diff_type;
+  typedef value_type argument_type;
-  typedef double			 	float_type;
-  typedef int					signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef value_type				argument_type;
 
   static inline
-  base_type real(argument_type x) { return x; }
+  real_type real(argument_type x) { return x; }
 
   static inline
-  base_type imag(argument_type x) { TVMET_UNUSED(x); return 0; }
+  real_type imag(argument_type x) { TVMET_UNUSED(x); return 0; }
 
   static inline
   value_type conj(argument_type x) { return x; }
 
-  static inline
-  base_type abs(argument_type x) { return std::abs(x); }
-
   static inline
   value_type sqrt(argument_type x) {
     return static_cast<value_type>(std::sqrt(static_cast<float_type>(x)));
   }
   
   static inline
-  base_type norm_1(argument_type x) { return traits_type::abs(x); }
+  value_type abs(argument_type x) {
-
+    return std::abs(x);
-  static inline
+  }
-  base_type norm_2(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_inf(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  bool equals(argument_type lhs, argument_type rhs) { return lhs == rhs; }
 
   enum { is_complex = false };
 
@@ -224,57 +110,34 @@ struct NumericTraits<int> {
   };
 };
 
-
 /**
  * \class NumericTraits<float> NumericTraits.h "tvmet/NumericTraits.h"
  * \brief Traits specialized for float.
  */
 template<>
 struct NumericTraits<float> {
-  typedef float					value_type;
+  typedef float value_type;
-  typedef value_type				base_type;
+  typedef value_type real_type;
-  typedef double 				sum_type;
+  typedef value_type float_type;
-  typedef float 				diff_type;
+  typedef value_type argument_type;
-  typedef float 				float_type;
-  typedef float					signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef value_type				argument_type;
 
   static inline
-  base_type real(argument_type x) { return x; }
+  real_type real(argument_type x) { return x; }
 
   static inline
-  base_type imag(argument_type x) { TVMET_UNUSED(x); return 0; }
+  real_type imag(argument_type x) { TVMET_UNUSED(x); return 0; }
 
   static inline
   value_type conj(argument_type x) { return x; }
 
   static inline
-  base_type abs(argument_type x) { return std::abs(x); }
+  value_type sqrt(argument_type x) {
+    return std::sqrt(x);
+  }
 
   static inline
-  value_type sqrt(argument_type x) { return std::sqrt(x); }
+  value_type abs(argument_type x) {
-
+    return std::abs(x);
-  static inline
-  base_type norm_1(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_2(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_inf(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  bool equals(argument_type lhs, argument_type rhs) {
-    static base_type sqrt_epsilon(
-      NumericTraits<base_type>::sqrt(
-        std::numeric_limits<base_type>::epsilon()));
-
-    return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon *
-      std::max(std::max(traits_type::norm_inf(lhs),
-			traits_type::norm_inf(rhs)),
-	       std::numeric_limits<base_type>::min());
   }
 
   enum { is_complex = false };
@@ -293,50 +156,28 @@ struct NumericTraits<float> {
  */
 template<>
 struct NumericTraits<double> {
-  typedef double 				value_type;
+  typedef double value_type;
-  typedef value_type				base_type;
+  typedef value_type real_type;
-  typedef double 				sum_type;
+  typedef value_type float_type;
-  typedef double				diff_type;
+  typedef value_type argument_type;
-  typedef double 				float_type;
-  typedef double				signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef value_type				argument_type;
 
   static inline
-  base_type real(argument_type x) { return x; }
+  real_type real(argument_type x) { return x; }
 
   static inline
-  base_type imag(argument_type x) { TVMET_UNUSED(x); return 0; }
+  real_type imag(argument_type x) { TVMET_UNUSED(x); return 0; }
 
   static inline
   value_type conj(argument_type x) { return x; }
 
   static inline
-  base_type abs(argument_type x) { return std::abs(x); }
+  value_type sqrt(argument_type x) {
+    return std::sqrt(x);
+  }
 
   static inline
-  value_type sqrt(argument_type x) { return std::sqrt(x); }
+  value_type abs(argument_type x) {
-
+    return std::abs(x);
-  static inline
-  base_type norm_1(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_2(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  base_type norm_inf(argument_type x) { return traits_type::abs(x); }
-
-  static inline
-  bool equals(argument_type lhs, argument_type rhs) {
-    static base_type sqrt_epsilon(
-      NumericTraits<base_type>::sqrt(
-        std::numeric_limits<base_type>::epsilon()));
-
-    return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon *
-      std::max(std::max(traits_type::norm_inf(lhs),
-			traits_type::norm_inf(rhs)),
-	       std::numeric_limits<base_type>::min());
   }
 
   enum { is_complex = false };
@@ -360,75 +201,30 @@ struct NumericTraits<double> {
  */
 template<>
 struct NumericTraits< std::complex<int> > {
-  typedef int					base_type;
+  typedef std::complex<int> value_type;
-  typedef std::complex<int>			value_type;
+  typedef int real_type;
-  typedef std::complex<long> 			sum_type;
+  typedef std::complex<float> float_type;
-  typedef std::complex<int>			diff_type;
+  typedef const value_type& argument_type;
-  typedef std::complex<float>			float_type;
-  typedef std::complex<int>			signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef const value_type&			argument_type;
 
   static inline
-  base_type real(argument_type z) { return std::real(z); }
+  real_type real(argument_type z) { return std::real(z); }
 
   static inline
-  base_type imag(argument_type z) { return std::imag(z); }
+  real_type imag(argument_type z) { return std::imag(z); }
 
   static inline
   value_type conj(argument_type z) { return std::conj(z); }
 
   static inline
-  base_type abs(argument_type z) {
+  real_type abs(argument_type z) {
-    base_type x = z.real();
+    // the use of ceil() guarantees e.g. that abs(real(x)) <= abs(x),
-    base_type y = z.imag();
+    // and that abs(x)==0 if and only if x==0.
-
+    return static_cast<value_type>(std::ceil(std::abs(static_cast<float_type>(x))));
-    // XXX probably case of overrun; header complex uses scaling
-    return static_cast<base_type>(NumericTraits<base_type>::sqrt(x * x + y * y));
-  }
-
-  static /* inline */
-  value_type sqrt(argument_type z) {
-    // borrowed and adapted from header complex
-    base_type x = z.real();
-    base_type y = z.imag();
-
-    if(x == base_type()) {
-	base_type t = NumericTraits<base_type>::sqrt(
-                        NumericTraits<base_type>::abs(y) / 2);
-	return value_type(t, y < base_type() ? -t : t);
-    }
-    else {
-      base_type t = NumericTraits<base_type>::sqrt(
-		      2 * (traits_type::abs(z)
-		            + NumericTraits<base_type>::abs(x)));
-      base_type u = t / 2;
-      return x > base_type()
-	? value_type(u, y / t)
-	: value_type(NumericTraits<base_type>::abs(y) / t, y < base_type() ? -u : u);
-    }
   }
 
   static inline
-  base_type norm_1(argument_type z) {
+  value_type sqrt(argument_type x) {
-    return NumericTraits<base_type>::abs((traits_type::real(z)))
+    return static_cast<value_type>(std::sqrt(static_cast<float_type>(x)));
-         + NumericTraits<base_type>::abs((traits_type::imag(z)));
-  }
-
-  static inline
-  base_type norm_2(argument_type z) { return traits_type::abs(z); }
-
-  static inline
-  base_type norm_inf(argument_type z) {
-    return std::max(NumericTraits<base_type>::abs(traits_type::real(z)),
-		    NumericTraits<base_type>::abs(traits_type::imag(z)));
-  }
-
-  static inline
-  bool equals(argument_type lhs, argument_type rhs) {
-    return (traits_type::real(lhs) == traits_type::real(rhs))
-        && (traits_type::imag(lhs) == traits_type::imag(rhs));
   }
 
   enum { is_complex = true };
@@ -447,56 +243,28 @@ struct NumericTraits< std::complex<int> > {
  */
 template<>
 struct NumericTraits< std::complex<float> > {
-  typedef float					base_type;
+  typedef std::complex<float> value_type;
-  typedef std::complex<float>			value_type;
+  typedef float real_type;
-  typedef std::complex<double>			sum_type;
+  typedef value_type float_type;
-  typedef std::complex<float>			diff_type;
+  typedef const value_type& argument_type;
-  typedef std::complex<float>			float_type;
-  typedef std::complex<float>			signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef const value_type&			argument_type;
 
   static inline
-  base_type real(argument_type z) { return std::real(z); }
+  real_type real(argument_type z) { return std::real(z); }
 
   static inline
-  base_type imag(argument_type z) { return std::imag(z); }
+  real_type imag(argument_type z) { return std::imag(z); }
 
   static inline
   value_type conj(argument_type z) { return std::conj(z); }
 
   static inline
-  base_type abs(argument_type z) { return std::abs(z); }
+  value_type sqrt(argument_type x) {
-
+    return std::sqrt(x);
-  static inline
-  value_type sqrt(argument_type z) { return std::sqrt(z); }
-
-  static inline
-  base_type norm_1(argument_type z) {
-    return NumericTraits<base_type>::abs((traits_type::real(z)))
-         + NumericTraits<base_type>::abs((traits_type::imag(z)));
   }
 
   static inline
-  base_type norm_2(argument_type z) { return traits_type::abs(z); }
+  value_type abs(argument_type x) {
-
+    return std::abs(x);
-  static inline
-  base_type norm_inf(argument_type z) {
-    return std::max(NumericTraits<base_type>::abs(traits_type::real(z)),
-	            NumericTraits<base_type>::abs(traits_type::imag(z)));
-  }
-
- static inline
-  bool equals(argument_type lhs, argument_type rhs) {
-    static base_type sqrt_epsilon(
-      NumericTraits<base_type>::sqrt(
-        std::numeric_limits<base_type>::epsilon()));
-
-    return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon *
-      std::max(std::max(traits_type::norm_inf(lhs),
-			traits_type::norm_inf(rhs)),
-	       std::numeric_limits<base_type>::min());
   }
 
   enum { is_complex = true };
@@ -515,56 +283,28 @@ struct NumericTraits< std::complex<float> > {
  */
 template<>
 struct NumericTraits< std::complex<double> > {
-  typedef double				base_type;
+  typedef std::complex<double> value_type;
-  typedef std::complex<double>			value_type;
+  typedef double real_type;
-  typedef std::complex<double>			sum_type;
+  typedef value_type float_type;
-  typedef std::complex<double>			diff_type;
+  typedef const value_type& argument_type;
-  typedef std::complex<double>			float_type;
-  typedef std::complex<double>			signed_type;
-
-  typedef NumericTraits<value_type>		traits_type;
-  typedef const value_type&			argument_type;
 
   static inline
-  base_type real(argument_type z) { return std::real(z); }
+  real_type real(argument_type z) { return std::real(z); }
 
   static inline
-  base_type imag(argument_type z) { return std::imag(z); }
+  real_type imag(argument_type z) { return std::imag(z); }
 
   static inline
   value_type conj(argument_type z) { return std::conj(z); }
 
   static inline
-  base_type abs(argument_type z) { return std::abs(z); }
+  value_type sqrt(argument_type x) {
-
+    return std::sqrt(x);
-  static inline
-  value_type sqrt(argument_type z) { return std::sqrt(z); }
-
-  static inline
-  base_type norm_1(argument_type z) {
-    return NumericTraits<base_type>::abs((traits_type::real(z)))
-         + NumericTraits<base_type>::abs((traits_type::imag(z)));
   }
 
   static inline
-  base_type norm_2(argument_type z) { return traits_type::abs(z); }
+  value_type abs(argument_type x) {
-
+    return std::abs(x);
-  static inline
-  base_type norm_inf(argument_type z) {
-    return std::max(NumericTraits<base_type>::abs(traits_type::real(z)),
-	            NumericTraits<base_type>::abs(traits_type::imag(z)));
-  }
-
- static inline
-  bool equals(argument_type lhs, argument_type rhs) {
-    static base_type sqrt_epsilon(
-      NumericTraits<base_type>::sqrt(
-        std::numeric_limits<base_type>::epsilon()));
-
-    return traits_type::norm_inf(lhs - rhs) < sqrt_epsilon *
-      std::max(std::max(traits_type::norm_inf(lhs),
-			traits_type::norm_inf(rhs)),
-	       std::numeric_limits<base_type>::min());
   }
 
   enum { is_complex = true };
@@ -576,16 +316,8 @@ struct NumericTraits< std::complex<double> > {
   };
 };
 
-
 #endif // defined(EIGEN_USE_COMPLEX)
 
-
 } // namespace tvmet
 
-
 #endif //  TVMET_NUMERIC_TRAITS_H
-
-
-// Local Variables:
-// mode:C++
-// End:

tvmet-1.7.1/include/tvmet/util/Random.h

@@ -0,0 +1,88 @@
+/* This file is part of Eigen, a C++ template library for linear algebra
+ * Copyright (C) 2006-2007 Benoit Jacob <jacob@math.jussieu.fr>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *
+ * $Id: SelfTest.cc,v 1.1 2004/04/24 11:55:15 opetzold Exp $
+ */
+
+#ifndef TVMET_UTIL_RANDOM_H
+#define TVMET_UTIL_RANDOM_H
+
+#ifdef __GNUC__
+# if __GNUC__>=4
+#  define EIGEN_WITH_GCC_4_OR_LATER
+# endif
+#endif
+
+#include <cstdlib>
+
+#ifdef EIGEN_USE_COMPLEX
+#include <complex>
+#endif
+
+namespace tvmet {
+
+namespace util {
+
+/** Stores in x a random int between -RAND_MAX/2 and RAND_MAX/2 */
+inline void pickRandom( int & x )
+{
+    x = rand() - RAND_MAX / 2;
+}
+
+
+/** Stores in x a random float between -1.0 and 1.0 */
+inline void pickRandom( float & x )
+{
+    x = 2.0f * rand() / RAND_MAX - 1.0f;
+}
+
+/** Stores in x a random double between -1.0 and 1.0 */
+inline void pickRandom( double & x )
+{
+    x = 2.0 * rand() / RAND_MAX - 1.0;
+}
+
+#ifdef EIGEN_USE_COMPLEX
+/** Stores in the real and imaginary parts of x
+  * random values between -1.0 and 1.0 */
+template<typename T> void pickRandom( std::complex<T> & x )
+{
+#ifdef EIGEN_WITH_GCC_4_OR_LATER
+    pickRandom( x.real() );
+    pickRandom( x.imag() );
+#else // workaround by David Faure for MacOS 10.3 and GCC 3.3, commit 630812
+    T r = x.real();
+    T i = x.imag();
+    pickRandom( r );
+    pickRandom( i );
+    x = std::complex<T>(r,i);
+#endif
+}
+#endif // EIGEN_USE_COMPLEX
+
+template<typename T> T someRandom()
+{
+  T t;
+  pickRandom(t);
+  return t;
+}
+
+} // namespace util
+
+} // namespace tvmet
+
+#endif // TVMET_UTIL_RANDOM_H