[libc] Add float type and flag for nearest_integer to enable SSE4.2.

[lntue]

Summary

Add float type and flag for nearest integer to automatically test with
and without SSE4.2 flag.

[sivachandra]

libc/src/__support/FPUtil/nearest_integer.h

29

Can you explain where this overload is to be used? I am asking because I see couple of problems here:

1. The constants that you used are of double type. So, the RHS expression on line 33 will get evaluated as a double expression and will be equivalent to r = x for a large range of x values.
2. Even if you add a suffix of f to the constants, the behavior is still dependent on the value of FLT_EVAL_METHOD: https://en.cppreference.com/w/c/types/limits/FLT_EVAL_METHOD. Perhaps this aspect is not of a major concern on platforms like x86_64 and aarch64 anyway.

Going by #1 at least, I am surprised that this fallback implementation is serving as a drop-in-replacement to the x86_64 and aarch64 specializations. If a drop-in-replacement is in fact the intention, then the testing strategy should be this:

1. Unit test this fallback implementation separately - testing with the user of this fallback will amount to a kind of integration test.
2. Unit test the users with either the fallback or the preferred path. There is not much gain in unit testing the users for both the preferred path and the fallback path.
3. If testing of the plumbing of the fallback/preferred path is desired, then testing for that should be done with explicit intent.

One more point that strikes me here is the comment you have above, "... in case of a tie, might pick a random one among 2 closest integers when the rounding mode is not FE_TONEAREST." If the users' algorithms are really tolerant to such "randomness", then perhaps we should give this fallback a name which captures that randomness aspect. Also, if there is only one user of this fallback, then may be this fallback should live with that user and not as utility?

[lntue]

libc/src/__support/FPUtil/nearest_integer.h

29

Thanks for spotting the data type issue!

About the intention of the change, this is to be used for improving exp*f functions as experimenting in https://reviews.llvm.org/D130008. The main thing that these functions are used to to replace the idiom for rounding to nearest integer that we've been using:

int k = static_cast<int>(x < 0 ? x - 0.5f : x + 0.5f);
float kf = static_cast<float>(k);

This by itself also has all the problems that you mentioned in the comment, such as FLT_EVAL_METHOD, possible different answers for different rounding modes, etc. Maybe I shouldn't use the word random. It well-defined, it's just too wordy to describe its behavior explicitly in all cases.

About testing, this function is similar to fputil::multiply_add and fputil::polyeval that might give slightly different answers with/without special instructions, and that the fallback path is never taken on aarch64. So it's suitable for a flag control in the same way:

• The end goal for each math function is that the final outputs are identical regardless of which path nearest_integer takes, so technically, we should require both preferred path and fallback path to be tested, unless explicitly stated otherwise.

So to answer to your questions: this function is kind of similar to polyeval / multiply_add, and will be used by several functions, hence being shared here.

[lntue]

Update comments and float literals.

[sivachandra]

libc/src/__support/FPUtil/nearest_integer.h

29

Thanks for spotting the data type issue!

About the intention of the change, this is to be used for improving exp*f functions as experimenting in https://reviews.llvm.org/D130008. The main thing that these functions are used to to replace the idiom for rounding to nearest integer that we've been using:

int k = static_cast<int>(x < 0 ? x - 0.5f : x + 0.5f);
float kf = static_cast<float>(k);

For my knowledge, can you explain what the problem is with this idiom and how the proposed implementation is solving it?

This by itself also has all the problems that you mentioned in the comment, such as FLT_EVAL_METHOD, possible different answers for different rounding modes, etc. Maybe I shouldn't use the word random. It well-defined, it's just too wordy to describe its behavior explicitly in all cases.

Considering that this fallback is not equivalent to the behavior of the native instructions, and that it was previously erroneous for some inputs, can I conclude that exp*f you are experimenting with is tolerant to such incorrectness?

About testing, this function is similar to fputil::multiply_add and fputil::polyeval that might give slightly different answers with/without special instructions, and that the fallback path is never taken on aarch64. So it's suitable for a flag control in the same way:

• The end goal for each math function is that the final outputs are identical regardless of which path nearest_integer takes, so technically, we should require both preferred path and fallback path to be tested, unless explicitly stated otherwise.

So to answer to your questions: this function is kind of similar to polyeval / multiply_add, and will be used by several functions, hence being shared here.

