mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-18 14:34:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
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.
This commit is contained in:
Gael Guennebaud 2016-10-24 15:50:27 +02:00
parent e8e56c7642
commit 53c77061f0
4 changed files with 62 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Eigen/src/Core/CwiseNullaryOp.h
View File

@ -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)

25
Eigen/src/Core/functors/NullaryFunctors.h
View File

@ -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
};

8
doc/snippets/DenseBase_LinSpacedInt.cpp Normal file
View File

@ -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;

38
test/nullary.cpp
View File

@ -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);
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);
// random access version
m = VectorType::LinSpaced(size,low,high);
VERIFY_IS_APPROX(m,n);
// 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) );
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);
VERIFY( internal::isApprox(m(m.size()-1),high) );
}
// sequential access version
m = VectorType::LinSpaced(Sequential,size,low,high);
VERIFY_IS_APPROX(m,n);
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);
}