Attachment #960 for bug #1777

This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen

View | Details | Raw Unified | Return to bug 1777
Collapse All | Expand All




  y = pmadd(y, r, cst_cephes_exp_p5);
  y = pmadd(y, r2, r);
  y = padd(y, cst_1);

  // Return 2^m * exp(r).
  return pmax(pldexp(y,m), _x);
}

// make it the default path for scalar float
template<>
float pexp(const float& a) { return pexp_float(a); }

template <typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
EIGEN_UNUSED
Packet pexp_double(const Packet _x)
{
  Packet x = _x;

  const Packet cst_1 = pset1<Packet>(1.0);

  x = pdiv(px, psub(qx, px));
  x = pmadd(cst_2, x, cst_1);

  // Construct the result 2^n * exp(g) = e * x. The max is used to catch
  // non-finite values in the input.
  return pmax(pldexp(x,fx), _x);
}

// make it the default path for scalar double
template<>
double pexp(const double& a) { return pexp_double(a); }

// The following code is inspired by the following stack-overflow answer:
//   https://stackoverflow.com/questions/30463616/payne-hanek-algorithm-implementation-in-c/30465751#30465751
// It has been largely optimized:
//  - By-pass calls to frexp.
//  - Aligned loads of required 96 bits of 2/pi. This is accomplished by
//    (1) balancing the mantissa and exponent to the required bits of 2/pi are
//    aligned on 8-bits, and (2) replicating the storage of the bits of 2/pi.
//  - Avoid a branch in rounding and extraction of the remaining fractional part.




/** \internal
  * \brief Template functor to compute the logistic function of a scalar
  * \sa class CwiseUnaryOp, ArrayBase::logistic()
  */
template <typename T>
struct scalar_logistic_op {
  EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(const T& x) const {
    return packetOp(x);
    return one / (one + numext::exp(-x));
  }

  template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
  Packet packetOp(const Packet& x) const {
    const Packet one = pset1<Packet>(T(1));
    return pdiv(one, padd(one, pexp(pnegate(x))));
  }
};

  *  logistic is interpolated because it was easier to make the fit converge.
  *
  */

template <>
struct scalar_logistic_op<float> {
  EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float operator()(const float& x) const {
    return packetOp(x);
    else if (x > 18.0f) return 1.0f;
    else return 1.0f / (1.0f + numext::exp(-x));
  }

  template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
  Packet packetOp(const Packet& _x) const {
    // Clamp the inputs to the range [-18, 18] since anything outside
    // this range is 0.0f or 1.0f in single-precision.
    const Packet x = pmax(pmin(_x, pset1<Packet>(18.0)), pset1<Packet>(-18.0));





  if (PacketTraits::HasTanh) {
    // NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details.
    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
    packet_helper<internal::packet_traits<Scalar>::HasTanh,Packet> h;
    h.store(data2, internal::ptanh(h.load(data1)));
    VERIFY((numext::isnan)(data2[0]));
  }

  {
    internal::scalar_logistic_op<Scalar> logistic;
    for (int i=0; i<size; ++i)
    {
      data1[i] = internal::random<Scalar>(-20,20);
    }
    internal::pstore(data2, logistic.packetOp(internal::pload<Packet>(data1)));
    for (int i=0; i<PacketSize; ++i) {
      VERIFY_IS_APPROX(data2[i],logistic(data1[i]));
      #ifdef EIGEN_VECTORIZE // don't check for exactness when using the i387 FPU
      VERIFY_IS_EQUAL(data2[i],logistic(data1[i]));
      #endif
    }
  }

#if EIGEN_HAS_C99_MATH
  {
    data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
    packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h;
    h.store(data2, internal::plgamma(h.load(data1)));
    VERIFY((numext::isnan)(data2[0]));
  }
  if (internal::packet_traits<Scalar>::HasErf) {

    runall<Scalar,PacketType>::run();
  }
};

EIGEN_DECLARE_TEST(packetmath)
{
  g_first_pass = true;
  for(int i = 0; i < g_repeat; i++) {

    CALL_SUBTEST_1( runner<float>::run() );
    CALL_SUBTEST_2( runner<double>::run() );
    CALL_SUBTEST_3( runner<int>::run() );
    CALL_SUBTEST_4( runner<std::complex<float> >::run() );
    CALL_SUBTEST_5( runner<std::complex<double> >::run() );
    CALL_SUBTEST_6(( packetmath<half,internal::packet_traits<half>::type>() ));
    g_first_pass = false;
  }

Return to bug 1777

Lines 228-243 Packet pexp_float(const Packet _x) Link Here

(-)a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h (+8 lines)
228	y = pmadd(y, r, cst_cephes_exp_p5);	228	y = pmadd(y, r, cst_cephes_exp_p5);
229	y = pmadd(y, r2, r);	229	y = pmadd(y, r2, r);
230	y = padd(y, cst_1);	230	y = padd(y, cst_1);
231		231
232	// Return 2^m * exp(r).	232	// Return 2^m * exp(r).
233	return pmax(pldexp(y,m), _x);	233	return pmax(pldexp(y,m), _x);
234	}	234	}
235		235
		236	// make it the default path for scalar float
		237	template<>
		238	float pexp(const float& a) { return pexp_float(a); }
		239
