* fix issues with "long double" type (useful to enforce the use of x87 registers)

* extend the documentation on "extending Eigen"

5 files changed, 152 insertions, 14 deletions

diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index a523db4c3..1a6d8014e 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -34,6 +34,10 @@ template<typename T> inline T ei_random_amplitude()
   else return static_cast<T>(10);
 }
 
+/**************
+***   int   ***
+**************/
+
 template<> inline int precision<int>() { return 0; }
 inline int ei_real(int x)  { return x; }
 inline int ei_imag(int)    { return 0; }
@@ -74,6 +78,10 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision<int>())
   return a <= b;
 }
 
+/**************
+*** float   ***
+**************/
+
 template<> inline float precision<float>() { return 1e-5f; }
 inline float ei_real(float x)  { return x; }
 inline float ei_imag(float)    { return 0.f; }
@@ -81,10 +89,10 @@ inline float ei_conj(float x)  { return x; }
 inline float ei_abs(float x)   { return std::abs(x); }
 inline float ei_abs2(float x)  { return x*x; }
 inline float ei_sqrt(float x)  { return std::sqrt(x); }
-inline float ei_exp(float x)  { return std::exp(x); }
-inline float ei_log(float x)  { return std::log(x); }
-inline float ei_sin(float x)  { return std::sin(x); }
-inline float ei_cos(float x)  { return std::cos(x); }
+inline float ei_exp(float x)   { return std::exp(x); }
+inline float ei_log(float x)   { return std::log(x); }
+inline float ei_sin(float x)   { return std::sin(x); }
+inline float ei_cos(float x)   { return std::cos(x); }
 inline float ei_pow(float x, float y)  { return std::pow(x, y); }
 
 template<> inline float ei_random(float a, float b)
@@ -115,6 +123,10 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float
   return a <= b || ei_isApprox(a, b, prec);
 }
 
+/**************
+*** double  ***
+**************/
+
 template<> inline double precision<double>() { return 1e-11; }
 inline double ei_real(double x)  { return x; }
 inline double ei_imag(double)    { return 0.; }
@@ -122,10 +134,10 @@ inline double ei_conj(double x)  { return x; }
 inline double ei_abs(double x)   { return std::abs(x); }
 inline double ei_abs2(double x)  { return x*x; }
 inline double ei_sqrt(double x)  { return std::sqrt(x); }
-inline double ei_exp(double x)  { return std::exp(x); }
-inline double ei_log(double x)  { return std::log(x); }
-inline double ei_sin(double x)  { return std::sin(x); }
-inline double ei_cos(double x)  { return std::cos(x); }
+inline double ei_exp(double x)   { return std::exp(x); }
+inline double ei_log(double x)   { return std::log(x); }
+inline double ei_sin(double x)   { return std::sin(x); }
+inline double ei_cos(double x)   { return std::cos(x); }
 inline double ei_pow(double x, double y)  { return std::pow(x, y); }
 
 template<> inline double ei_random(double a, double b)
@@ -156,6 +168,10 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision<do
   return a <= b || ei_isApprox(a, b, prec);
 }
 
+/*********************
+*** complex<float> ***
+*********************/
+
 template<> inline float precision<std::complex<float> >() { return precision<float>(); }
 inline float ei_real(const std::complex<float>& x) { return std::real(x); }
 inline float ei_imag(const std::complex<float>& x) { return std::imag(x); }
@@ -185,6 +201,10 @@ inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>&
 }
 // ei_isApproxOrLessThan wouldn't make sense for complex numbers
 
+/**********************
+*** complex<double> ***
+**********************/
+
 template<> inline double precision<std::complex<double> >() { return precision<double>(); }
 inline double ei_real(const std::complex<double>& x) { return std::real(x); }
 inline double ei_imag(const std::complex<double>& x) { return std::imag(x); }
@@ -214,4 +234,43 @@ inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double
 }
 // ei_isApproxOrLessThan wouldn't make sense for complex numbers
 
