[compiler-rt] [builtins] Add logb/logbf/logbl methods to compiler-rt to avoid libm dependencies when possible.

[rupprecht]

Summary

The complex division builtins (div?c3) use logb methods from libm to scale numbers during division and avoid rounding issues. However, these come from libm, meaning anyone that uses --rtlib=compiler-rt also has to include -lm. Implement logb* methods for standard ieee 754 floats so we can avoid -lm on those platforms, falling back to the old behavior (using either logb() or __builtin_logb()) when not supported.

These new methods are defined internally as __compiler_rt_logb so as not to conflict with the libm definitions in any way.

This fixes just the libm methods mentioned in PR32279 and PR28652. libc is still required, although that seems to not be an issue.

Note: this is proposed as an alternative to just adding -lm: D49330.

[rupprecht]

Ping -- Eli, do you have time to review this, or add someone else that can take a look if you're busy?

[efriedma]

Maybe just make this a static inline function in fp_lib.h?

Maybe worth implementing something like ilogb instead, to avoid an extra int->float->int conversion? I guess it's not that important.

Fixes PR32279 and PR28652.

Is this really enough to fix PR32279? That report mentions some other symbols (scalbn*, fmaxl).

lib/builtins/ppc/divtc3.c

29–32

This is wrong. The original code is calling logb, not logbl.

[rupprecht]

In D49514#1172432, @efriedma wrote:

Maybe just make this a static inline function in fp_lib.h?

Thanks, done. That also makes these symbols not exported by the library, which is nice.

Maybe worth implementing something like ilogb instead, to avoid an extra int->float->int conversion? I guess it's not that important.

That sounds like a good idea, although the change is a little more invasive (e.g. all the div methods check against logb == infinity). I'd prefer to do that as a followup change if that's OK.

Fixes PR32279 and PR28652.

Is this really enough to fix PR32279? That report mentions some other symbols (scalbn*, fmaxl).

Hmm, this fixes it enough for us, but you're right this isn't providing scalbn*/fmax*. We must be adding a -lc automatically somewhere, which provides those methods. I'll update the patch description that this is only a partial fix.

lib/builtins/ppc/divtc3.c

29–32

This is wrong. The original code is calling logb, not logbl.

Yes, you're right.
I must have used logbl here because this is providing complex division for long doubles (c.f. non-ppc divtc3 uses logbl). But despite that, the method used is (non-long) doubles, due to ppc's floating point format here. I've fixed to use logb and left a comment to avoid that confusion.

[rupprecht]

• Inline methods to fp_lib.h and fix the ppc logb call

[efriedma]

lib/builtins/divxc3.c

16 ↗