236	template <typename Packet>	240	template <typename Packet>
237	EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS	241	EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
238	EIGEN_UNUSED	242	EIGEN_UNUSED
239	Packet pexp_double(const Packet _x)	243	Packet pexp_double(const Packet _x)
240	{	244	{
241	Packet x = _x;	245	Packet x = _x;
242		246
243	const Packet cst_1 = pset1<Packet>(1.0);	247	const Packet cst_1 = pset1<Packet>(1.0);
Lines 296-311 Packet pexp_double(const Packet _x) Link Here
296	x = pdiv(px, psub(qx, px));	300	x = pdiv(px, psub(qx, px));
297	x = pmadd(cst_2, x, cst_1);	301	x = pmadd(cst_2, x, cst_1);
298		302
299	// Construct the result 2^n * exp(g) = e * x. The max is used to catch	303	// Construct the result 2^n * exp(g) = e * x. The max is used to catch
300	// non-finite values in the input.	304	// non-finite values in the input.
301	return pmax(pldexp(x,fx), _x);	305	return pmax(pldexp(x,fx), _x);
302	}	306	}
303		307
		308	// make it the default path for scalar double
		309	template<>
		310	double pexp(const double& a) { return pexp_double(a); }
		311
304	// The following code is inspired by the following stack-overflow answer:	312	// The following code is inspired by the following stack-overflow answer:
305	// https://stackoverflow.com/questions/30463616/payne-hanek-algorithm-implementation-in-c/30465751#30465751	313	// https://stackoverflow.com/questions/30463616/payne-hanek-algorithm-implementation-in-c/30465751#30465751
306	// It has been largely optimized:	314	// It has been largely optimized:
307	// - By-pass calls to frexp.	315	// - By-pass calls to frexp.
308	// - Aligned loads of required 96 bits of 2/pi. This is accomplished by	316	// - Aligned loads of required 96 bits of 2/pi. This is accomplished by
309	// (1) balancing the mantissa and exponent to the required bits of 2/pi are	317	// (1) balancing the mantissa and exponent to the required bits of 2/pi are
310	// aligned on 8-bits, and (2) replicating the storage of the bits of 2/pi.	318	// aligned on 8-bits, and (2) replicating the storage of the bits of 2/pi.
311	// - Avoid a branch in rounding and extraction of the remaining fractional part.	319	// - Avoid a branch in rounding and extraction of the remaining fractional part.

Lines 566-581 template<typename Scalar,typename Packet Link Here

(-)a/test/packetmath.cpp (-1 / +15 lines)
566	if (PacketTraits::HasTanh) {	566	if (PacketTraits::HasTanh) {
567	// NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details.	567	// NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details.
568	data1[0] = std::numeric_limits<Scalar>::quiet_NaN();	568	data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
569	packet_helper<internal::packet_traits<Scalar>::HasTanh,Packet> h;	569	packet_helper<internal::packet_traits<Scalar>::HasTanh,Packet> h;
570	h.store(data2, internal::ptanh(h.load(data1)));	570	h.store(data2, internal::ptanh(h.load(data1)));
571	VERIFY((numext::isnan)(data2[0]));	571	VERIFY((numext::isnan)(data2[0]));
572	}	572	}
573		573
		574	{
		575	internal::scalar_logistic_op<Scalar> logistic;
		576	for (int i=0; i<size; ++i)
		577	{
		578	data1[i] = internal::random<Scalar>(-20,20);
		579	}
		580	internal::pstore(data2, logistic.packetOp(internal::pload<Packet>(data1)));
		581	for (int i=0; i<PacketSize; ++i) {
		582	VERIFY_IS_APPROX(data2[i],logistic(data1[i]));
		583	#ifdef EIGEN_VECTORIZE // don't check for exactness when using the i387 FPU
		584	VERIFY_IS_EQUAL(data2[i],logistic(data1[i]));
		585	#endif
		586	}
		587	}
		588
