Use __builtin_clz to find leading 1 in generic sqrt (where possible)

[cratonica]

Summary

__builtin_clz requires just a single instruction on x86 and arm, so this is a performance improvement.

[cratonica]

Fix some declaration ordering

[cratonica]

Attempting to fix patch

[lntue]

libc/src/__support/FPUtil/generic/sqrt.h

65–66

Can you do a simple perf test to see if using clz for 64 bits is faster than binary search? Something like:

uint64_t hi_bits = static_cast<uint64_t>(mantissa >> 64);
int shift  = hi_bits ? (clz(hi_bits) - 15) : (clz(static_cast<uint64_t>(mantissa)) + 49);
exponent -= shift;
mantissa <<= shift;

[cratonica]

libc/src/__support/FPUtil/generic/sqrt.h

65–66

1. Unfortunately there aren't any differential or performance tests for sqrt or sqrtl, so I'll need to add those first in a separate PR. It won't be too much work, just clone the ones for sqrtf.
2. Also, I don't actually have a machine on which I can test 128-bit floats -- my machine uses the 80-bit x87 format, and aarch64 uses 64-bit for long double.

[cratonica]

libc/src/__support/FPUtil/generic/sqrt.h

65–66

Following up here:

re: #1: It seems as though we only have perf tests for the float32 variants of math functions due to the logarithmic increase in the domain required for 64-bit inputs. Trying to run an exhaustive performance test using float64 never completed even after an hour of waiting, and I can't imagine an exhaustive test for 128-bit inputs would complete even after days. So I'm going to write a one-off performance test that terminates after 2^24 iterations to test this, but I won't be checking it in

re: #2: The x87 8-bit variant uses a 64-bit mantissa, which means that __bulting_clzll is will still work after truncation, so this is trivial to implement (and I have confirmed a slight performance increase here using the method mentioned above). I was incorrect that aarch64 uses 64-bit floats for long double, and I have access to some hardware with an aarch64 Cortex-A53 that I can run these performance test on with the changes you mentioned. If performance is improved, then I will update the patch accordingly.

[cratonica]

libc/src/__support/FPUtil/generic/sqrt.h

65–66

Results for 128-bit floats (long double) from the aarch64 Cortex-A53 core in denormal range:

clz: 28894915309ns
binary search: 29253929397ns

So, just over a 1% performance improvement, which is in-line with what I'm seeing on the 32-bit float sqrtf function.

Therefore, I'm patching in that change (as well as the x87 80-bit specialization).

[cratonica]

Add optimizations for x87 80-bit floats and aarch64 128-bit floats

[lntue]

Thanks for doing the performance tests for all the data types!

[cratonica]

Thank you! As I don't have repo access (nor probably should I yet), please submit per https://llvm.org/docs/Phabricator.html#committing-someone-s-change-from-phabricator if you would (sivachandra did this for me on D117684).

libc/src/__support/FPUtil/generic/sqrt.h

 };
 
 #if defined(SPECIAL_X86_LONG_DOUBLE)
 template <> struct SpecialLongDouble<long double> {
   static constexpr bool VALUE = true;
 };
 #endif // SPECIAL_X86_LONG_DOUBLE
 