[lntue]

libc/src/__support/FPUtil/nearest_integer.h

29

Thanks for spotting the data type issue!

About the intention of the change, this is to be used for improving exp*f functions as experimenting in https://reviews.llvm.org/D130008. The main thing that these functions are used to to replace the idiom for rounding to nearest integer that we've been using:

int k = static_cast<int>(x < 0 ? x - 0.5f : x + 0.5f);
float kf = static_cast<float>(k);

For my knowledge, can you explain what the problem is with this idiom and how the proposed implementation is solving it?

This idiom first converting a floating point type to an integer type, and then converting back to a floating point. CPUs seem to hate this ping pong back and forth between floating point and integer registers, significantly reduce the throughput and increase latency with its on the critical path, which is mostly floating point computations. The key point that makes this implementation improves the performance is that the entire computation is within floating point type and not too many branches.

This by itself also has all the problems that you mentioned in the comment, such as FLT_EVAL_METHOD, possible different answers for different rounding modes, etc. Maybe I shouldn't use the word random. It well-defined, it's just too wordy to describe its behavior explicitly in all cases.

Considering that this fallback is not equivalent to the behavior of the native instructions, and that it was previously erroneous for some inputs, can I conclude that exp*f you are experimenting with is tolerant to such incorrectness?

Yes, it can kind of tolerate these different inputs (I wouldn't call it incorrect, since technically they are all correct, just different decisions / rounding modes when there is a tie), probably few exceptional values need to be updated similar to sinf.

libc/cmake/modules/LLVMLibCFlagRules.cmake

else()
       list(APPEND fq_dep_no_flag_list ${fq_dep_name})
     endif()
   endforeach(dep)
   set(${output_list} ${fq_dep_no_flag_list} PARENT_SCOPE)
 endfunction(get_fq_dep_list_without_flag)
 
 # Special flags
 set(FMA_OPT_FLAG "FMA_OPT")
+set(ROUND_OPT_FLAG "ROUND_OPT")
 
 # Skip FMA_OPT flag for targets that don't support fma.
 if(NOT(LIBC_TARGET_ARCHITECTURE_IS_X86 AND (LIBC_CPU_FEATURES MATCHES "FMA")))
   set(SKIP_FLAG_EXPANSION_FMA_OPT TRUE)
 endif()
+
+# Skip ROUND_OPT flag for targets that don't support SSE 4.2.
+if(NOT(LIBC_TARGET_ARCHITECTURE_IS_X86 AND (LIBC_CPU_FEATURES MATCHES "SSE4_2")))
+  set(SKIP_FLAG_EXPANSION_ROUND_OPT TRUE)
+endif()

libc/cmake/modules/LLVMLibCObjectRules.cmake

 set(OBJECT_LIBRARY_TARGET_TYPE "OBJECT_LIBRARY")
 
 function(_get_common_compile_options output_var flags)
   list(FIND flags ${FMA_OPT_FLAG} fma)
   if(${fma} LESS 0)
     list(FIND flags "${FMA_OPT_FLAG}__ONLY" fma)
   endif()
   if((${fma} GREATER -1) AND (LIBC_CPU_FEATURES MATCHES "FMA"))
     set(ADD_FMA_FLAG TRUE)
   endif()
 
+  list(FIND flags ${ROUND_OPT_FLAG} round)
+  if(${round} LESS 0)
+    list(FIND flags "${ROUND_OPT_FLAG}__ONLY" round)
+  endif()
+  if((${round} GREATER -1) AND (LIBC_CPU_FEATURES MATCHES "SSE4_2"))
+    set(ADD_SSE4_2_FLAG TRUE)
+  endif()
+
   set(compile_options ${LIBC_COMPILE_OPTIONS_DEFAULT} ${ARGN})
   if(NOT ${LIBC_TARGET_OS} STREQUAL "windows")
     set(compile_options ${compile_options} -fpie -ffreestanding -fno-builtin)
   endif()
   if(LLVM_COMPILER_IS_GCC_COMPATIBLE)
     list(APPEND compile_options "-fno-exceptions")
     list(APPEND compile_options "-fno-unwind-tables")
     list(APPEND compile_options "-fno-asynchronous-unwind-tables")
     list(APPEND compile_options "-fno-rtti")
     if(ADD_FMA_FLAG)
       list(APPEND compile_options "-mfma")
     endif()
+    if(ADD_SSE4_2_FLAG)
+      list(APPEND compile_options "-msse4.2")
+    endif()
   elseif(MSVC)
     list(APPEND compile_options "/EHs-c-")
     list(APPEND compile_options "/GR-")
     if(ADD_FMA_FLAG)
       list(APPEND compile_options "/arch:AVX2")
     endif()
   endif()
   set(${output_var} ${compile_options} PARENT_SCOPE)

libc/src/__support/FPUtil/CMakeLists.txt

 )
 
 add_header_library(
   nearest_integer
   HDRS
     nearest_integer.h
   DEPENDS
     libc.src.__support.common
+  FLAGS
+    ROUND_OPT
 )
 
 add_subdirectory(generic)

