fix linking issue with manage_caching_sizes_second_if_negative

(transplanted from 8dd3ae282d
)

Eigen/SVD

@@ -13,9 +13,9 @@ namespace Eigen {
   *
   *
   *
-  * This module provides SVD decomposition for (currently) real matrices.
+  * This module provides SVD decomposition for matrices (both real and complex).
   * This decomposition is accessible via the following MatrixBase method:
-  *  - MatrixBase::svd()
+  *  - MatrixBase::jacobiSvd()
   *
   * \code
   * #include <Eigen/SVD>

Eigen/src/Core/Block.h

@@ -94,7 +94,7 @@ struct traits<Block<XprType, BlockRows, BlockCols, InnerPanel, HasDirectAccess>
     MaskPacketAccessBit = (InnerSize == Dynamic || (InnerSize % packet_traits<Scalar>::size) == 0)
                        && (InnerStrideAtCompileTime == 1)
                         ? PacketAccessBit : 0,
-    MaskAlignedBit = (InnerPanel && (OuterStrideAtCompileTime!=Dynamic) && (((OuterStrideAtCompileTime * sizeof(Scalar)) % 16) == 0)) ? AlignedBit : 0,
+    MaskAlignedBit = (InnerPanel && (OuterStrideAtCompileTime!=Dynamic) && (((OuterStrideAtCompileTime * int(sizeof(Scalar))) % 16) == 0)) ? AlignedBit : 0,
     FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1) ? LinearAccessBit : 0,
     FlagsLvalueBit = is_lvalue<XprType>::value ? LvalueBit : 0,
     FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0,

Eigen/src/Core/Map.h

@@ -102,7 +102,7 @@ struct traits<Map<PlainObjectType, MapOptions, StrideType> >
                            || HasNoOuterStride
                            || ( OuterStrideAtCompileTime!=Dynamic
                            && ((static_cast<int>(sizeof(Scalar))*OuterStrideAtCompileTime)%16)==0 ) ),
-    Flags0 = TraitsBase::Flags,
+    Flags0 = TraitsBase::Flags & (~NestByRefBit),
     Flags1 = IsAligned ? (int(Flags0) | AlignedBit) : (int(Flags0) & ~AlignedBit),
     Flags2 = (bool(HasNoStride) || bool(PlainObjectType::IsVectorAtCompileTime))
            ? int(Flags1) : int(Flags1 & ~LinearAccessBit),
@@ -120,7 +120,6 @@ template<typename PlainObjectType, int MapOptions, typename StrideType> class Ma
   public:
 
     typedef MapBase<Map> Base;
-
     EIGEN_DENSE_PUBLIC_INTERFACE(Map)
 
     typedef typename Base::PointerType PointerType;
@@ -181,7 +180,6 @@ template<typename PlainObjectType, int MapOptions, typename StrideType> class Ma
       PlainObjectType::Base::_check_template_params();
     }
 
-
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Map)
 
   protected:

Eigen/src/Core/products/GeneralBlockPanelKernel.h

@@ -30,19 +30,16 @@ namespace internal {
 template<typename _LhsScalar, typename _RhsScalar, bool _ConjLhs=false, bool _ConjRhs=false>
 class gebp_traits;
 
+inline std::ptrdiff_t manage_caching_sizes_second_if_negative(std::ptrdiff_t a, std::ptrdiff_t b)
+{
+  return a<=0 ? b : a;
+}
+
 /** \internal */
 inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1=0, std::ptrdiff_t* l2=0)
 {
-  static std::ptrdiff_t m_l1CacheSize = 0;
-  static std::ptrdiff_t m_l2CacheSize = 0;
-  if(m_l1CacheSize==0)
-  {
-    m_l1CacheSize = queryL1CacheSize();
-    m_l2CacheSize = queryTopLevelCacheSize();
-
-    if(m_l1CacheSize<=0) m_l1CacheSize = 8 * 1024;
-    if(m_l2CacheSize<=0) m_l2CacheSize = 1 * 1024 * 1024;
-  }
+  static std::ptrdiff_t m_l1CacheSize = manage_caching_sizes_second_if_negative(queryL1CacheSize(),8 * 1024);
+  static std::ptrdiff_t m_l2CacheSize = manage_caching_sizes_second_if_negative(queryTopLevelCacheSize(),1*1024*1024);
 
   if(action==SetAction)
   {

Eigen/src/Core/util/Macros.h

@@ -1,3 +1,4 @@
+
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
@@ -28,7 +29,7 @@
 
 #define EIGEN_WORLD_VERSION 3
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 4
+#define EIGEN_MINOR_VERSION 5
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \

Eigen/src/Core/util/XprHelper.h

@@ -127,7 +127,7 @@ class compute_matrix_flags
             ((Options&DontAlign)==0)
         && (
 #if EIGEN_ALIGN_STATICALLY
-             ((!is_dynamic_size_storage) && (((MaxCols*MaxRows*sizeof(Scalar)) % 16) == 0))
+             ((!is_dynamic_size_storage) && (((MaxCols*MaxRows*int(sizeof(Scalar))) % 16) == 0))
 #else
              0
 #endif

Eigen/src/Eigenvalues/EigenSolver.h

@@ -291,7 +291,7 @@ template<typename _MatrixType> class EigenSolver
 
     ComputationInfo info() const
     {
-      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+      eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
       return m_realSchur.info();
     }
 
@@ -339,7 +339,7 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
   EigenvectorsType matV(n,n);
   for (Index j=0; j<n; ++j)
   {
-    if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))))
+    if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n)
     {
       // we have a real eigen value
       matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
@@ -570,10 +570,13 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
 
         }
       }