+
+/******************
+*** long double ***
+******************/
+
+template<> inline long double precision<long double>() { return precision<double>(); }
+inline long double ei_real(long double x)  { return x; }
+inline long double ei_imag(long double)    { return 0.; }
+inline long double ei_conj(long double x)  { return x; }
+inline long double ei_abs(long double x)   { return std::abs(x); }
+inline long double ei_abs2(long double x)  { return x*x; }
+inline long double ei_sqrt(long double x)  { return std::sqrt(x); }
+inline long double ei_exp(long double x)   { return std::exp(x); }
+inline long double ei_log(long double x)   { return std::log(x); }
+inline long double ei_sin(long double x)   { return std::sin(x); }
+inline long double ei_cos(long double x)   { return std::cos(x); }
+inline long double ei_pow(long double x, long double y)  { return std::pow(x, y); }
+
+template<> inline long double ei_random(long double a, long double b)
+{
+  return ei_random<double>(a,b);
+}
+template<> inline long double ei_random()
+{
+  return ei_random<double>(-ei_random_amplitude<double>(), ei_random_amplitude<double>());
+}
+inline bool ei_isMuchSmallerThan(long double a, long double b, long double prec = precision<long double>())
+{
+  return ei_abs(a) <= ei_abs(b) * prec;
+}
+inline bool ei_isApprox(long double a, long double b, long double prec = precision<long double>())
+{
+  return ei_abs(a - b) <= std::min(ei_abs(a), ei_abs(b)) * prec;
+}
+inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec = precision<long double>())
+{
+  return a <= b || ei_isApprox(a, b, prec);
+}
+
 #endif // EIGEN_MATHFUNCTIONS_H
diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h
index b1ac341db..842c7bce2 100644
--- a/Eigen/src/Core/NumTraits.h
+++ b/Eigen/src/Core/NumTraits.h
@@ -31,7 +31,8 @@
   *
   * \param T the numeric type about which this class provides data. Recall that Eigen allows
   *          only the following types for \a T: \c int, \c float, \c double,
-  *          \c std::complex<float>, \c std::complex<double>.
+  *          \c std::complex<float>, \c std::complex<double>, and \c long \c double (especially
+  *          useful to enforce x87 arithmetics when SSE is the default).
   *
   * The provided data consists of:
   * \li A typedef \a Real, giving the "real part" type of \a T. If \a T is already real,
@@ -120,8 +121,8 @@ template<> struct NumTraits<long double>
     IsComplex = 0,
     HasFloatingPoint = 1,
     ReadCost = 1,
-    AddCost = 2,
-    MulCost = 2
+    AddCost = 1,
+    MulCost = 1
   };
 };
 
