Fix the setting of the rounding mode

luhenry · luhenry · commit 800cc558640b · 2024-02-14T07:15:20.000Z
It would previously read the rounding mode, but wouldn't set it when it
was not TONEAREST.

That would work in most cases because TONEAREST is the default rounding
mode in a process, and none of our tests would set it differently
anyway.
diff --git a/include/rvvlm.h b/include/rvvlm.h
@@ -2,6 +2,8 @@
 //
 // SPDX-License-Identifier: Apache-2.0
 
+#pragma once
+
 #include <stdbool.h>
 #include <stdint.h>
 
@@ -35,18 +37,41 @@ union sui64_fp64 {
 #define UNIT_STRIDE 1
 #define GENERAL_STRIDE 2
 
+#ifndef CSR_FRM
+#define CSR_FRM       0x002
+#endif
+
+#ifndef FE_TONEAREST
+#define FE_TONEAREST  0x000
+#endif
+
+
+#define read_frm()                                                             \
+({                                                                             \
+  unsigned long __value;                                                       \
+  __asm__ __volatile__ (                                                       \
+    "frrm %0" : "=r" (__value) :: "memory");                                   \
+  __value;                                                                     \
+})
+
+#define write_frm(value)                                                        \
+({                                                                             \
+  unsigned long __value;                                                       \
+  __asm__ __volatile__ (                                                       \
+    "fsrm %0, %1" : "=r" (__value) : "r" (CSR_FRM) : "memory");                \
+  __value;                                                                     \
+})
+
 #define SET_ROUNDTONEAREST                                                     \
-  int original_frm;                                                            \
-  bool need_to_restore;                                                        \
-  do {                                                                         \
-    (original_frm) = fegetround();                                             \
-    need_to_restore = (original_frm != FE_TONEAREST);                          \
-  } while (0)
+  int __original_frm = read_frm();                                             \
+  if (__original_frm != FE_TONEAREST) {                                        \
+    write_frm(FE_TONEAREST);                                                   \
+  }
 
 #define RESTORE_FRM                                                            \
   do {                                                                         \
-    if (need_to_restore) {                                                     \
-      fesetround((original_frm));                                              \
+    if (__original_frm != FE_TONEAREST) {                                      \
+      write_frm(__original_frm);                                               \
     }                                                                          \
   } while (0)
 
diff --git a/include/rvvlm_acoshD.inc.h b/include/rvvlm_acoshD.inc.h
@@ -10,8 +10,6 @@
 #define F_VER1 RVVLM_ACOSHDI_STD
 #endif
 
-#include <fenv.h>
-
 // Acosh(x)  is defined for x >= 1 by the formula log(x + sqrt(x*x - 1))
 // and for the log function log(2^n z), we uses the expansion in terms of atanh
 // n log(2) + 2 atanh((z-1)/(z+1))
diff --git a/include/rvvlm_asincosD.inc.h b/include/rvvlm_asincosD.inc.h
@@ -2,8 +2,6 @@
 //
 // SPDX-License-Identifier: Apache-2.0
 
-#include <fenv.h>
-
 #if defined(COMPILE_FOR_ASIN)
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_ASIND_FIXEDPT
@@ -116,8 +114,6 @@ static_assert(false, "Must specify asin, acos, asinpi or acospi" __FILE__);
     }                                                                          \
   } while (0)
 
-#include <fenv.h>
-
 // For asin/acos, the computation is of the form  Const +/- (r + r*s*poly(s))
 // This version computes this entire expression in fixed point by converting
 // r and s into fixed point.
diff --git a/include/rvvlm_asinhcoshD.inc.h b/include/rvvlm_asinhcoshD.inc.h
@@ -18,8 +18,6 @@
 #endif
 #endif
 