574	#if EIGEN_HAS_C99_MATH	589	#if EIGEN_HAS_C99_MATH
575	{	590	{
576	data1[0] = std::numeric_limits<Scalar>::quiet_NaN();	591	data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
577	packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h;	592	packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h;
578	h.store(data2, internal::plgamma(h.load(data1)));	593	h.store(data2, internal::plgamma(h.load(data1)));
579	VERIFY((numext::isnan)(data2[0]));	594	VERIFY((numext::isnan)(data2[0]));
580	}	595	}
581	if (internal::packet_traits<Scalar>::HasErf) {	596	if (internal::packet_traits<Scalar>::HasErf) {
Lines 960-976 struct runner<Scalar,PacketType,false,fa Link Here
960	runall<Scalar,PacketType>::run();	975	runall<Scalar,PacketType>::run();
961	}	976	}
962	};	977	};
963		978
964	EIGEN_DECLARE_TEST(packetmath)	979	EIGEN_DECLARE_TEST(packetmath)
965	{	980	{
966	g_first_pass = true;	981	g_first_pass = true;
967	for(int i = 0; i < g_repeat; i++) {	982	for(int i = 0; i < g_repeat; i++) {
968
969	CALL_SUBTEST_1( runner<float>::run() );	983	CALL_SUBTEST_1( runner<float>::run() );
970	CALL_SUBTEST_2( runner<double>::run() );	984	CALL_SUBTEST_2( runner<double>::run() );
971	CALL_SUBTEST_3( runner<int>::run() );	985	CALL_SUBTEST_3( runner<int>::run() );
972	CALL_SUBTEST_4( runner<std::complex<float> >::run() );	986	CALL_SUBTEST_4( runner<std::complex<float> >::run() );
973	CALL_SUBTEST_5( runner<std::complex<double> >::run() );	987	CALL_SUBTEST_5( runner<std::complex<double> >::run() );
974	CALL_SUBTEST_6(( packetmath<half,internal::packet_traits<half>::type>() ));	988	CALL_SUBTEST_6(( packetmath<half,internal::packet_traits<half>::type>() ));
975	g_first_pass = false;	989	g_first_pass = false;
976	}	990	}

Lines 889-906 struct functor_traits<scalar_sign_op<Sca Link Here

(-)a/Eigen/src/Core/functors/UnaryFunctors.h (-5 / +2 lines)
889	/** \internal	889	/** \internal
890	* \brief Template functor to compute the logistic function of a scalar	890	* \brief Template functor to compute the logistic function of a scalar
891	* \sa class CwiseUnaryOp, ArrayBase::logistic()	891	* \sa class CwiseUnaryOp, ArrayBase::logistic()
892	*/	892	*/
893	template <typename T>	893	template <typename T>
894	struct scalar_logistic_op {	894	struct scalar_logistic_op {
895	EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)	895	EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)
896	EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(const T& x) const {	896	EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE T operator()(const T& x) const {
897	const T one = T(1);	897	return packetOp(x);
898	return one / (one + numext::exp(-x));
899	}	898	}
900		899
901	template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE	900	template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
902	Packet packetOp(const Packet& x) const {	901	Packet packetOp(const Packet& x) const {
903	const Packet one = pset1<Packet>(T(1));	902	const Packet one = pset1<Packet>(T(1));
904	return pdiv(one, padd(one, pexp(pnegate(x))));	903	return pdiv(one, padd(one, pexp(pnegate(x))));
905	}	904	}
906	};	905	};
Lines 914-932 struct scalar_logistic_op { Link Here
914	* logistic is interpolated because it was easier to make the fit converge.	913	* logistic is interpolated because it was easier to make the fit converge.
915	*	914	*
916	*/	915	*/
917		916
918	template <>	917	template <>
919	struct scalar_logistic_op<float> {	918	struct scalar_logistic_op<float> {
920	EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)	919	EIGEN_EMPTY_STRUCT_CTOR(scalar_logistic_op)
921	EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float operator()(const float& x) const {	920	EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE float operator()(const float& x) const {
922	if (x < -18.0f) return 0.0f;	921	return packetOp(x);
923	else if (x > 18.0f) return 1.0f;
924	else return 1.0f / (1.0f + numext::exp(-x));
925	}	922	}
926		923
927	template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE	924	template <typename Packet> EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
928	Packet packetOp(const Packet& _x) const {	925	Packet packetOp(const Packet& _x) const {
929	// Clamp the inputs to the range [-18, 18] since anything outside	926	// Clamp the inputs to the range [-18, 18] since anything outside
930	// this range is 0.0f or 1.0f in single-precision.	927	// this range is 0.0f or 1.0f in single-precision.
931	const Packet x = pmax(pmin(_x, pset1<Packet>(18.0)), pset1<Packet>(-18.0));	928	const Packet x = pmax(pmin(_x, pset1<Packet>(18.0)), pset1<Packet>(-18.0));
932		929