diff --git a/doc/CustomizingEigen.dox b/doc/CustomizingEigen.dox
index 75981b891..b5cdf144f 100644
--- a/doc/CustomizingEigen.dox
+++ b/doc/CustomizingEigen.dox
@@ -1,10 +1,13 @@
 namespace Eigen {
 
-/** \page CustomizingEigen
+/** \page CustomizingEigen Customizing/Extending Eigen
 
-<h1>Customizing Eigen</h1>
+Eigen2 can be extended in several way, for instance, by defining global methods, \ref ExtendingMatrixBase "by adding custom methods to MatrixBase", adding support to \ref CustomScalarType "custom types" etc.
 
-Eigen2 can be extended in several way, for instance, by defining global methods, \ref ExtendingMatrixBase "by adding custom methods to MatrixBase", etc.
+\b Table \b of \b contents
+  - \ref ExtendingMatrixBase
+  - \ref CustomScalarType
+  - \ref PreprocessorDirectives
 
 \section ExtendingMatrixBase Extending MatrixBase
 
@@ -66,11 +69,77 @@ Then one can the following declaration in the config.h or whatever prerequisites
 #define EIGEN_MATRIXBASE_PLUGIN "MatrixBaseAddons.h"
 \endcode
 
+
+
+\section CustomScalarType Using custom scalar types
+
+By default, Eigen currently supports the following scalar types: \c int, \c float, \c double, \c std::complex<float>, \c std::complex<double>, \c long \c double, \c long \c long \c int (64 bits integers), and \c bool. The \c long \c double is especially useful on x86-64 systems or when the SSE2 instruction set is enabled because it enforces the use of x87 registers with extended accuracy.
+
+In order to add support for a custom type \c T you need:
+ 1 - make sure the common operator (+,-,*,/,etc.) are supported by the type \c T
+ 2 - add a specialization of struct Eigen::NumTraits<T> (see \ref class NumTraits)
+ 3 - define a couple of math functions for your type such as: ei_sqrt, ei_abs, etc...
+     (see the file Eigen/src/Core/MathFunctions.h)
+
+Here is a concrete example adding support for the Adolc's \c adouble type. <a href="http://www.math.tu-dresden.de/~adol-c/">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives.
+
+\code
+#ifndef ADLOCSUPPORT_H
+#define ADLOCSUPPORT_H
+
+#define ADOLC_TAPELESS
+#include <adolc/adouble.h>
+#include <Eigen/Core>
+
+namespace Eigen {
+
+template<> struct NumTraits<adtl::adouble>
+{
+  typedef adtl::adouble Real;
+  typedef adtl::adouble FloatingPoint;
+  enum {
+    IsComplex = 0,
+    HasFloatingPoint = 1,
+    ReadCost = 1,
+    AddCost = 1,
+    MulCost = 1
+  };
+};
+
+}
+
+// the Adolc's type adouble is defined in the adtl namespace
+// therefore, the following ei_* functions *must* be defined
+// in the same namespace
+namespace adtl {
+
+  inline const adouble& ei_conj(const adouble& x)  { return x; }
+  inline const adouble& ei_real(const adouble& x)  { return x; }
+  inline adouble ei_imag(const adouble&)    { return 0.; }
+  inline adouble ei_abs(const adouble&  x)  { return fabs(x); }
+  inline adouble ei_abs2(const adouble& x)  { return x*x; }
+  inline adouble ei_sqrt(const adouble& x)  { return sqrt(x); }
+  inline adouble ei_exp(const adouble&  x)  { return exp(x); }
+  inline adouble ei_log(const adouble&  x)  { return log(x); }
+  inline adouble ei_sin(const adouble&  x)  { return sin(x); }
+  inline adouble ei_cos(const adouble&  x)  { return cos(x); }
+  inline adouble ei_pow(const adouble& x, adouble y)  { return pow(x, y); }
+
+}
+
+#endif // ADLOCSUPPORT_H
+\endcode
+
+
+
 \section PreprocessorDirectives Preprocessor directives
 
+The value of the following preprocessor tokens can be overwritten by defining them before including any Eigen's headers.
  - \b EIGEN_DONT_VECTORIZE disables explicit vectorization when defined.
  - \b EIGEN_UNROLLING_LIMIT defines the maximal instruction counts to enable meta unrolling of loops. Set it to zero to disable unrolling. The default is 100.
  - \b EIGEN_TUNE_FOR_L2_CACHE_SIZE represents the maximal size in Bytes of L2 blocks. Since several blocks have to stay concurently in L2 cache, this value should correspond to at most 1/4 of the size of L2 cache.
+ - \b EIGEN_NO_STATIC_ASSERT replaces compile time static assertions by runtime assertions
+ - \b EIGEN_MATRIXBASE_PLUGIN see \ref ExtendingMatrixBase
 
 */
 
diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp
index 769e021af..40e999c87 100644
--- a/test/basicstuff.cpp
+++ b/test/basicstuff.cpp
@@ -111,5 +111,6 @@ void test_basicstuff()
     CALL_SUBTEST( basicStuff(MatrixXi(8, 12)) );
     CALL_SUBTEST( basicStuff(MatrixXcd(20, 20)) );
     CALL_SUBTEST( basicStuff(Matrix<float, 100, 100>()) );
+    CALL_SUBTEST( basicStuff(Matrix<long double,Dynamic,Dynamic>(10,10)) );
   }
 }
diff --git a/test/main.h b/test/main.h
index 5079fdb5b..aeba4429d 100644
--- a/test/main.h
+++ b/test/main.h
@@ -169,6 +169,7 @@ template<> inline float test_precision<float>() { return 1e-4f; }
 template<> inline double test_precision<double>() { return 1e-6; }
 template<> inline float test_precision<std::complex<float> >() { return test_precision<float>(); }
 template<> inline double test_precision<std::complex<double> >() { return test_precision<double>(); }
+template<> inline long double test_precision<long double>() { return 1e-6; }
 
 inline bool test_ei_isApprox(const int& a, const int& b)
 { return ei_isApprox(a, b, test_precision<int>()); }
@@ -201,6 +202,13 @@ inline bool test_ei_isApprox(const std::complex<double>& a, const std::complex<d
 inline bool test_ei_isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b)
 { return ei_isMuchSmallerThan(a, b, test_precision<std::complex<double> >()); }
 
+inline bool test_ei_isApprox(const long double& a, const long double& b)
+{ return ei_isApprox(a, b, test_precision<long double>()); }
+inline bool test_ei_isMuchSmallerThan(const long double& a, const long double& b)
+{ return ei_isMuchSmallerThan(a, b, test_precision<long double>()); }
+inline bool test_ei_isApproxOrLessThan(const long double& a, const long double& b)
+{ return ei_isApproxOrLessThan(a, b, test_precision<long double>()); }
+
 template<typename Derived1, typename Derived2>
 inline bool test_ei_isApprox(const MatrixBase<Derived1>& m1,
                    const MatrixBase<Derived2>& m2)