pl/math: Fix fp exceptions in Neon sinhf and sinh

Both routines previously relied on the vector expm1(f) routine exposed
by the library, which depended on WANT_SIMD_EXCEPT for its fenv
behaviour, however both routines were expected to always trigger fp
exceptions correctly. To remedy this, both routines now use an inlined
helper for expm1 (reused from vector tanhf in the case of sinhf), and
special-case small input as well as large when WANT_SIMD_EXCEPT is
enabled.

1 files changed, 27 insertions, 15 deletions

diff --git a/pl/math/v_sinhf_2u3.c b/pl/math/v_sinhf_2u3.c
index ce2fe0e..a54c178 100644
--- a/pl/math/v_sinhf_2u3.c
+++ b/pl/math/v_sinhf_2u3.c
@@ -5,17 +5,25 @@
  */
 
 #include "v_math.h"
-#include "mathlib.h"
 #include "pl_sig.h"
 #include "pl_test.h"
 
 #if V_SUPPORTED
 
+#include "v_expm1f_inline.h"
+
 #define AbsMask 0x7fffffff
 #define Half 0x3f000000
-#define Expm1OFlowLimit                                                        \
-  0x42b17218 /* 0x1.62e43p+6, 2^7*ln2, minimum value for which expm1f          \
-		overflows.  */
+#define BigBound                                                               \
+  0x42b0c0a7 /* 0x1.61814ep+6, above which expm1f helper overflows.  */
+#define TinyBound                                                              \
+  0x2fb504f4 /* 0x1.6a09e8p-32, below which expm1f underflows.  */
+
+static NOINLINE VPCS_ATTR v_f32_t
+special_case (v_f32_t x)
+{
+  return v_call_f32 (sinhf, x, x, v_u32 (-1));
+}
 
 /* Approximation for vector single-precision sinh(x) using expm1.
    sinh(x) = (exp(x) - exp(-x)) / 2.
@@ -29,28 +37,32 @@ VPCS_ATTR v_f32_t V_NAME (sinhf) (v_f32_t x)
   v_u32_t sign = ix & ~AbsMask;
   v_f32_t halfsign = v_as_f32_u32 (sign | Half);
 
-  v_u32_t special = v_cond_u32 (iax >= Expm1OFlowLimit);
+#if WANT_SIMD_EXCEPT
+  v_u32_t special = v_cond_u32 ((iax - TinyBound) >= (BigBound - TinyBound));
+#else
+  v_u32_t special = v_cond_u32 (iax >= BigBound);
+#endif
+
   /* Fall back to the scalar variant for all lanes if any of them should trigger
      an exception.  */
   if (unlikely (v_any_u32 (special)))
-    return v_call_f32 (sinhf, x, x, v_u32 (-1));
+    return special_case (x);
 
   /* Up to the point that expm1f overflows, we can use it to calculate sinhf
      using a slight rearrangement of the definition of asinh. This allows us to
      retain acceptable accuracy for very small inputs.  */
-  v_f32_t t = V_NAME (expm1f) (ax);
+  v_f32_t t = expm1f_inline (ax);
   return (t + t / (t + 1)) * halfsign;
 }
 VPCS_ALIAS
 
 PL_SIG (V, F, 1, sinh, -10.0, 10.0)
 PL_TEST_ULP (V_NAME (sinhf), 1.76)
-/* TODO: reinstate PL_TEST_EXPECT_FENV here once fp exceptions are triggered
-   correctly.  */
-PL_TEST_INTERVAL (V_NAME (sinhf), 0, 0x1.62e43p+6, 100000)
-PL_TEST_INTERVAL (V_NAME (sinhf), -0, -0x1.62e43p+6, 100000)
-PL_TEST_INTERVAL (V_NAME (sinhf), 0x1.62e43p+6, 0x1.65a9fap+6, 100)
-PL_TEST_INTERVAL (V_NAME (sinhf), -0x1.62e43p+6, -0x1.65a9fap+6, 100)
-PL_TEST_INTERVAL (V_NAME (sinhf), 0x1.65a9fap+6, inf, 100)
-PL_TEST_INTERVAL (V_NAME (sinhf), -0x1.65a9fap+6, -inf, 100)
+PL_TEST_EXPECT_FENV (V_NAME (sinhf), WANT_SIMD_EXCEPT)
+PL_TEST_INTERVAL (V_NAME (sinhf), 0, TinyBound, 1000)
+PL_TEST_INTERVAL (V_NAME (sinhf), -0, -TinyBound, 1000)
+PL_TEST_INTERVAL (V_NAME (sinhf), TinyBound, BigBound, 100000)
+PL_TEST_INTERVAL (V_NAME (sinhf), -TinyBound, -BigBound, 100000)
+PL_TEST_INTERVAL (V_NAME (sinhf), BigBound, inf, 1000)
+PL_TEST_INTERVAL (V_NAME (sinhf), -BigBound, -inf, 1000)
 #endif