(On Diff #157339)

Could you use a define that isn't QUAD_PRECISION here? Or just skip changing this file for now? This is the x86 80-bit long double.

lib/builtins/ppc/divtc3.c

8

Please change the comment so it doesn't refer to QUAD_PRECISION; that's specifically an IEEE 128-bit float, which isn't relevant here.

[rupprecht]

lib/builtins/divxc3.c

16 ↗

(On Diff #157339)

Could you use a define that isn't QUAD_PRECISION here? Or just skip changing this file for now? This is the x86 80-bit long double.

80-bit long double doesn't fit the abstractions used in fp_lib.h, so I'll just skip for now.

[rupprecht]

• Exclude 80-bit floating point and fix ppc comment

[rupprecht]

• Un-revert divtc3 (128-bit fp), revert divxc3 (80-bit fp)

lib/builtins/divxc3.c

16 ↗

(On Diff #157339)

Just realized I reverted the wrong one -- fixed so that I'm skipping divxc3 in this revision, but including divtc3 (both generic and ppc versions)

[efriedma]

LGTM, but I'd feel more comfortable if a second set of eyes looked at this.

[rupprecht]

Ping for @howard.hinnant @compnerd -- mind taking a quick look for a second LGTM?

[rupprecht]

Another ping for @howard.hinnant or @compnerd

[howard.hinnant]

• msg-22467-327.txt265 BDownload
• signature.asc848 BDownload

[rupprecht]

@compnerd, can you take a look? Or add another reviewer comfortable approving if you can't review?

[echristo]

Adding Steve here in case he wants to opine.

[compnerd]

lib/builtins/fp_lib.h

274

I don't see why we need __inline, I'm pretty sure that we compile in C99 mode, so inline should be fine. Doesn't this require GNU semantics though? Did you try building the windows target with this change?

302

Is this taken from somewhere else or is this code written from scratch?

323

Why the fallback to libm?

Additionally, the names are different, so is there a reason to not just define the variants unconditionally?

[rupprecht]

lib/builtins/fp_lib.h

274

This is consistent with all the other methods in this file.

Actually, it looks like the change from inline to __inline was from your revision: https://reviews.llvm.org/rL249953

302

This is written from scratch.

323

Why the fallback to libm?

This generic implementation only makes sense for ieee754, but some of the compiler-rt methods use other floating point types (e.g. divxc3 uses x86 80-bit long doubles). Hence, this implementation is incomplete and is using libm for now. (Note: this is not adding a dependency on libm, this is just calling out an existing one).

Additionally, the names are different, so is there a reason to not just define the variants unconditionally?

I originally had it split up like that (history link -- https://reviews.llvm.org/D49514?vs=on&id=156161&whitespace=ignore-most#toc -- things are split up into compiler_rt_logb.c etc.), but Eli suggested making it inline in an earlier comment. I think this way is a little nicer, but I don't have a strong opinion.

[rupprecht]

• Elaborate libm dependency in comment

[compnerd]

This seems reasonable to me, however, I would want to give @scanon a few (work) days to see if he has any comments on the actual implementation of logb.

[scanon]

I suspect that we'd rather use ilogb in the long term, but as a like-for-like replacement this looks OK.

[echristo]

Sounds like everyone is in favor. Hitting the accept button for Jordan :)

[itaraban]

It seems that http://lab.llvm.org:8011/builders/clang-x64-ninja-win7 buildbot is failing because of your commit.

[hans]

I've disabled the tests on Windows in r343095.

lib/builtins/divdc3.c

 /* ===-- divdc3.c - Implement __divdc3 -------------------------------------===
  *
  *                     The LLVM Compiler Infrastructure
  *
  * This file is dual licensed under the MIT and the University of Illinois Open
  * Source Licenses. See LICENSE.TXT for details.
  *
  * ===----------------------------------------------------------------------===
  *
  * This file implements __divdc3 for the compiler_rt library.
  *
  * ===----------------------------------------------------------------------===
  */
 
+#define DOUBLE_PRECISION
+#include "fp_lib.h"
 #include "int_lib.h"
 #include "int_math.h"
 
 /* Returns: the quotient of (a + ib) / (c + id) */
 
 COMPILER_RT_ABI Dcomplex
 __divdc3(double __a, double __b, double __c, double __d)
 {
     int __ilogbw = 0;
-    double __logbw = crt_logb(crt_fmax(crt_fabs(__c), crt_fabs(__d)));
+    double __logbw = __compiler_rt_logb(crt_fmax(crt_fabs(__c), crt_fabs(__d)));
     if (crt_isfinite(__logbw))
     {
         __ilogbw = (int)__logbw;
         __c = crt_scalbn(__c, -__ilogbw);
         __d = crt_scalbn(__d, -__ilogbw);
     }
     double __denom = __c * __c + __d * __d;
     Dcomplex z;

lib/builtins/divsc3.c

 /*===-- divsc3.c - Implement __divsc3 -------------------------------------===
  *
  *                     The LLVM Compiler Infrastructure
  *
  * This file is dual licensed under the MIT and the University of Illinois Open
  * Source Licenses. See LICENSE.TXT for details.
  *
  * ===----------------------------------------------------------------------===
  *
  * This file implements __divsc3 for the compiler_rt library.
  *
  *===----------------------------------------------------------------------===
  */
 
+#define SINGLE_PRECISION
+#include "fp_lib.h"
 #include "int_lib.h"
 #include "int_math.h"
 
 /* Returns: the quotient of (a + ib) / (c + id) */
 
 COMPILER_RT_ABI Fcomplex
 __divsc3(float __a, float __b, float __c, float __d)
 {
     int __ilogbw = 0;
-    float __logbw = crt_logbf(crt_fmaxf(crt_fabsf(__c), crt_fabsf(__d)));
+    float __logbw =
+        __compiler_rt_logbf(crt_fmaxf(crt_fabsf(__c), crt_fabsf(__d)));
     if (crt_isfinite(__logbw))
     {
         __ilogbw = (int)__logbw;
         __c = crt_scalbnf(__c, -__ilogbw);
         __d = crt_scalbnf(__d, -__ilogbw);
     }
     float __denom = __c * __c + __d * __d;
     Fcomplex z;

lib/builtins/divtc3.c

 /*===-- divtc3.c - Implement __divtc3 -------------------------------------===
  *
  *                     The LLVM Compiler Infrastructure
  *
  * This file is dual licensed under the MIT and the University of Illinois Open
  * Source Licenses. See LICENSE.TXT for details.
  *
  * ===----------------------------------------------------------------------===
  *
  * This file implements __divtc3 for the compiler_rt library.
  *
  *===----------------------------------------------------------------------===
  */
 
+#define QUAD_PRECISION
+#include "fp_lib.h"
 #include "int_lib.h"
 #include "int_math.h"
 
 /* Returns: the quotient of (a + ib) / (c + id) */
 
 COMPILER_RT_ABI Lcomplex
 __divtc3(long double __a, long double __b, long double __c, long double __d)
 {
     int __ilogbw = 0;
-    long double __logbw = crt_logbl(crt_fmaxl(crt_fabsl(__c), crt_fabsl(__d)));
+    long double __logbw =
+        __compiler_rt_logbl(crt_fmaxl(crt_fabsl(__c), crt_fabsl(__d)));
     if (crt_isfinite(__logbw))
     {
         __ilogbw = (int)__logbw;
         __c = crt_scalbnl(__c, -__ilogbw);
         __d = crt_scalbnl(__d, -__ilogbw);
     }
     long double __denom = __c * __c + __d * __d;
     Lcomplex z;

lib/builtins/fp_lib.h

 
 #ifndef FP_LIB_HEADER
 #define FP_LIB_HEADER
 
 #include <stdint.h>
 #include <stdbool.h>
 #include <limits.h>
 #include "int_lib.h"
+#include "int_math.h"
 
 // x86_64 FreeBSD prior v9.3 define fixed-width types incorrectly in
 // 32-bit mode.
 #if defined(__FreeBSD__) && defined(__i386__)
 # include <sys/param.h>
 # if __FreeBSD_version < 903000  // v9.3
 #  define uint64_t unsigned long long
 #  define int64_t long long
else if (count < 2*typeWidth) {
         *lo = *hi >> (count - typeWidth) | sticky;
         *hi = 0;
     } else {
         const bool sticky = *hi | *lo;
         *lo = sticky;
         *hi = 0;
     }
 }
+
+// Implements logb methods (logb, logbf, logbl) for IEEE-754. This avoids
+// pulling in a libm dependency from compiler-rt, but is not meant to replace
+// it (i.e. code calling logb() should get the one from libm, not this), hence
+// the __compiler_rt prefix.
+static __inline fp_t __compiler_rt_logbX(fp_t x) {
+  rep_t rep = toRep(x);
+  int exp = (rep & exponentMask) >> significandBits;
+
+  // Abnormal cases:
+  // 1) +/- inf returns +inf; NaN returns NaN
+  // 2) 0.0 returns -inf
+  if (exp == maxExponent) {
+    if (((rep & signBit) == 0) || (x != x)) {
+      return x;  // NaN or +inf: return x
+    } else {
+      return -x;  // -inf: return -x
+    }
+  } else if (x == 0.0) {
+    // 0.0: return -inf
+    return fromRep(infRep | signBit);
+  }
+
+  if (exp != 0) {
+    // Normal number
+    return exp - exponentBias;  // Unbias exponent
+  } else {
+    // Subnormal number; normalize and repeat
+    rep &= absMask;
+    const int shift = 1 - normalize(&rep);
+    exp = (rep & exponentMask) >> significandBits;
+    return exp - exponentBias - shift;  // Unbias exponent
+  }
+}
+#endif
+
+#if defined(SINGLE_PRECISION)
+static __inline fp_t __compiler_rt_logbf(fp_t x) {
+  return __compiler_rt_logbX(x);
+}
+#elif defined(DOUBLE_PRECISION)
+static __inline fp_t __compiler_rt_logb(fp_t x) {
+  return __compiler_rt_logbX(x);
+}
+#elif defined(QUAD_PRECISION)
+  #if defined(CRT_LDBL_128BIT)
+static __inline fp_t __compiler_rt_logbl(fp_t x) {
+  return __compiler_rt_logbX(x);
+}
+  #else
+// Rely on the libm implementation, for now.
+static __inline long double __compiler_rt_logbl(long double x) {
+  return crt_logbl(x);
+}
+  #endif
 #endif
 
 #endif // FP_LIB_HEADER

lib/builtins/int_math.h

 #define crt_fmaxl(x, y) __max((x), (y))
 #else
 #define crt_fmax(x, y) __builtin_fmax((x), (y))
 #define crt_fmaxf(x, y) __builtin_fmaxf((x), (y))
 #define crt_fmaxl(x, y) __builtin_fmaxl((x), (y))
 #endif
 
 #if defined(_MSC_VER) && !defined(__clang__)
-#define crt_logb(x) logb((x))
-#define crt_logbf(x) logbf((x))
 #define crt_logbl(x) logbl((x))
 #else
-#define crt_logb(x) __builtin_logb((x))
-#define crt_logbf(x) __builtin_logbf((x))
 #define crt_logbl(x) __builtin_logbl((x))
 #endif
 
 #if defined(_MSC_VER) && !defined(__clang__)
 #define crt_scalbn(x, y) scalbn((x), (y))
 #define crt_scalbnf(x, y) scalbnf((x), (y))
 #define crt_scalbnl(x, y) scalbnl((x), (y))
 #else
 #define crt_scalbn(x, y) __builtin_scalbn((x), (y))
 #define crt_scalbnf(x, y) __builtin_scalbnf((x), (y))
 #define crt_scalbnl(x, y) __builtin_scalbnl((x), (y))
 #endif
 
 #endif /* INT_MATH_H */

lib/builtins/ppc/divtc3.c

 /* This file is distributed under the University of Illinois Open Source
  * License. See LICENSE.TXT for details.
  */
 
 #include "DD.h"
 #include "../int_math.h"
+// Use DOUBLE_PRECISION because the soft-fp method we use is logb (on the upper
+// half of the long doubles), even though this file defines complex division for
+// 128-bit floats.
+#define DOUBLE_PRECISION
+#include "../fp_lib.h"
 
 #if !defined(CRT_INFINITY) && defined(HUGE_VAL)
 #define CRT_INFINITY HUGE_VAL
 #endif /* CRT_INFINITY */
 
 #define makeFinite(x) { \
     (x).s.hi = crt_copysign(crt_isinf((x).s.hi) ? 1.0 : 0.0, (x).s.hi); \
     (x).s.lo = 0.0;                                                     \
   }
 
 long double _Complex
 __divtc3(long double a, long double b, long double c, long double d)
 {
   DD cDD = { .ld = c };
   DD dDD = { .ld = d };
   
   int ilogbw = 0;
-  const double logbw = crt_logb(crt_fmax(crt_fabs(cDD.s.hi), crt_fabs(dDD.s.hi) ));
-  
-  if (crt_isfinite(logbw))
+  const double logbw = __compiler_rt_logb(
+    crt_fmax(crt_fabs(cDD.s.hi), crt_fabs(dDD.s.hi)));
+        if (crt_isfinite(logbw))
   {
     ilogbw = (int)logbw;
     
     cDD.s.hi = crt_scalbn(cDD.s.hi, -ilogbw);
     cDD.s.lo = crt_scalbn(cDD.s.lo, -ilogbw);
     dDD.s.hi = crt_scalbn(dDD.s.hi, -ilogbw);
     dDD.s.lo = crt_scalbn(dDD.s.lo, -ilogbw);
   }

test/builtins/Unit/compiler_rt_logb_test.c

+// RUN: %clang_builtins %s %librt -lm -o %t && %run %t
+//===-- compiler_rt_logb_test.c - Test __compiler_rt_logb -----------------===//
+//
+//                     The LLVM Compiler Infrastructure
+//
+// This file is dual licensed under the MIT and the University of Illinois Open
+// Source Licenses. See LICENSE.TXT for details.
+//
+//===----------------------------------------------------------------------===//
+//
+// This file checks __compiler_rt_logb from the compiler_rt library for
+// conformance against libm.
+//
+//===----------------------------------------------------------------------===//
+
+#define DOUBLE_PRECISION
+#include <math.h>
+#include <stdio.h>
+#include "fp_lib.h"
+
+int test__compiler_rt_logb(fp_t x) {
+  fp_t crt_value = __compiler_rt_logb(x);
+  fp_t libm_value = logb(x);
+  // Compare actual rep, e.g. to avoid NaN != the same NaN
+  if (toRep(crt_value) != toRep(libm_value)) {
+    printf("error: in __compiler_rt_logb(%a [%lX]) = %a [%lX] !=  %a [%lX]\n",
+           x, toRep(x), crt_value, toRep(crt_value), libm_value,
+           toRep(libm_value));
+    return 1;
+  }
+  return 0;
+}
+
+double cases[] = {
+    1.e-6, -1.e-6, NAN, -NAN, INFINITY, -INFINITY, -1,
+    -0.0,  0.0,    1,   -2,   2,        -0.5,      0.5,
+};
+
+int main() {
+  const unsigned N = sizeof(cases) / sizeof(cases[0]);
+  unsigned i;
+  for (i = 0; i < N; ++i) {
+    if (test__compiler_rt_logb(cases[i])) return 1;
+  }
+
+  // Test a moving 1 bit, especially to handle denormal values.
+  // Test the negation as well.
+  rep_t x = signBit;
+  while (x) {
+    if (test__compiler_rt_logb(fromRep(x))) return 1;
+    if (test__compiler_rt_logb(fromRep(signBit ^ x))) return 1;
+    x >>= 1;
+  }
+  // Also try a couple moving ones
+  x = signBit | (signBit >> 1) | (signBit >> 2);
+  while (x) {
+    if (test__compiler_rt_logb(fromRep(x))) return 1;
+    if (test__compiler_rt_logb(fromRep(signBit ^ x))) return 1;
+    x >>= 1;
+  }
+
+  return 0;
+}

test/builtins/Unit/compiler_rt_logbf_test.c

+// RUN: %clang_builtins %s %librt -lm -o %t && %run %t
+//===-- compiler_rt_logbf_test.c - Test __compiler_rt_logbf ---------------===//
+//
+//                     The LLVM Compiler Infrastructure
+//
+// This file is dual licensed under the MIT and the University of Illinois Open
+// Source Licenses. See LICENSE.TXT for details.
+//
+//===----------------------------------------------------------------------===//
+//
+// This file checks __compiler_rt_logbf from the compiler_rt library for
+// conformance against libm.
+//
+//===----------------------------------------------------------------------===//
+
+#define SINGLE_PRECISION
+#include <math.h>
+#include <stdio.h>
+#include "fp_lib.h"
+
+int test__compiler_rt_logbf(fp_t x) {
+  fp_t crt_value = __compiler_rt_logbf(x);
+  fp_t libm_value = logbf(x);
+  // Compare actual rep, e.g. to avoid NaN != the same NaN
+  if (toRep(crt_value) != toRep(libm_value)) {
+    printf("error: in __compiler_rt_logb(%a [%X]) = %a [%X] !=  %a [%X]\n", x,
+           toRep(x), crt_value, toRep(crt_value), libm_value,
+           toRep(libm_value));
+    return 1;
+  }
+  return 0;
+}
+
+double cases[] = {
+    1.e-6, -1.e-6, NAN, -NAN, INFINITY, -INFINITY, -1,
+    -0.0,  0.0,    1,   -2,   2,        -0.5,      0.5,
+};
+
+int main() {
+  const unsigned N = sizeof(cases) / sizeof(cases[0]);
+  unsigned i;
+  for (i = 0; i < N; ++i) {
+    if (test__compiler_rt_logbf(cases[i])) return 1;
+  }
+
+  // Test a moving 1 bit, especially to handle denormal values.
+  // Test the negation as well.
+  rep_t x = signBit;
+  while (x) {
+    if (test__compiler_rt_logbf(fromRep(x))) return 1;
+    if (test__compiler_rt_logbf(fromRep(signBit ^ x))) return 1;
+    x >>= 1;
+  }
+  // Also try a couple moving ones
+  x = signBit | (signBit >> 1) | (signBit >> 2);
+  while (x) {
+    if (test__compiler_rt_logbf(fromRep(x))) return 1;
+    if (test__compiler_rt_logbf(fromRep(signBit ^ x))) return 1;
+    x >>= 1;
+  }
+
+  return 0;
+}

test/builtins/Unit/compiler_rt_logbl_test.c

+// RUN: %clang_builtins %s %librt -lm -o %t && %run %t
+//===-- compiler_rt_logbl_test.c - Test __compiler_rt_logbl ---------------===//
+//
+//                     The LLVM Compiler Infrastructure
+//
+// This file is dual licensed under the MIT and the University of Illinois Open
+// Source Licenses. See LICENSE.TXT for details.
+//
+//===----------------------------------------------------------------------===//
+//
+// This file checks __compiler_rt_logbl from the compiler_rt library for
+// conformance against libm.
+//
+//===----------------------------------------------------------------------===//
+
+#define QUAD_PRECISION
+#include <math.h>
+#include <stdio.h>
+#include "fp_lib.h"
+#include "int_lib.h"
+
+#if defined(CRT_HAS_128BIT) && defined(CRT_LDBL_128BIT)
+
+int test__compiler_rt_logbl(fp_t x) {
+  fp_t crt_value = __compiler_rt_logbl(x);
+  fp_t libm_value = logbl(x);
+  // Compare actual rep, e.g. to avoid NaN != the same NaN
+  if (toRep(crt_value) != toRep(libm_value)) {
+    // Split expected values into two for printf
+    twords x_t, crt_value_t, libm_value_t;
+    x_t.all = toRep(x);
+    crt_value_t.all = toRep(crt_value);
+    libm_value_t.all = toRep(libm_value);
+    printf(
+        "error: in __compiler_rt_logb(%a [%llX %llX]) = %a [%llX %llX] !=  %a "
+        "[%llX %llX]\n",
+        x, x_t.s.high, x_t.s.low, crt_value, crt_value_t.s.high,
+        crt_value_t.s.low, libm_value, libm_value_t.s.high, libm_value_t.s.low);
+    return 1;
+  }
+  return 0;
+}
+
+double cases[] = {
+    1.e-6, -1.e-6, NAN, -NAN, INFINITY, -INFINITY, -1,
+    -0.0,  0.0,    1,   -2,   2,        -0.5,      0.5,
+};
+
+#endif
+
+int main() {
+#if defined(CRT_HAS_128BIT) && defined(CRT_LDBL_128BIT)
+  const unsigned N = sizeof(cases) / sizeof(cases[0]);
+  unsigned i;
+  for (i = 0; i < N; ++i) {
+    if (test__compiler_rt_logbl(cases[i])) return 1;
+  }
+
+  // Test a moving 1 bit, especially to handle denormal values.
+  // Test the negation as well.
+  rep_t x = signBit;
+  while (x) {
+    if (test__compiler_rt_logbl(fromRep(x))) return 1;
+    if (test__compiler_rt_logbl(fromRep(signBit ^ x))) return 1;
+    x >>= 1;
+  }
+  // Also try a couple moving ones
+  x = signBit | (signBit >> 1) | (signBit >> 2);
+  while (x) {
+    if (test__compiler_rt_logbl(fromRep(x))) return 1;
+    if (test__compiler_rt_logbl(fromRep(signBit ^ x))) return 1;
+    x >>= 1;
+  }
+#else
+  printf("skipped\n");
+#endif
+
+  return 0;
+}