From 3232e73c446a3a70fd2fcce4eabc0564b08312a7 Mon Sep 17 00:00:00 2001
From: Georg-Johann Lay <avr@gjlay.de>
Date: Tue, 14 Nov 2023 12:05:19 +0100
Subject: [PATCH] LibF7: sinh: Fix loss of precision due to cancellation for
 small values.

libgcc/config/avr/libf7/
	* libf7-const.def [F7MOD_sinh_]: Add MiniMax polynomial.
	* libf7.c (f7_sinh): Use it instead of (exp(x) - exp(-x)) / 2
	when |x| < 0.5 to avoid loss of precision due to cancellation.
---
 libgcc/config/avr/libf7/libf7-const.def | 10 ++++++++++
 libgcc/config/avr/libf7/libf7.c         | 17 +++++++++++++++++
 2 files changed, 27 insertions(+)

diff --git a/libgcc/config/avr/libf7/libf7-const.def b/libgcc/config/avr/libf7/libf7-const.def
index 0e4c4d8701eb..f772adb1262b 100644
--- a/libgcc/config/avr/libf7/libf7-const.def
+++ b/libgcc/config/avr/libf7/libf7-const.def
@@ -194,5 +194,15 @@ F7_CONST_DEF (X, 1, 0xc7,0xb5,0x6a,0xf8,0x0e,0x32,0x07, -37)
 F7_CONST_DEF (pi_low,0, 0xd3,0x13,0x19,0x8a,0x2e,0x03,0x70, 1 - F7_MANT_BITS-2)
 #endif
 
+#elif defined (F7MOD_sinh_)
+// MiniMax for sinh(q)/q, q = sqrt(x) for q in [0, 0.2505].  Quality pQ10 > 70.
+// 0.99999999999999998094379 + 0.16666666666667217765428 x + 0.0083333333330755574996361 x^2 + 1.9841270281701916844502e-4 x^3 + 2.7556979384534689658282e-6 x^4 + 2.5172859336028750964929e-8 x^5
+F7_CONST_DEF (X, 0, 0xff,0xff,0xff,0xff,0xff,0xff,0xff, -1)
+F7_CONST_DEF (X, 0, 0xaa,0xaa,0xaa,0xaa,0xaa,0xb0,0xdf, -3)
+F7_CONST_DEF (X, 0, 0x88,0x88,0x88,0x88,0x76,0x64,0xdb, -7)
+F7_CONST_DEF (X, 0, 0xd0,0x0d,0x01,0x1d,0x88,0x4c,0xed, -13)
+F7_CONST_DEF (X, 0, 0xb8,0xee,0x87,0xb4,0x30,0xf0,0xa1, -19)
+F7_CONST_DEF (X, 0, 0xd8,0x3b,0xb3,0xfd,0x9e,0x6c,0xcf, -26)
+
 #endif
 #endif /* ! IN_LIBF7_H && ! F7MOD_const_ */
diff --git a/libgcc/config/avr/libf7/libf7.c b/libgcc/config/avr/libf7/libf7.c
index da2a4b61b746..bf1cd140bd4a 100644
--- a/libgcc/config/avr/libf7/libf7.c
+++ b/libgcc/config/avr/libf7/libf7.c
@@ -2022,9 +2022,26 @@ void f7_sinhcosh (f7_t *cc, const f7_t *aa, bool do_sinh)
 
 
 #ifdef F7MOD_sinh_
+
+#define ARRAY_NAME coeff_sinh
+#include "libf7-array.def"
+#undef ARRAY_NAME
+
 F7_WEAK
 void f7_sinh (f7_t *cc, const f7_t *aa)
 {
+  if (aa->expo <= -2)
+    {
+      // For small values, exp(A) - exp(-A) suffers from cancellation, hence
+      // use a MiniMax polynomial for |A| < 0.5.
+      f7_t xx7, *xx = &xx7;
+      f7_t hh7, *hh = &hh7;
+      f7_square (xx, aa);
+      f7_horner (hh, xx, n_coeff_sinh, coeff_sinh, NULL);
+      f7_mul (cc, aa, hh);
+      return;
+    }
+
   f7_sinhcosh (cc, aa, true);
 }
 #endif // F7MOD_sinh_
-- 
GitLab