+      
+      // We handled a pair of complex conjugate eigenvalues, so need to skip them both
+      n--;
     }
     else
     {
-      eigen_assert("Internal bug in EigenSolver"); // this should not happen
+      eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
     }
   }
 

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -220,7 +220,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     const MatrixType& eigenvectors() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec;
     }
@@ -243,7 +242,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     const RealVectorType& eigenvalues() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       return m_eivalues;
     }
 
@@ -268,7 +266,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     MatrixType operatorSqrt() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
     }
@@ -294,7 +291,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     MatrixType operatorInverseSqrt() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
-      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
     }

Eigen/src/Geometry/Hyperplane.h

@@ -225,7 +225,7 @@ public:
       normal() = mat * normal();
     else
     {
-      eigen_assert("invalid traits value in Hyperplane::transform()");
+      eigen_assert(0 && "invalid traits value in Hyperplane::transform()");
     }
     return *this;
   }

Eigen/src/Geometry/Quaternion.h

@@ -682,7 +682,7 @@ QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& oth
   Scalar scale0;
   Scalar scale1;
 
-  if (absD>=one)
+  if(absD>=one)
   {
     scale0 = Scalar(1) - t;
     scale1 = t;
@@ -695,9 +695,8 @@ QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& oth
 
     scale0 = internal::sin( ( Scalar(1) - t ) * theta) / sinTheta;
     scale1 = internal::sin( ( t * theta) ) / sinTheta;
-    if (d<0)
-      scale1 = -scale1;
   }
+  if(d<0) scale1 = -scale1;
 
   return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs());
 }

Eigen/src/Geometry/Rotation2D.h

@@ -89,7 +89,7 @@ public:
 
   /** Concatenates two rotations */
   inline Rotation2D& operator*=(const Rotation2D& other)
-  { return m_angle += other.m_angle; return *this; }
+  { m_angle += other.m_angle; return *this; }
 
   /** Applies the rotation to a 2D vector */
   Vector2 operator* (const Vector2& vec) const

Eigen/src/LU/arch/Inverse_SSE.h

@@ -55,7 +55,7 @@ struct compute_inverse_size4<Architecture::SSE, float, MatrixType, ResultType>
   
   static void run(const MatrixType& matrix, ResultType& result)
   {
-    EIGEN_ALIGN16 const  int _Sign_PNNP[4] = { 0x00000000, 0x80000000, 0x80000000, 0x00000000 };
+    EIGEN_ALIGN16 const unsigned int _Sign_PNNP[4] = { 0x00000000, 0x80000000, 0x80000000, 0x00000000 };
 
     // Load the full matrix into registers
     __m128 _L1 = matrix.template packet<MatrixAlignment>( 0);

Eigen/src/SVD/JacobiSVD.h

@@ -708,6 +708,13 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
 };
 } // end namespace internal
 
+/** \svd_module
+  *
+  * \return the singular value decomposition of \c *this computed by two-sided
+  * Jacobi transformations.
+  *
+  * \sa class JacobiSVD
+  */
 template<typename Derived>
 JacobiSVD<typename MatrixBase<Derived>::PlainObject>
 MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const

doc/C09_TutorialSparse.dox