-template <typename T>
-static inline void normalize(int &exponent,
-                             typename FPBits<T>::UIntType &mantissa);
-
-template <> inline void normalize<float>(int &exponent, uint32_t &mantissa) {
-  // Use binary search to shift the leading 1 bit.
-  // With MantissaWidth<float> = 23, it will take
-  // ceil(log2(23)) = 5 steps checking the mantissa bits as followed:
-  // Step 1: 0000 0000 0000 XXXX XXXX XXXX
-  // Step 2: 0000 00XX XXXX XXXX XXXX XXXX
-  // Step 3: 000X XXXX XXXX XXXX XXXX XXXX
-  // Step 4: 00XX XXXX XXXX XXXX XXXX XXXX
-  // Step 5: 0XXX XXXX XXXX XXXX XXXX XXXX
-  constexpr int NSTEPS = 5; // = ceil(log2(MantissaWidth))
-  constexpr uint32_t BOUNDS[NSTEPS] = {1 << 12, 1 << 18, 1 << 21, 1 << 22,
-                                       1 << 23};
-  constexpr int SHIFTS[NSTEPS] = {12, 6, 3, 2, 1};
-
-  for (int i = 0; i < NSTEPS; ++i) {
-    if (mantissa < BOUNDS[i]) {
-      exponent -= SHIFTS[i];
-      mantissa <<= SHIFTS[i];
-    }
-  }
+// The following overloads are matched based on what is accepted by
+// __builtin_clz* rather than using the exactly-sized aliases from stdint.h.
+// This way, we can avoid making any assumptions about integer sizes and let the
+// compiler match for us.
+template <typename T> static inline int clz(T val);
+template <> inline int clz<unsigned int>(unsigned int val) {
+  return __builtin_clz(val);
+template <> inline int clz<unsigned long int>(unsigned long int val) {
+  return __builtin_clzl(val);
+template <> inline int clz<unsigned long long int>(unsigned long long int val) {
+  return __builtin_clzll(val);
 }
 
-template <> inline void normalize<double>(int &exponent, uint64_t &mantissa) {
-  // Use binary search to shift the leading 1 bit similar to float.
-  // With MantissaWidth<double> = 52, it will take
-  // ceil(log2(52)) = 6 steps checking the mantissa bits.
-  constexpr int NSTEPS = 6; // = ceil(log2(MantissaWidth))
-  constexpr uint64_t BOUNDS[NSTEPS] = {1ULL << 26, 1ULL << 39, 1ULL << 46,
-                                       1ULL << 49, 1ULL << 51, 1ULL << 52};
-  constexpr int SHIFTS[NSTEPS] = {27, 14, 7, 4, 2, 1};
-
-  for (int i = 0; i < NSTEPS; ++i) {
-    if (mantissa < BOUNDS[i]) {
-      exponent -= SHIFTS[i];
-      mantissa <<= SHIFTS[i];
-    }
-  }
+template <typename T>
+static inline void normalize(int &exponent,
+                             typename FPBits<T>::UIntType &mantissa) {
+  const int shift =
+      clz(mantissa) - (8 * sizeof(mantissa) - 1 - MantissaWidth<T>::VALUE);
+  exponent -= shift;
+  mantissa <<= shift;
 }
 
 #ifdef LONG_DOUBLE_IS_DOUBLE
 template <>
 inline void normalize<long double>(int &exponent, uint64_t &mantissa) {
   normalize<double>(exponent, mantissa);
 }
 #elif !defined(SPECIAL_X86_LONG_DOUBLE)
 template <>
 inline void normalize<long double>(int &exponent, __uint128_t &mantissa) {
-  // Use binary search to shift the leading 1 bit similar to float.
-  // With MantissaWidth<long double> = 112, it will take
-  // ceil(log2(112)) = 7 steps checking the mantissa bits.
-  constexpr int NSTEPS = 7; // = ceil(log2(MantissaWidth))
-  constexpr __uint128_t BOUNDS[NSTEPS] = {
-      __uint128_t(1) << 56,  __uint128_t(1) << 84,  __uint128_t(1) << 98,
-      __uint128_t(1) << 105, __uint128_t(1) << 109, __uint128_t(1) << 111,
-      __uint128_t(1) << 112};
-  constexpr int SHIFTS[NSTEPS] = {57, 29, 15, 8, 4, 2, 1};
-
-  for (int i = 0; i < NSTEPS; ++i) {
-    if (mantissa < BOUNDS[i]) {
-      exponent -= SHIFTS[i];
-      mantissa <<= SHIFTS[i];
-    }
-  }
+  const uint64_t hi_bits = static_cast<uint64_t>(mantissa >> 64);
+  const int shift = hi_bits ? (clz(hi_bits) - 15)
+                            : (clz(static_cast<uint64_t>(mantissa)) + 49);
+  exponent -= shift;
+  mantissa <<= shift;
 }
 #endif
 
 } // namespace internal
 
 // Correctly rounded IEEE 754 SQRT for all rounding modes.
 // Shift-and-add algorithm.
 template <typename T>

libc/src/__support/FPUtil/generic/sqrt_80_bit_long_double.h

 #include "src/__support/FPUtil/FPBits.h"
 #include "src/__support/FPUtil/PlatformDefs.h"
 
 namespace __llvm_libc {
 namespace fputil {
 namespace x86 {
 
 inline void normalize(int &exponent, __uint128_t &mantissa) {
-  // Use binary search to shift the leading 1 bit similar to float.
-  // With MantissaWidth<long double> = 63, it will take
-  // ceil(log2(63)) = 6 steps checking the mantissa bits.
-  constexpr int NSTEPS = 6; // = ceil(log2(MantissaWidth))
-  constexpr __uint128_t BOUNDS[NSTEPS] = {
-      __uint128_t(1) << 32, __uint128_t(1) << 48, __uint128_t(1) << 56,
-      __uint128_t(1) << 60, __uint128_t(1) << 62, __uint128_t(1) << 63};
-  constexpr int SHIFTS[NSTEPS] = {32, 16, 8, 4, 2, 1};
-
-  for (int i = 0; i < NSTEPS; ++i) {
-    if (mantissa < BOUNDS[i]) {
-      exponent -= SHIFTS[i];
-      mantissa <<= SHIFTS[i];
-    }
-  }
+  const int shift =
+      __builtin_clzll(static_cast<uint64_t>(mantissa)) -
+      (8 * sizeof(uint64_t) - 1 - MantissaWidth<long double>::VALUE);
+  exponent -= shift;
+  mantissa <<= shift;
 }
 
 // if constexpr statement in sqrt.h still requires x86::sqrt to be declared
 // even when it's not used.
 static inline long double sqrt(long double x);
 
 // Correctly rounded SQRT for all rounding modes.
 // Shift-and-add algorithm.