From 0a9270a27f48bea87c5bd3f0f9c759da66fb45a3 Mon Sep 17 00:00:00 2001
From: Joe Ramsay <Joe.Ramsay@arm.com>
Date: Thu, 22 Dec 2022 16:20:22 +0000
Subject: 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.
---
 pl/math/v_expm1f_inline.h | 49 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 49 insertions(+)
 create mode 100644 pl/math/v_expm1f_inline.h

(limited to 'pl/math/v_expm1f_inline.h')

diff --git a/pl/math/v_expm1f_inline.h b/pl/math/v_expm1f_inline.h
new file mode 100644
index 0000000..ef9e934
--- /dev/null
+++ b/pl/math/v_expm1f_inline.h
@@ -0,0 +1,49 @@
+/*
+ * Helper for single-precision routines which calculate exp(x) - 1 and do not
+ * need special-case handling
+ *
+ * Copyright (c) 2022, Arm Limited.
+ * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+ */
+
+#ifndef PL_MATH_V_EXPM1F_INLINE_H
+#define PL_MATH_V_EXPM1F_INLINE_H
+
+#include "v_math.h"
+#include "math_config.h"
+#include "estrinf.h"
+
+#define One 0x3f800000
+#define Shift v_f32 (0x1.8p23f)
+#define InvLn2 v_f32 (0x1.715476p+0f)
+#define MLn2hi v_f32 (-0x1.62e4p-1f)
+#define MLn2lo v_f32 (-0x1.7f7d1cp-20f)
+
+#define C(i) v_f32 (__expm1f_poly[i])
+
+static inline v_f32_t
+expm1f_inline (v_f32_t x)
+{
+  /* Helper routine for calculating exp(x) - 1.
+     Copied from v_expm1f_1u6.c, with all special-case handling removed - the
+     calling routine should handle special values if required.  */
+
+  /* Reduce argument: f in [-ln2/2, ln2/2], i is exact.  */
+  v_f32_t j = v_fma_f32 (InvLn2, x, Shift) - Shift;
+  v_s32_t i = v_to_s32_f32 (j);
+  v_f32_t f = v_fma_f32 (j, MLn2hi, x);
+  f = v_fma_f32 (j, MLn2lo, f);
+
+  /* Approximate expm1(f) with polynomial P, expm1(f) ~= f + f^2 * P(f).
+     Uses Estrin scheme, where the main __v_expm1f routine uses Horner.  */
+  v_f32_t f2 = f * f;
+  v_f32_t p = ESTRIN_4 (f, f2, f2 * f2, C);
+  p = v_fma_f32 (f2, p, f);
+
+  /* t = 2^i.  */
+  v_f32_t t = v_as_f32_u32 (v_as_u32_s32 (i << 23) + One);
+  /* expm1(x) ~= p * t + (t - 1).  */
+  return v_fma_f32 (p, t, t - 1);
+}
+
+#endif // PL_MATH_V_EXPM1F_INLINE_H
-- 
cgit v1.2.3


From f0f80b8a19b2593491847ed87456694d789f6f80 Mon Sep 17 00:00:00 2001
From: Joe Ramsay <Joe.Ramsay@arm.com>
Date: Fri, 6 Jan 2023 09:10:57 +0000
Subject: pl/math: Update copyright years

All files in pl/math updated to 2023.
---
 pl/math/v_expm1f_inline.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'pl/math/v_expm1f_inline.h')

diff --git a/pl/math/v_expm1f_inline.h b/pl/math/v_expm1f_inline.h
index ef9e934..c261941 100644
--- a/pl/math/v_expm1f_inline.h
+++ b/pl/math/v_expm1f_inline.h
@@ -2,7 +2,7 @@
  * Helper for single-precision routines which calculate exp(x) - 1 and do not
  * need special-case handling
  *
- * Copyright (c) 2022, Arm Limited.
+ * Copyright (c) 2022-2023, Arm Limited.
  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
  */
 
-- 
cgit v1.2.3