@@ -4,7 +4,7 @@ namespace Eigen {
     \ingroup Tutorial
 
 \li \b Previous: \ref TutorialGeometry
-\li \b Next: TODO
+\li \b Next: \ref TutorialMapClass
 
 \b Table \b of \b contents \n
   - \ref TutorialSparseIntro
@@ -254,7 +254,7 @@ if(!lu_of_A.solve(b,&x)) {
 
 See also the class SparseLLT, class SparseLU, and class SparseLDLT.
 
-\li \b Next: TODO
+\li \b Next: \ref TutorialMapClass
 
 */
 

doc/C10_TutorialMapClass.dox

@@ -0,0 +1,96 @@
+namespace Eigen {
+
+/** \page TutorialMapClass Tutorial page 10 - Interfacing with C/C++ arrays and external libraries: the %Map class
+
+\ingroup Tutorial
+
+\li \b Previous: \ref TutorialSparse
+\li \b Next: \ref TODO
+
+This tutorial page explains how to work with "raw" C++ arrays.  This can be useful in a variety of contexts, particularly when "importing" vectors and matrices from other libraries into Eigen.
+
+\b Table \b of \b contents
+  - \ref TutorialMapIntroduction
+  - \ref TutorialMapTypes
+  - \ref TutorialMapUsing
+  - \ref TutorialMapPlacementNew
+
+\section TutorialMapIntroduction Introduction
+
+Occasionally you may have a pre-defined array of numbers that you want to use within Eigen as a vector or matrix. While one option is to make a copy of the data, most commonly you probably want to re-use this memory as an Eigen type. Fortunately, this is very easy with the Map class.
+
+\section TutorialMapTypes Map types and declaring Map variables
+
+A Map object has a type defined by its Eigen equivalent:
+\code
+Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >
+\endcode
+Note that, in this default case, a Map requires just a single template parameter.  
+
+To construct a Map variable, you need two other pieces of information: a pointer to the region of memory defining the array of coefficients, and the desired shape of the matrix or vector.  For example, to define a matrix of \c float with sizes determined at compile time, you might do the following:
+\code
+Map<MatrixXf> mf(pf,rows,columns);
+\endcode
+where \c pf is a \c float \c * pointing to the array of memory.  A fixed-size read-only vector of integers might be declared as
+\code
+Map<const Vector4i> mi(pi);
+\endcode
+where \c pi is an \c int \c *. In this case the size does not have to be passed to the constructor, because it is already specified by the Matrix/Array type.
+
+Note that Map does not have a default constructor; you \em must pass a pointer to intialize the object. However, you can work around this requirement (see \ref TutorialMapPlacementNew).
+
+Map is flexible enough to accomodate a variety of different data representations.  There are two other (optional) template parameters:
+\code
+Map<typename MatrixType,
+    int MapOptions,
+    typename StrideType>
+\endcode
+\li \c MapOptions specifies whether the pointer is \c #Aligned, or \c #Unaligned.  The default is \c #Unaligned.
+\li \c StrideType allows you to specify a custom layout for the memory array, using the Stride class.  One example would be to specify that the data array is organized in row-major format:
+<table class="example">
+<tr><th>Example:</th><th>Output:</th></tr>
+<tr>
+<td>\include Tutorial_Map_rowmajor.cpp </td>
+<td>\verbinclude Tutorial_Map_rowmajor.out </td>
+</table>
+However, Stride is even more flexible than this; for details, see the documentation for the Map and Stride classes.
+
+\section TutorialMapUsing Using Map variables
+
+You can use a Map object just like any other Eigen type:
+<table class="example">
+<tr><th>Example:</th><th>Output:</th></tr>
+<tr>
+<td>\include Tutorial_Map_using.cpp </td>
+<td>\verbinclude Tutorial_Map_using.out </td>
+</table>
+
+However, when writing functions taking Eigen types, it is important to realize that a Map type is \em not identical to its Dense equivalent.  See \ref TopicFunctionTakingEigenTypesMultiarguments for details.
+
+\section TutorialMapPlacementNew Changing the mapped array
+
+It is possible to change the array of a Map object after declaration, using the C++ "placement new" syntax:
+<table class="example">
+<tr><th>Example:</th><th>Output:</th></tr>
+<tr>
+<td>\include Map_placement_new.cpp </td>
+<td>\verbinclude Map_placement_new.out </td>
+</table>
+Despite appearances, this does not invoke the memory allocator, because the syntax specifies the location for storing the result.
+
+This syntax makes it possible to declare a Map object without first knowing the mapped array's location in memory:
+\code
+Map<Matrix3f> A(NULL);  // don't try to use this matrix yet!
+VectorXf b(n_matrices);
+for (int i = 0; i < n_matrices; i++)
+{
+  new (&A) Map<Matrix3f>(get_matrix_pointer(i));
+  b(i) = A.trace();
+}
+\endcode
+
+\li \b Next: \ref TODO
+
+*/
+
+}

doc/I13_FunctionsTakingEigenTypes.dox

@@ -47,9 +47,9 @@ void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
 Prints the maximum coefficient of the array or array-expression.
 \code
 template <typename Derived>
-void print_max(const ArrayBase<Derived>& a, const ArrayBase<Derived>& b)
+void print_max_coeff(const ArrayBase<Derived> &a)
 {
-  std::cout << "max: " << (a.max(b)).maxCoeff() << std::endl;
+  std::cout << "max: " << a.maxCoeff() << std::endl;
 }
 \endcode
 <b> %MatrixBase Example </b><br/><br/>
@@ -63,6 +63,21 @@ void print_inv_cond(const MatrixBase<Derived>& a)
   std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
 }
 \endcode
+<b> Multiple templated arguments example </b><br/><br/>
+Calculate the Euclidean distance between two points.
+\code
+template <typename DerivedA,typename DerivedB>
+typename DerivedA::Scalar squaredist(const MatrixBase<DerivedA>& p1,const MatrixBase<DerivedB>& p2)
+{
+  return (p1-p2).squaredNorm();
+}
+\endcode
+Notice that we used two template parameters, one per argument. This permits the function to handle inputs of different types, e.g.,
+\code
+squaredist(v1,2*v2)
+\endcode
+where the first argument \c v1 is a vector and the second argument \c 2*v2 is an expression.
+<br/><br/>
 
 These examples are just intended to give the reader a first impression of how functions can be written which take a plain and constant Matrix or Array argument. They are also intended to give the reader an idea about the most common base classes being the optimal candidates for functions. In the next section we will look in more detail at an example and the different ways it can be implemented, while discussing each implementation's problems and advantages. For the discussion below, Matrix and Array as well as MatrixBase and ArrayBase can be exchanged and all arguments still hold.
 
@@ -128,6 +143,8 @@ The implementation above does now not only work with temporary expressions but i
 
 \b Note: The const cast hack will only work with templated functions. It will not work with the MatrixXf implementation because it is not possible to cast a Block expression to a Matrix reference!
 
+
+
 \section TopicResizingInGenericImplementations How to resize matrices in generic implementations?
 
 One might think we are done now, right? This is not completely true because in order for our covariance function to be generically applicable, we want the follwing code to work

doc/Overview.dox

@@ -27,6 +27,7 @@ For a first contact with Eigen, the best place is to have a look at the \ref Get
     - \ref TutorialReductionsVisitorsBroadcasting
     - \ref TutorialGeometry
     - \ref TutorialSparse
+    - \ref TutorialMapClass
   - \ref QuickRefPage
   - <b>Advanced topics</b>
     - \ref TopicAliasing

doc/QuickReference.dox

@@ -645,11 +645,11 @@ m3 -= s1 * m3.adjoint() * m1.selfadjointView<Eigen::Lower>();\endcode
 </td></tr>
 <tr><td>
 Rank 1 and rank K update: \n
-\f$ upper(M_1) \mathrel{{+}{=}} s_1 M_2^* M_2 \f$ \n
-\f$ lower(M_1) \mathbin{{-}{=}} M_2 M_2^* \f$
+\f$ upper(M_1) \mathrel{{+}{=}} s_1 M_2 M_2^* \f$ \n
+\f$ lower(M_1) \mathbin{{-}{=}} M_2^* M_2 \f$
 </td><td>\n \code
 M1.selfadjointView<Eigen::Upper>().rankUpdate(M2,s1);
-m1.selfadjointView<Eigen::Lower>().rankUpdate(m2.adjoint(),-1); \endcode
+M1.selfadjointView<Eigen::Lower>().rankUpdate(M2.adjoint(),-1); \endcode
 </td></tr>
 <tr><td>
 Rank 2 update: (\f$ M \mathrel{{+}{=}} s u v^* + s v u^* \f$)

doc/snippets/Tutorial_Map_rowmajor.cpp

@@ -0,0 +1,7 @@
+int array[8];
+for(int i = 0; i < 8; ++i) array[i] = i;
+cout << "Column-major:\n" << Map<Matrix<int,2,4> >(array) << endl;
+cout << "Row-major:\n" << Map<Matrix<int,2,4,RowMajor> >(array) << endl;
+cout << "Row-major using stride:\n" <<
+  Map<Matrix<int,2,4>, Unaligned, Stride<1,4> >(array) << endl;
+

doc/snippets/Tutorial_Map_using.cpp

@@ -0,0 +1,21 @@
+typedef Matrix<float,1,Dynamic> MatrixType;
+typedef Map<MatrixType> MapType;
+typedef Map<const MatrixType> MapTypeConst;   // a read-only map
+const int n_dims = 5;
+  
+MatrixType m1(n_dims), m2(n_dims);
+m1.setRandom();
+m2.setRandom();
+float *p = &m2(0);  // get the address storing the data for m2
+MapType m2map(p,m2.size());   // m2map shares data with m2
+MapTypeConst m2mapconst(p,m2.size());  // a read-only accessor for m2
+
+cout << "m1: " << m1 << endl;
+cout << "m2: " << m2 << endl;
+cout << "Squared euclidean distance: " << (m1-m2).squaredNorm() << endl;
+cout << "Squared euclidean distance, using map: " <<
+  (m1-m2map).squaredNorm() << endl;
+m2map(3) = 7;   // this will change m2, since they share the same array
+cout << "Updated m2: " << m2 << endl;
+cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << endl;
+/* m2mapconst(2) = 5; */   // this yields a compile-time error

scripts/eigen_gen_docs

@@ -15,7 +15,7 @@ mkdir build -p
 
 #step 2 : upload
 # (the '/' at the end of path are very important, see rsync documentation)
-rsync -az build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-3.0/ || { echo "upload failed"; exit 1; }
+rsync -az --no-p build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-3.0/ || { echo "upload failed"; exit 1; }
 
 echo "Uploaded successfully"
 

test/geo_quaternion.cpp

@@ -28,6 +28,35 @@
 #include <Eigen/LU>
 #include <Eigen/SVD>
 
+template<typename T> T bounded_acos(T v)
+{
+  using std::acos;
+  using std::min;
+  using std::max;
+  return acos((max)(T(-1),(min)(v,T(1))));
+}
+
+template<typename QuatType> void check_slerp(const QuatType& q0, const QuatType& q1)
+{
+  typedef typename QuatType::Scalar Scalar;
+  typedef Matrix<Scalar,3,1> VectorType;
+  typedef AngleAxis<Scalar> AA;
+
+  Scalar largeEps = test_precision<Scalar>();
+
+  Scalar theta_tot = AA(q1*q0.inverse()).angle();
+  if(theta_tot>M_PI)
+    theta_tot = 2.*M_PI-theta_tot;
+  for(Scalar t=0; t<=1.001; t+=0.1)
+  {
+    QuatType q = q0.slerp(t,q1);
+    Scalar theta = AA(q*q0.inverse()).angle();
+    VERIFY(internal::abs(q.norm() - 1) < largeEps);
+    if(theta_tot==0)  VERIFY(theta_tot==0);
+    else              VERIFY(internal::abs(theta/theta_tot - t) < largeEps);
+  }
+}
+
 template<typename Scalar, int Options> void quaternion(void)
 {
   /* this test covers the following files:
@@ -36,6 +65,7 @@ template<typename Scalar, int Options> void quaternion(void)
 
   typedef Matrix<Scalar,3,3> Matrix3;
   typedef Matrix<Scalar,3,1> Vector3;
+  typedef Matrix<Scalar,4,1> Vector4;
   typedef Quaternion<Scalar,Options> Quaternionx;
   typedef AngleAxis<Scalar> AngleAxisx;
 
@@ -50,7 +80,8 @@ template<typename Scalar, int Options> void quaternion(void)
           v2 = Vector3::Random(),
           v3 = Vector3::Random();
 
-  Scalar a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI));
+  Scalar  a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI)),
+          b = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI));
 
   // Quaternion: Identity(), setIdentity();
   Quaternionx q1, q2;
@@ -124,6 +155,22 @@ template<typename Scalar, int Options> void quaternion(void)
   // test bug 369 - improper alignment.
   Quaternionx *q = new Quaternionx;
   delete q;
+
+  q1 = AngleAxisx(a, v0.normalized());
+  q2 = AngleAxisx(b, v1.normalized());
+  check_slerp(q1,q2);
+
+  q1 = AngleAxisx(b, v1.normalized());
+  q2 = AngleAxisx(b+M_PI, v1.normalized());
+  check_slerp(q1,q2);
+
+  q1 = AngleAxisx(b,  v1.normalized());
+  q2 = AngleAxisx(-b, -v1.normalized());
+  check_slerp(q1,q2);
+
+  q1.coeffs() = Vector4::Random().normalized();
+  q2.coeffs() = -q1.coeffs();
+  check_slerp(q1,q2);
 }
 
 template<typename Scalar> void mapQuaternion(void){