From 38d704def8b6799b14c319d6a67c671374daccc3 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Sun, 23 Dec 2018 16:13:24 +0100
Subject: [PATCH] Make sure that psin/pcos return number in [-1,1] for large
 inputs (though sin/cos on large entries is quite useless because it's
 inaccurate)

---
 .../arch/Default/GenericPacketMathFunctions.h    |  6 ++++++
 test/packetmath.cpp                              | 16 ++++++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
index 83fed95de..80bcc077d 100644
--- a/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
+++ b/Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
@@ -289,9 +289,15 @@ Packet psincos_float(const Packet& _x)
   const Packet cst_coscof_p1        = pset1<Packet>(-1.388731625493765E-003f);
   const Packet cst_coscof_p2        = pset1<Packet>( 4.166664568298827E-002f);
   const Packet cst_cephes_FOPI      = pset1<Packet>( 1.27323954473516f); // 4 / M_PI
+  const Packet cst_sincos_max_arg   = pset1<Packet>( 13176795.0f);  // Approx. (2**24) / (4/Pi).
 
   Packet x = pabs(_x);
 
+  // Prevent sin/cos from generating values larger than 1.0 in magnitude
+  // for very large arguments by setting x to 0.0.
+  Packet small_or_nan_mask = pcmp_lt_or_nan(x, cst_sincos_max_arg);
+  x = pand(x, small_or_nan_mask);
+
   // Scale x by 4/Pi to find x's octant.
   Packet y = pmul(x, cst_cephes_FOPI);
 
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
index 7e46b01de..bf0312a73 100644
--- a/test/packetmath.cpp
+++ b/test/packetmath.cpp
@@ -564,6 +564,22 @@ template<typename Scalar,typename Packet> void packetmath_real()
       VERIFY((numext::isnan)(data2[0]));
       VERIFY((numext::isnan)(data2[1]));
     }
+    if(PacketTraits::HasCos)
+    {
+      packet_helper<PacketTraits::HasCos,Packet> h;
+      for(Scalar k = 1; k<Scalar(1000)/std::numeric_limits<Scalar>::epsilon(); k*=2) {
+        data1[0] = k*Scalar(EIGEN_PI) * internal::random<Scalar>(0.8,1.2);
+        data1[1] = (k+1)*Scalar(EIGEN_PI) * internal::random<Scalar>(0.8,1.2);
+        h.store(data2, internal::pcos(h.load(data1)));
+        VERIFY(data2[0]<=Scalar(1.) && data2[0]>=Scalar(-1.));
+        VERIFY(data2[1]<=Scalar(1.) && data2[1]>=Scalar(-1.));
+        data1[0] = (2*k+1)*Scalar(EIGEN_PI)/2 * internal::random<Scalar>(0.8,1.2);
+        data1[1] = (2*k+3)*Scalar(EIGEN_PI)/2 * internal::random<Scalar>(0.8,1.2);
+        h.store(data2, internal::psin(h.load(data1)));
+        VERIFY(data2[0]<=Scalar(1.) && data2[0]>=Scalar(-1.));
+        VERIFY(data2[1]<=Scalar(1.) && data2[1]>=Scalar(-1.));
+      }
+    }
   }
 }