bug #698: rewrite LinSpaced for integer scalar types to avoid overflow and guarantee an even spacing when possible.

Otherwise, the "high" bound is implicitly lowered to the largest value allowing for an even distribution.
This changeset also disable vectorization for this integer path.

4 files changed, 62 insertions, 31 deletions

diff --git a/Eigen/src/Core/CwiseNullaryOp.h b/Eigen/src/Core/CwiseNullaryOp.h
index 1caf37bb7..757320ace 100644
--- a/Eigen/src/Core/CwiseNullaryOp.h
+++ b/Eigen/src/Core/CwiseNullaryOp.h
@@ -264,6 +264,16 @@ DenseBase<Derived>::LinSpaced(Sequential_t, const Scalar& low, const Scalar& hig
   * Example: \include DenseBase_LinSpaced.cpp
   * Output: \verbinclude DenseBase_LinSpaced.out
   *
+  * For integer scalar types, an even spacing is possible if and only if the length of the range,
+  * i.e., \c high-low is a scalar multiple of \c size-1, or if \c size is a scalar multiple of the
+  * number of values \c high-low+1 (meaning each value can be repeated the same number of time).
+  * If one of these two considions is not satisfied, then \c high is lowered to the largest value
+  * satisfying one of this constraint.
+  * Here are some examples:
+  *
+  * Example: \include DenseBase_LinSpacedInt.cpp
+  * Output: \verbinclude DenseBase_LinSpacedInt.out
+  *
   * \sa setLinSpaced(Index,const Scalar&,const Scalar&), LinSpaced(Sequential_t,Index,const Scalar&,const Scalar&,Index), CwiseNullaryOp
   */
 template<typename Derived>
@@ -377,7 +387,10 @@ PlainObjectBase<Derived>::setConstant(Index rows, Index cols, const Scalar& val)
   * Example: \include DenseBase_setLinSpaced.cpp
   * Output: \verbinclude DenseBase_setLinSpaced.out
   *
-  * \sa CwiseNullaryOp
+  * For integer scalar types, do not miss the explanations on the definition
+  * of \link LinSpaced(Index,const Scalar&,const Scalar&) even spacing \endlink.
+  *
+  * \sa LinSpaced(Index,const Scalar&,const Scalar&), CwiseNullaryOp
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(Index newSize, const Scalar& low, const Scalar& high)
@@ -389,12 +402,15 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(Index newSize, con
 /**
   * \brief Sets a linearly spaced vector.
   *
-  * The function fills *this with equally spaced values in the closed interval [low,high].
+  * The function fills \c *this with equally spaced values in the closed interval [low,high].
   * When size is set to 1, a vector of length 1 containing 'high' is returned.
   *
   * \only_for_vectors
   *
-  * \sa setLinSpaced(Index, const Scalar&, const Scalar&), CwiseNullaryOp
+  * For integer scalar types, do not miss the explanations on the definition
+  * of \link LinSpaced(Index,const Scalar&,const Scalar&) even spacing \endlink.
+  *
+  * \sa LinSpaced(Index,const Scalar&,const Scalar&), setLinSpaced(Index, const Scalar&, const Scalar&), CwiseNullaryOp
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(const Scalar& low, const Scalar& high)
diff --git a/Eigen/src/Core/functors/NullaryFunctors.h b/Eigen/src/Core/functors/NullaryFunctors.h
index a2154d3b5..0c00f4661 100644
--- a/Eigen/src/Core/functors/NullaryFunctors.h
+++ b/Eigen/src/Core/functors/NullaryFunctors.h
@@ -99,25 +99,24 @@ template <typename Scalar, typename Packet>
 struct linspaced_op_impl<Scalar,Packet,/*RandomAccess*/true,/*IsInteger*/true>
 {
   linspaced_op_impl(const Scalar& low, const Scalar& high, Index num_steps) :
-    m_low(low), m_length(high-low), m_divisor(convert_index<Scalar>(num_steps==1?1:num_steps-1)), m_interPacket(plset<Packet>(0))
+    m_low(low),
+    m_multiplier((high-low)/convert_index<Scalar>(num_steps<=1 ? 1 : num_steps-1)),
+    m_divisor(convert_index<Scalar>(num_steps+high-low)/(high-low+1)),
+    m_use_divisor((high-low+1)<num_steps)
   {}
 
   template<typename IndexType>
   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
-  const Scalar operator() (IndexType i) const {
-    return m_low + (m_length*Scalar(i))/m_divisor;
+  const Scalar operator() (IndexType i) const
+  {
+    if(m_use_divisor) return m_low + convert_index<Scalar>(i)/m_divisor;
+    else              return m_low + convert_index<Scalar>(i)*m_multiplier;
   }
 
-  template<typename IndexType>
-  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
-  const Packet packetOp(IndexType i) const {
-    return internal::padd(pset1<Packet>(m_low), pdiv(pmul(pset1<Packet>(m_length), padd(pset1<Packet>(Scalar(i)),m_interPacket)),
-                                                     pset1<Packet>(m_divisor))); }
-
   const Scalar m_low;
-  const Scalar m_length;
-  const Scalar  m_divisor;
-  const Packet m_interPacket;
+  const Scalar m_multiplier;
+  const Scalar m_divisor;
+  const bool m_use_divisor;
 };
 
 // ----- Linspace functor ----------------------------------------------------------------
@@ -131,7 +130,7 @@ template <typename Scalar, typename PacketType, bool RandomAccess> struct functo
   enum
   {
     Cost = 1,
-    PacketAccess =   packet_traits<Scalar>::HasSetLinear
+    PacketAccess =   (!NumTraits<Scalar>::IsInteger) && packet_traits<Scalar>::HasSetLinear
                   && ((!NumTraits<Scalar>::IsInteger) || packet_traits<Scalar>::HasDiv),
     IsRepeatable = true
   };
diff --git a/doc/snippets/DenseBase_LinSpacedInt.cpp b/doc/snippets/DenseBase_LinSpacedInt.cpp
new file mode 100644
index 000000000..0d7ae068e
--- /dev/null
+++ b/doc/snippets/DenseBase_LinSpacedInt.cpp
@@ -0,0 +1,8 @@
+cout << "Even spacing inputs:" << endl;
+cout << VectorXi::LinSpaced(8,1,4).transpose() << endl;
+cout << VectorXi::LinSpaced(8,1,8).transpose() << endl;
+cout << VectorXi::LinSpaced(8,1,15).transpose() << endl;
+cout << "Uneven spacing inputs:" << endl;
+cout << VectorXi::LinSpaced(8,1,7).transpose() << endl;
+cout << VectorXi::LinSpaced(8,1,9).transpose() << endl;
+cout << VectorXi::LinSpaced(8,1,16).transpose() << endl;
diff --git a/test/nullary.cpp b/test/nullary.cpp
index 35f24de47..162e84210 100644
--- a/test/nullary.cpp
+++ b/test/nullary.cpp
@@ -57,23 +57,31 @@ void testVectorType(const VectorType& base)
     VERIFY_IS_APPROX(m,n);
   }
 
-  VectorType n(size);
-  for (int i=0; i<size; ++i)
-    n(i) = size==1 ? low : (low + ((high-low)*Scalar(i))/(size-1));
-  VERIFY_IS_APPROX(m,n);
-
-  // random access version
-  m = VectorType::LinSpaced(size,low,high);
-  VERIFY_IS_APPROX(m,n);
+  if((!NumTraits<Scalar>::IsInteger) || ((high-low)>=size && (Index(high-low)%(size-1))==0) || (Index(high-low+1)<size && (size%Index(high-low+1))==0))
+  {
+    VectorType n(size);
+    if((!NumTraits<Scalar>::IsInteger) || (high-low>=size))
+      for (int i=0; i<size; ++i)
+        n(i) = size==1 ? low : (low + ((high-low)*Scalar(i))/(size-1));
+    else
+      for (int i=0; i<size; ++i)
+        n(i) = size==1 ? low : low + Scalar((double(high-low+1)*double(i))/double(size));
+    VERIFY_IS_APPROX(m,n);
 
-  VERIFY( internal::isApprox(m(m.size()-1),high) );
-  VERIFY( size==1 || internal::isApprox(m(0),low) );
+    // random access version
+    m = VectorType::LinSpaced(size,low,high);
+    VERIFY_IS_APPROX(m,n);
+    VERIFY( internal::isApprox(m(m.size()-1),high) );
+    VERIFY( size==1 || internal::isApprox(m(0),low) );
 
-  // sequential access version
-  m = VectorType::LinSpaced(Sequential,size,low,high);
-  VERIFY_IS_APPROX(m,n);
+    // sequential access version
+    m = VectorType::LinSpaced(Sequential,size,low,high);
+    VERIFY_IS_APPROX(m,n);
+    VERIFY( internal::isApprox(m(m.size()-1),high) );
+  }
 
-  VERIFY( internal::isApprox(m(m.size()-1),high) );
+  Scalar tol_factor = (high>=0) ? (1+NumTraits<Scalar>::dummy_precision()) : (1-NumTraits<Scalar>::dummy_precision());
+  VERIFY( m(m.size()-1) <= high*tol_factor );
   VERIFY( size==1 || internal::isApprox(m(0),low) );
 
   // check whether everything works with row and col major vectors
@@ -96,7 +104,7 @@ void testVectorType(const VectorType& base)
   VERIFY_IS_APPROX( ScalarMatrix::LinSpaced(1,low,high), ScalarMatrix::Constant(high) );
 
   // regression test for bug 526 (linear vectorized transversal)
-  if (size > 1) {
+  if (size > 1 && (!NumTraits<Scalar>::IsInteger)) {
     m.tail(size-1).setLinSpaced(low, high);
     VERIFY_IS_APPROX(m(size-1), high);
   }