-#include <fenv.h>
-
 // Acosh(x)  is defined for x >= 1 by the formula log(x + sqrt(x*x - 1))
 // Asinh(x) is defined for all finite x by the formula log(x + sqrt(x*x + 1))
 // Acosh is always positive, and Asinh(-x) = -Asinh(x). Thus we in general
diff --git a/include/rvvlm_atan2D.inc.h b/include/rvvlm_atan2D.inc.h
@@ -2,8 +2,6 @@
 //
 // SPDX-License-Identifier: Apache-2.0
 
-#include <fenv.h>
-
 #if defined(COMPILE_FOR_ATAN2)
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_ATAN2D_FIXEDPT
diff --git a/include/rvvlm_atanD.inc.h b/include/rvvlm_atanD.inc.h
@@ -2,8 +2,6 @@
 //
 // SPDX-License-Identifier: Apache-2.0
 
-#include <fenv.h>
-
 #if defined(COMPILE_FOR_ATAN)
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_ATAND_FIXEDPT
@@ -71,8 +69,6 @@ static_assert(false, "Must specify atan or atanpi" __FILE__);
     }                                                                          \
   } while (0)
 
-#include <fenv.h>
-
 // For atan, atan(x) ~=~ r + r*s*poly(s), r = x and s = r*r for |x| < 1
 // and atan(x) = pi/2 - atan(1/x) for |x| >= 1
 // Thus atan(x) = (pi/2 or 0) +/- (r + r*s*poly(s)), where r is x or 1/x, s is
diff --git a/include/rvvlm_atanhD.inc.h b/include/rvvlm_atanhD.inc.h
@@ -10,8 +10,6 @@
 #define F_VER1 RVVLM_ATANHDI_MIXED
 #endif
 
-#include <fenv.h>
-
 // Atanh(x) is defined only for |x| <= 1. As atanh(-x) = -atanh(x), the
 // main computation works with |x|.
 // For |x| > 1 and x being a sNaN, the invalid signal has to be generated
diff --git a/include/rvvlm_cbrtD.inc.h b/include/rvvlm_cbrtD.inc.h
@@ -2,16 +2,12 @@
 //
 // SPDX-License-Identifier: Apache-2.0
 
-#include <fenv.h>
-
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_CBRTD_ITER
 #else
 #define F_VER1 RVVLM_CBRTDI_ITER
 #endif
 
