[libc-commits] [libc] [libc][math] Add shared functions to check exceptions for exp* functions. (PR #202503)

[via libc-commits]

================
@@ -0,0 +1,84 @@
+//===-- Check exceptions for exp functions ----------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_CHECK_EXP_EXCEPTIONS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_CHECK_EXP_EXCEPTIONS_H
+
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace check {
+
+namespace exp_internal {
+
+template <typename T> struct Bounds;
+
+template <> struct Bounds<float> {
+  // Smallest value that will cause overflow, generated from Sollya:
+  // > float_max = 2^127 * (2 - 2^-23);
+  // > upper = round(log(float_max), SG, RU);
+  // > printfloat(upper);
+  static constexpr float UPPER = 0x1.62e43p6f;
+  static constexpr uint32_t UPPER_BITS = 0x42b1'7218;
+  // Largest value that will cause underflow, generated from Sollya:
+  // > float_min = 2^-126;
+  // > lower = round(log(float_min), SG, RD);
+  // > printfloat(lower);
+  static constexpr float LOWER = -0x1.5d58ap6f;
+  static constexpr uint32_t LOWER_BITS = 0xc2ae'ac50;
+};
+
+template <> struct Bounds<double> {
+  // Smallest value that will cause overflow, generated from Sollya:
+  // > double_max = 2^1023 * (2 - 2^-52);
+  // > upper = round(log(double_max), D, RU);
+  // > printdouble(upper);
+  static constexpr double UPPER = 0x1.62e42fefa39fp9;
+  static constexpr uint64_t UPPER_BITS = 0x4086'2e42'fefa'39f0;
+  // Largest value that will cause underflow, generated from Sollya:
+  // > double_min = 2^-1022;
+  // > lower = round(log(double_min), D, RD);
+  // > printfloat(lower);
+  static constexpr double LOWER = -0x1.6232bdd7abcd3p9;
+  static constexpr uint64_t LOWER_BITS = 0xc086'232b'dd7a'bcd3;
+};
+
+} // namespace exp_internal
+
+template <typename T> LIBC_INLINE int exp_exceptions(T x, int rounding_mode) {
+  using FPBits = typename fputil::FPBits<T>;
+  using StorageType = typename FPBits::StorageType;
+
+  FPBits x_bits(x);
+  if (x_bits.is_signaling_nan())
+    return FE_INVALID;
+  if (x_bits.is_inf_or_nan() || x_bits.is_zero())
+    return 0;
+  StorageType x_u = x_bits.uintval();
+  if (x_u >= exp_internal::Bounds<T>::UPPER_BITS && x_bits.is_pos())
+    return (rounding_mode == FE_TONEAREST || rounding_mode == FE_UPWARD)
----------------
lntue wrote:

> > In this example, the result is overflow with `FE_UPWARD` but not overflow with `FE_TONEAREST` because one is resulted in Inf, and the other is finite?
> 
> That is completely consistent with what I have been saying though.
> 
> Where infinity is an option, if a large result rounds to a finite value under `FE_TONEAREST`, then (by definition) the error between the finite result value and the true value is at most 0.5ulp of the true value.
> 
> Rounding from a finite true value to infinity has an error greater than 0.5ulp.

So the point for that example is that, the standard says that the overflow is raised regardless of the rounding modes.  But in here, you don't want to raise overflow exception for `FE_TONEAREST`, it is within 0.5 ulp of the true value, while overflow should be raised for `FE_UPWARD`.

I guess, the general question would be whether we want to raise exceptions based on pre-rounding / real values, or post-rounding values.  For underflow, the differences are not too much, but in overflow, the differences between pre/post rounding could be very large.

For most of libm implementations, using plain arithmetic does have performance advantage over touching floating-point environment registers, sometime quite significant.

https://github.com/llvm/llvm-project/pull/202503