diff --git a/Eigen/src/Core/functors/NullaryFunctors.h b/Eigen/src/Core/functors/NullaryFunctors.h index 128013aba..0000ea1f1 100644 --- a/Eigen/src/Core/functors/NullaryFunctors.h +++ b/Eigen/src/Core/functors/NullaryFunctors.h @@ -43,12 +43,17 @@ template struct linspaced_op_impl { linspaced_op_impl(const Scalar& low, const Scalar& high, Index num_steps) : - m_low(low), m_high(high), m_size1(num_steps==1 ? 1 : num_steps-1), m_step(num_steps==1 ? Scalar() : (high-low)/Scalar(num_steps-1)), m_interPacket(plset(0)) + m_low(low), m_high(high), m_size1(num_steps==1 ? 1 : num_steps-1), m_step(num_steps==1 ? Scalar() : (high-low)/Scalar(num_steps-1)), + m_interPacket(plset(0)), + m_flip(std::abs(high) EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (IndexType i) const { - return (i==m_size1)? m_high : (m_low + i*m_step); + if(m_flip) + return (i==0)? m_low : (m_high - (m_size1-i)*m_step); + else + return (i==m_size1)? m_high : (m_low + i*m_step); } template @@ -56,11 +61,22 @@ struct linspaced_op_impl { // Principle: // [low, ..., low] + ( [step, ..., step] * ( [i, ..., i] + [0, ..., size] ) ) - Packet pi = padd(pset1(Scalar(i)),m_interPacket); - Packet res = padd(pset1(m_low), pmul(pset1(m_step), pi)); - if(i==m_size1-unpacket_traits::size+1) - res = pinsertlast(res, m_high); - return res; + if(m_flip) + { + Packet pi = padd(pset1(Scalar(i-m_size1)),m_interPacket); + Packet res = padd(pset1(m_high), pmul(pset1(m_step), pi)); + if(i==0) + res = pinsertfirst(res, m_low); + return res; + } + else + { + Packet pi = padd(pset1(Scalar(i)),m_interPacket); + Packet res = padd(pset1(m_low), pmul(pset1(m_step), pi)); + if(i==m_size1-unpacket_traits::size+1) + res = pinsertlast(res, m_high); + return res; + } } const Scalar m_low; @@ -68,6 +84,7 @@ struct linspaced_op_impl const Index m_size1; const Scalar m_step; const Packet m_interPacket; + const bool m_flip; }; template diff --git a/test/nullary.cpp b/test/nullary.cpp index 6d16bd4d2..351d26e74 100644 --- a/test/nullary.cpp +++ b/test/nullary.cpp @@ -33,10 +33,38 @@ bool equalsIdentity(const MatrixType& A) } +template +void check_extremity_accuracy(const VectorType &v, const typename VectorType::Scalar &low, const typename VectorType::Scalar &high) +{ + typedef typename VectorType::Scalar Scalar; + typedef typename VectorType::RealScalar RealScalar; + + RealScalar prec = internal::is_same::value ? NumTraits::dummy_precision()*10 : NumTraits::dummy_precision()/10; + Index size = v.size(); + + if(size<20) + return; + + for (int i=0; isize-6) + { + Scalar ref = (low*RealScalar(size-i-1))/RealScalar(size-1) + (high*RealScalar(i))/RealScalar(size-1); + if(std::abs(ref)>1) + { + if(!internal::isApprox(v(i), ref, prec)) + std::cout << v(i) << " != " << ref << " ; relative error: " << std::abs((v(i)-ref)/ref) << " ; required precision: " << prec << " ; range: " << low << "," << high << " ; i: " << i << "\n"; + VERIFY(internal::isApprox(v(i), (low*RealScalar(size-i-1))/RealScalar(size-1) + (high*RealScalar(i))/RealScalar(size-1), prec)); + } + } + } +} + template void testVectorType(const VectorType& base) { typedef typename VectorType::Scalar Scalar; + typedef typename VectorType::RealScalar RealScalar; const Index size = base.size(); @@ -47,6 +75,9 @@ void testVectorType(const VectorType& base) // check low==high if(internal::random(0.f,1.f)<0.05f) low = high; + // check abs(low) >> abs(high) + else if(size>2 && std::numeric_limits::max_exponent10>0 && internal::random(0.f,1.f)<0.1f) + low = -internal::random(1,2) * RealScalar(std::pow(RealScalar(10),std::numeric_limits::max_exponent10/2)); const Scalar step = ((size == 1) ? 1 : (high-low)/(size-1)); @@ -60,6 +91,8 @@ void testVectorType(const VectorType& base) for (int i=0; i::IsInteger) || ((high-low)>=size && (Index(high-low)%(size-1))==0) || (Index(high-low+1)::IsInteger) + CALL_SUBTEST( check_extremity_accuracy(m, low, high) ); } VERIFY( m(m.size()-1) <= high ); @@ -154,10 +189,10 @@ void test_nullary() CALL_SUBTEST_3( testMatrixType(MatrixXf(internal::random(1,300),internal::random(1,300))) ); for(int i = 0; i < g_repeat*10; i++) { - CALL_SUBTEST_4( testVectorType(VectorXd(internal::random(1,300))) ); + CALL_SUBTEST_4( testVectorType(VectorXd(internal::random(1,30000))) ); CALL_SUBTEST_5( testVectorType(Vector4d()) ); // regression test for bug 232 CALL_SUBTEST_6( testVectorType(Vector3d()) ); - CALL_SUBTEST_7( testVectorType(VectorXf(internal::random(1,300))) ); + CALL_SUBTEST_7( testVectorType(VectorXf(internal::random(1,30000))) ); CALL_SUBTEST_8( testVectorType(Vector3f()) ); CALL_SUBTEST_8( testVectorType(Vector4f()) ); CALL_SUBTEST_8( testVectorType(Matrix()) );