-#include <fenv.h>
-
 #define EXCEPTION_HANDLING_CBRT(vx, special_args, vy_special, n_adjust, vlen)  \
   do {                                                                         \
     VUINT vclass = __riscv_vfclass((vx), (vlen));                              \
diff --git a/include/rvvlm_expD.inc.h b/include/rvvlm_expD.inc.h
@@ -86,8 +86,6 @@
 static_assert(false, "Must specify base of exponential" __FILE__);
 #endif
 
-#include <fenv.h>
-
 // Version 1 is reduction to standard primary interval.
 // Reduced argument is represented as one FP64 variable.
 void F_VER1(API) {
diff --git a/include/rvvlm_expm1D.inc.h b/include/rvvlm_expm1D.inc.h
@@ -32,8 +32,6 @@
 #define X_MAX 0x1.65p+9
 #define X_MIN -0x1.5p+5
 
-#include <fenv.h>
-
 // We use the EPsim version of expD to compute expm1
 void F_VER1(API) {
   size_t vlen;
diff --git a/include/rvvlm_log1pD.inc.h b/include/rvvlm_log1pD.inc.h
@@ -36,8 +36,6 @@
 #define LOG2_HI 0x1.62e42fefa39efp-1
 #define LOG2_LO 0x1.abc9e3b39803fp-56
 
-#include <fenv.h>
-
 // Version 1 uses a 128-entry LUT
 void F_VER1(API) {
   size_t vlen;
diff --git a/include/rvvlm_logD.inc.h b/include/rvvlm_logD.inc.h
@@ -72,8 +72,6 @@
 static_assert(false, "Must specify base of logarithm" __FILE__);
 #endif
 
-#include <fenv.h>
-
 // Version 1 uses a 128-entry LUT
 void F_VER1(API) {
   size_t vlen;
diff --git a/include/rvvlm_powD.inc.h b/include/rvvlm_powD.inc.h
@@ -181,8 +181,6 @@ static const double negtwo_to_65 = -0x1.0p65;
 #define F_VER1 RVVLM_POWDI_TBL
 #endif
 
-#include <fenv.h>
-
 // Version 1 is reduction to standard primary interval.
 // Reduced argument is represented as one FP64 variable.
 void F_VER1(API) {
diff --git a/include/rvvlm_sinandcosD.inc.h b/include/rvvlm_sinandcosD.inc.h
@@ -4,8 +4,6 @@
 
 #include "rvvlm_trigD.h"
 
-#include <fenv.h>
-
 #if defined(COMPILE_FOR_SINCOS)
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_SINCOSD_STD
@@ -22,8 +20,6 @@
 static_assert(false, "Must specify sincos or sincospi" __FILE__);
 #endif
 
-#include <fenv.h>
-
 // This versions reduces argument to [-pi/4, pi/4] and computes sin(r) and
 // cos(r)
 void F_VER1(API) {
diff --git a/include/rvvlm_sincosD.inc.h b/include/rvvlm_sincosD.inc.h
@@ -4,8 +4,6 @@
 
 #include "rvvlm_trigD.h"
 
-#include <fenv.h>
-
 #if defined(COMPILE_FOR_SIN)
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_SIND_MERGED
@@ -34,8 +32,6 @@
 static_assert(false, "Must specify sin, sinpi, cos or cospi" __FILE__);
 #endif
 
-#include <fenv.h>
-
 // This versions reduces argument to [-pi/4, pi/4] and computes sin(r) or cos(r)
 // by merging the appropriate coefficients into a vector register
 void F_VER1(API) {
diff --git a/include/rvvlm_sinhcoshD.inc.h b/include/rvvlm_sinhcoshD.inc.h
@@ -4,8 +4,6 @@
 
 #include "rvvlm_hyperbolicsD.h"
 
-#include <fenv.h>
-
 #if defined(COMPILE_FOR_SINH)
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_SINHD_STD
@@ -26,8 +24,6 @@
 static_assert(false, "Must specify sinh or cosh" __FILE__);
 #endif
 
-#include <fenv.h>
-
 // This versions reduces argument to [-log2/2, log2/2]
 // Exploit common expressions exp(R) and exp(-R), and uses purely
 // floating point method to preserve precision
diff --git a/include/rvvlm_tanD.inc.h b/include/rvvlm_tanD.inc.h
@@ -4,8 +4,6 @@
 
 #include "rvvlm_trigD.h"
 
-#include <fenv.h>
-
 #if defined(COMPILE_FOR_TAN)
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_TAND_MERGED
@@ -22,8 +20,6 @@
 static_assert(false, "Must specify tan or tanpi" __FILE__);
 #endif
 
-#include <fenv.h>
-
 // This versions reduces argument to [-pi/4, pi/4] and computes sin(r) or cos(r)
 // tan(x) is either sin(r)/cos(r) or -cos(r)/sin(r)
 void F_VER1(API) {
diff --git a/include/rvvlm_tanhD.inc.h b/include/rvvlm_tanhD.inc.h
@@ -4,16 +4,12 @@
 
 #include "rvvlm_hyperbolicsD.h"
 
-#include <fenv.h>
-
 #if (STRIDE == UNIT_STRIDE)
 #define F_VER1 RVVLM_TANHD_STD
 #else
 #define F_VER1 RVVLM_TANHDI_STD
 #endif
 
-#include <fenv.h>
-
 // This versions reduces argument to [-log2/2, log2/2]
 // Exploit common expressions exp(R) and exp(-R),
 // and uses purely floating-point computation