libc/src/__support/FPUtil/aarch64/nearest_integer.h

 
 #if !defined(LLVM_LIBC_ARCH_AARCH64)
 #error "Invalid include"
 #endif
 
 namespace __llvm_libc {
 namespace fputil {
 
+static inline float nearest_integer(float x) {
+  float result;
+  __asm__ __volatile__("frintn %s0, %s1\n\t" : "=w"(result) : "w"(x));
+  return result;
+}
+
 static inline double nearest_integer(double x) {
   double result;
   __asm__ __volatile__("frintn %d0, %d1\n\t" : "=w"(result) : "w"(x));
   return result;
 }
 
 } // namespace fputil
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_AARCH64_NEAREST_INTEGER_H

libc/src/__support/FPUtil/nearest_integer.h

 #include "x86_64/nearest_integer.h"
 #elif defined(LLVM_LIBC_ARCH_AARCH64)
 #include "aarch64/nearest_integer.h"
 #else
 
 namespace __llvm_libc {
 namespace fputil {
 
-// This is a fast implementation for rounding to a nearest integer that, in case
-// of a tie, might pick a random one among 2 closest integers when the rounding
-// mode is not FE_TONEAREST.
+// This is a fast implementation for rounding to a nearest integer that.
 //
 // Notice that for AARCH64 and x86-64 with SSE4.2 support, we will use their
 // corresponding rounding instruction instead.  And in those cases, the results
 // are rounded to the nearest integer, tie-to-even.
+static inline float nearest_integer(float x) {
+  if (x < 0x1p24f && x > -0x1p24f) {
+    float r = x < 0 ? (x - 0x1.0p23f) + 0x1.0p23f : (x + 0x1.0p23f) - 0x1.0p23f;
+    float diff = x - r;
+    // The expression above is correct for the default rounding mode, round-to-
+    // nearest, tie-to-even.  For other rounding modes, it might be off by 1,
+    // which is corrected below.
+    if (unlikely(diff > 0.5f))
+      return r + 1.0f;
+    if (unlikely(diff < -0.5f))
+      return r - 1.0f;
+    return r;
+  }
+  return x;
+}
+
 static inline double nearest_integer(double x) {
   if (x < 0x1p53 && x > -0x1p53) {
     double r = x < 0 ? (x - 0x1.0p52) + 0x1.0p52 : (x + 0x1.0p52) - 0x1.0p52;
     double diff = x - r;
     // The expression above is correct for the default rounding mode, round-to-
     // nearest, tie-to-even.  For other rounding modes, it might be off by 1,
     // which is corrected below.
     if (unlikely(diff > 0.5))

libc/src/__support/FPUtil/x86_64/nearest_integer.h

 #error "SSE4.2 instruction set is not supported"
 #endif
 
 #include <immintrin.h>
 
 namespace __llvm_libc {
 namespace fputil {
 
+static inline float nearest_integer(float x) {
+  __m128 xmm = _mm_set_ss(x); // NOLINT
+  __m128 ymm =
+      _mm_round_ss(xmm, xmm, _MM_ROUND_NEAREST | _MM_FROUND_NO_EXC); // NOLINT
+  return ymm[0];
+}
+
 static inline double nearest_integer(double x) {
   __m128d xmm = _mm_set_sd(x); // NOLINT
   __m128d ymm =
       _mm_round_sd(xmm, xmm, _MM_ROUND_NEAREST | _MM_FROUND_NO_EXC); // NOLINT
   return ymm[0];
 }
 
 } // namespace fputil
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_NEAREST_INTEGER_H