[compiler-rt][builtins] Scaled Integer log10()

[leonardchan]

Summary

This patch proposes adding a log10() function that works on scaled integers. The only 2 arguments are the integer itself and its scale. The result is also a scaled integer with the same scale as the input. If the true result cannot be precisely represented, it is rounded down towards negative infinity.

I think this would be a good function to complement fixed point types, although we use scaled integers here since the types are still hidden behind a flag. This function also isn't something specified in the Embedded C spec, though we would like to use it in Fuchsia with fixed point types.

[leonardchan]

*ping*

[bjope]

I don't know much about compiler-rt, and neither what the procedure is for adding new builtins.
Maybe try to find some more reviewers (someone that has been doing/reviewing this kind of additions before).

[rjmccall]

I don't review compiler-rt patches generally. Can you explain what would actually use this function? Is there supposed to be a new Clang builtin for it?

[saugustine]

I can't speak to how useful this would be and for whom, but I'm not opposed to including it if Fuschia has a use for it.

Typically, compiler-support libraryfunctions are all lowercase, and just use a single letter for each type--"i" for integer--prefixed by a "size of type" character: 's' for "single", 'd' for double length as in this example (it's a little archaic), I don't know of a common one for "scaled" though.

Then each one is also numbered. The exact mechanism to choose the number escapes me at the moment, but I'm consulting with someone who should know.

So this should probably be named something more like __log10scaleddiX, while I look for what X should be.

Also, has this been run through clang-format? I can never remember the rules on line-length.

[saugustine]

The number is the argument count plus one for return value. So this function should be named

__log10scaleddi3

[leonardchan]

In D62088#1522446, @rjmccall wrote:

I don't review compiler-rt patches generally. Can you explain what would actually use this function? Is there supposed to be a new Clang builtin for it?

We will be using this as part of our WLAN library which represents energy units (decibels and watts) as scaled integers that we have custom classes for. Converting between these units and adding decibels mathematically involves using log10(). Currently, we just use float log10 and casting, but would like to use this for scaled integers and eventually fixed point types. The idea is that this would be a new clang builtin that we could use instead of a floating point log10().

[rjmccall]

In D62088#1531309, @leonardchan wrote:

In D62088#1522446, @rjmccall wrote:

I don't review compiler-rt patches generally. Can you explain what would actually use this function? Is there supposed to be a new Clang builtin for it?

We will be using this as part of our WLAN library which represents energy units (decibels and watts) as scaled integers that we have custom classes for. Converting between these units and adding decibels mathematically involves using log10(). Currently, we just use float log10 and casting, but would like to use this for scaled integers and eventually fixed point types. The idea is that this would be a new clang builtin that we could use instead of a floating point log10().

Okay. Will you need a signed variant of this? Is using a 64-bit argument reasonable for all targets, or should there be 32-bit and 64-bit (and maybe eventually 128-bit) variants?

[leonardchan]

In D62088#1523994, @saugustine wrote:

The number is the argument count plus one for return value. So this function should be named

__log10scaleddi3

Done and clang-formated

In D62088#1531310, @rjmccall wrote:

In D62088#1531309, @leonardchan wrote:

In D62088#1522446, @rjmccall wrote:

I don't review compiler-rt patches generally. Can you explain what would actually use this function? Is there supposed to be a new Clang builtin for it?

We will be using this as part of our WLAN library which represents energy units (decibels and watts) as scaled integers that we have custom classes for. Converting between these units and adding decibels mathematically involves using log10(). Currently, we just use float log10 and casting, but would like to use this for scaled integers and eventually fixed point types. The idea is that this would be a new clang builtin that we could use instead of a floating point log10().

Okay. Will you need a signed variant of this? Is using a 64-bit argument reasonable for all targets, or should there be 32-bit and 64-bit (and maybe eventually 128-bit) variants?

I don't think we need a signed version since ideally all inputs to log functions should be non-negative, otherwise the result is imaginary. This could also be reasonably extended to support other sizes. I only chose a 64 bit size for now since it's the largest type we use and any other types could be resized/rescaled into this. For 128-bit, this would be a little more tricky since the algorithm depends on wide multiplication (which we kinda do through x = (__uint128_t)x * x >> scale;), although I can't seem to find an existing compiler-rt multiplication function that does this.

[rjmccall]

In D62088#1531324, @leonardchan wrote:

In D62088#1523994, @saugustine wrote:

The number is the argument count plus one for return value. So this function should be named

__log10scaleddi3

Done and clang-formated

In D62088#1531310, @rjmccall wrote:

In D62088#1531309, @leonardchan wrote:

In D62088#1522446, @rjmccall wrote:

I don't review compiler-rt patches generally. Can you explain what would actually use this function? Is there supposed to be a new Clang builtin for it?

We will be using this as part of our WLAN library which represents energy units (decibels and watts) as scaled integers that we have custom classes for. Converting between these units and adding decibels mathematically involves using log10(). Currently, we just use float log10 and casting, but would like to use this for scaled integers and eventually fixed point types. The idea is that this would be a new clang builtin that we could use instead of a floating point log10().

Okay. Will you need a signed variant of this? Is using a 64-bit argument reasonable for all targets, or should there be 32-bit and 64-bit (and maybe eventually 128-bit) variants?

I don't think we need a signed version since ideally all inputs to log functions should be non-negative, otherwise the result is imaginary.

Well, I think the builtin has to be able to operate on signed scaled-integer / fixed-point types, and you'll have to specify the behavior for negative numbers. If you're comfortable implementing those semantics inline in the caller, then I agree you don't need signed functions in compiler-rt.

This could also be reasonably extended to support other sizes. I only chose a 64 bit size for now since it's the largest type we use and any other types could be resized/rescaled into this. For 128-bit, this would be a little more tricky since the algorithm depends on wide multiplication (which we kinda do through x = (__uint128_t)x * x >> scale;), although I can't seem to find an existing compiler-rt multiplication function that does this.

My concern about that is that it's going to limit this builtin to 64-bit targets — most 32-bit targets don't support __uint128_t (as an implementation concern), and there might be problems with using a 64-bit type in the parameter list.

As a separate concern, 128-bit multiplication might itself be implemented with a compiler-rt function on some targets, and I'm not sure how acceptable it is for compiler-rt functions to have this kind of cross-dependency.

[leonardchan]

Well, I think the builtin has to be able to operate on signed scaled-integer / fixed-point types, and you'll have to specify the behavior for negative numbers. If you're comfortable implementing those semantics inline in the caller, then I agree you don't need signed functions in compiler-rt.

Ok. I'll add it to the caller when adding the builtin clang function. We can handle usage with different types there.

This could also be reasonably extended to support other sizes. I only chose a 64 bit size for now since it's the largest type we use and any other types could be resized/rescaled into this. For 128-bit, this would be a little more tricky since the algorithm depends on wide multiplication (which we kinda do through x = (__uint128_t)x * x >> scale;), although I can't seem to find an existing compiler-rt multiplication function that does this.

My concern about that is that it's going to limit this builtin to 64-bit targets — most 32-bit targets don't support __uint128_t (as an implementation concern), and there might be problems with using a 64-bit type in the parameter list.

As a separate concern, 128-bit multiplication might itself be implemented with a compiler-rt function on some targets, and I'm not sure how acceptable it is for compiler-rt functions to have this kind of cross-dependency.

I updated this patch and made a separate 32 bit version for this. The 64 bit version now uses __uint128_t only if compiler rt supports it, otherwise it uses a wide multiplication that I borrowed from fp_lib.hand a custom funnel shift right.

[rjmccall]

In D62088#1531773, @leonardchan wrote:

Well, I think the builtin has to be able to operate on signed scaled-integer / fixed-point types, and you'll have to specify the behavior for negative numbers. If you're comfortable implementing those semantics inline in the caller, then I agree you don't need signed functions in compiler-rt.

Ok. I'll add it to the caller when adding the builtin clang function. We can handle usage with different types there.

Okay. Please mention the behavior for 0 in the documentation, since this function is not otherwise defined for 0.

This could also be reasonably extended to support other sizes. I only chose a 64 bit size for now since it's the largest type we use and any other types could be resized/rescaled into this. For 128-bit, this would be a little more tricky since the algorithm depends on wide multiplication (which we kinda do through x = (__uint128_t)x * x >> scale;), although I can't seem to find an existing compiler-rt multiplication function that does this.

My concern about that is that it's going to limit this builtin to 64-bit targets — most 32-bit targets don't support __uint128_t (as an implementation concern), and there might be problems with using a 64-bit type in the parameter list.

As a separate concern, 128-bit multiplication might itself be implemented with a compiler-rt function on some targets, and I'm not sure how acceptable it is for compiler-rt functions to have this kind of cross-dependency.

I updated this patch and made a separate 32 bit version for this. The 64 bit version now uses __uint128_t only if compiler rt supports it, otherwise it uses a wide multiplication that I borrowed from fp_lib.hand a custom funnel shift right.

Thanks. That sounds right to me, but I still don't feel qualified to actually review the approach. :)

[leonardchan]

In D62088#1537081, @rjmccall wrote:

In D62088#1531773, @leonardchan wrote:

Well, I think the builtin has to be able to operate on signed scaled-integer / fixed-point types, and you'll have to specify the behavior for negative numbers. If you're comfortable implementing those semantics inline in the caller, then I agree you don't need signed functions in compiler-rt.

Ok. I'll add it to the caller when adding the builtin clang function. We can handle usage with different types there.

Okay. Please mention the behavior for 0 in the documentation, since this function is not otherwise defined for 0.

Done. We just return INT64/32_MIN for input of 0.

This could also be reasonably extended to support other sizes. I only chose a 64 bit size for now since it's the largest type we use and any other types could be resized/rescaled into this. For 128-bit, this would be a little more tricky since the algorithm depends on wide multiplication (which we kinda do through x = (__uint128_t)x * x >> scale;), although I can't seem to find an existing compiler-rt multiplication function that does this.

My concern about that is that it's going to limit this builtin to 64-bit targets — most 32-bit targets don't support __uint128_t (as an implementation concern), and there might be problems with using a 64-bit type in the parameter list.

As a separate concern, 128-bit multiplication might itself be implemented with a compiler-rt function on some targets, and I'm not sure how acceptable it is for compiler-rt functions to have this kind of cross-dependency.

I updated this patch and made a separate 32 bit version for this. The 64 bit version now uses __uint128_t only if compiler rt supports it, otherwise it uses a wide multiplication that I borrowed from fp_lib.hand a custom funnel shift right.

Thanks. That sounds right to me, but I still don't feel qualified to actually review the approach. :)

No problem. @saugustine , would you feel comfortable reviewing this, or is there someone else you would recommend for compiler-rt?

[saugustine]

I am absolutely not the right person to review the math. I think we need to trust that the author has it correct. The logic looks reasonable and the tests look basically thorough. I am inclined to accept it.

[bjope]

Please add test cases for scale=0 and scale=width as I assume those need special handling (UB right now?).
And if scale=0 and scale=width needs special handling, then I guess scale=1 and scale=width-1 are new boundary values so I maybe it would be nice to have tests for those scales as well.

compiler-rt/lib/builtins/log10scaleddi3.c

65

So with scale=0 being a valid input this results in a negative shift count (UB).

67

So with scale=64 being a valid input this results in a too big shift count (UB).

[leonardchan]

In D62088#1537470, @bjope wrote:

Please add test cases for scale=0 and scale=width as I assume those need special handling (UB right now?).
And if scale=0 and scale=width needs special handling, then I guess scale=1 and scale=width-1 are new boundary values so I maybe it would be nice to have tests for those scales as well.

Added. For a scale of 0 we only return the integer result. For widths > 60 or 28, we return INT_MAX to represent an error since we need to represent 10 << scale in 64/32 bits for this to work. So the scale boundaries are [0,60] for the 64 bit version and [0,28] for the 32 bit version. We could technically increase this to go up to 63/31, although to get the precise result would require a bigger buffer (using more 128 bit ints) which would complicate things further.

For the 64 bit function, I realized that in order to get the precise result without a 128 bit int would require implementing division by 10 in 2 64 bit ints. I eventually found a way to do this and get the precise result, but am not sure if it should be included in this patch since it would be adding another large function. Perhaps it could be readded in a followup? For now, I removed the functions we fallback to if 128 bit ints aren't supported and only define the function if 128 bits are enabled.

[rjmccall]

In D62088#1549139, @leonardchan wrote:

In D62088#1537470, @bjope wrote:

Please add test cases for scale=0 and scale=width as I assume those need special handling (UB right now?).
And if scale=0 and scale=width needs special handling, then I guess scale=1 and scale=width-1 are new boundary values so I maybe it would be nice to have tests for those scales as well.

Added. For a scale of 0 we only return the integer result. For widths > 60 or 28, we return INT_MAX to represent an error since we need to represent 10 << scale in 64/32 bits for this to work. So the scale boundaries are [0,60] for the 64 bit version and [0,28] for the 32 bit version. We could technically increase this to go up to 63/31, although to get the precise result would require a bigger buffer (using more 128 bit ints) which would complicate things further.

For the 64 bit function, I realized that in order to get the precise result without a 128 bit int would require implementing division by 10 in 2 64 bit ints. I eventually found a way to do this and get the precise result, but am not sure if it should be included in this patch since it would be adding another large function. Perhaps it could be readded in a followup? For now, I removed the functions we fallback to if 128 bit ints aren't supported and only define the function if 128 bits are enabled.

Is there not a function for that already in compiler-rt? Or is the problem that it doesn't necessarily exist, e.g. if we're supporting a 64-bit scaled integer type but not a 128-bit integer type?

At any rate, I'd say it's better to have the operation with a big code-size cost than to not reliably have it.

[leonardchan]

Is there not a function for that already in compiler-rt? Or is the problem that it doesn't necessarily exist, e.g. if we're supporting a 64-bit scaled integer type but not a 128-bit integer type?

Yes, the extra functions are for if 64 bit scaled ints are supported but we aren't able to use 128 bit ints. There is __udivmodti4 which does 128 bit division, but only if 128 bit ints are enabled.

At any rate, I'd say it's better to have the operation with a big code-size cost than to not reliably have it.

Added.

[rjmccall]

Thanks! Still don't feel qualified to actually approve compiler-rt changes, though.

[leonardchan]

Hmm is there usually a protocol/etiquette for getting a new built-in submitted to compiler-rt?

Perhaps I could email the original creator and request a rubber stamp LGTM if he approves?

(my email reply might not have been sent)

[rjmccall]

My understanding is that Tim Northover sometimes reviews patches to compiler-rt, even if he's not an official owner; CC'ing him.

[leonardchan]

@t.p.northover Would you feel ok reviewing this?

[leonardchan]

Been a while since the last update, but would anyone still be comfortable reviewing this?

compiler-rt/lib/builtins/CMakeLists.txt

   floatuntidf.c
   floatuntisf.c
   int_util.c
+  log10scaledsi3.c
+  log10scaleddi3.c
   lshrdi3.c
   lshrti3.c
   moddi3.c

compiler-rt/lib/builtins/log10scaleddi3.c

+//===------- lib/log10scaleddi3.c - Scaled Integer log10() --------*- 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
+//
+//===----------------------------------------------------------------------===//
+//
+// This file implements log10() for 64 bit scaled integers, with a scale of at
+// most 60. This implementation is based on Clay. S. Turner's fast binary
+// logarithm algorithm. If the true result of log10() cannot be precisely stored
+// in the result type, the value is rounded down towards negative infinity.
+// Only up to 60 bits are allowed for the scale since 4 bits are required to
+// reppresent 10.
+//
+//===----------------------------------------------------------------------===//
+
+#include <assert.h>
+#include <stddef.h>
+
+#include "int_lib.h"
+
+#define BASE 10
+
+#if !defined(__SIZEOF_INT128__)
+#define loWord(a) (a & 0xffffffffU)
+#define hiWord(a) (a >> 32)
+
+// 64x64 -> 128 wide multiply for platforms that don't have such an operation.
+static void wideMultiply(uint64_t a, uint64_t b, uint64_t *hi, uint64_t *lo) {
+  // Each of the component 32x32 -> 64 products
+  const uint64_t plolo = loWord(a) * loWord(b);
+  const uint64_t plohi = loWord(a) * hiWord(b);
+  const uint64_t philo = hiWord(a) * loWord(b);
+  const uint64_t phihi = hiWord(a) * hiWord(b);
+  // Sum terms that contribute to lo in a way that allows us to get the carry
+  const uint64_t r0 = loWord(plolo);
+  const uint64_t r1 = hiWord(plolo) + loWord(plohi) + loWord(philo);
+  *lo = r0 + (r1 << 32);
+  // Sum terms contributing to hi with the carry from lo
+  *hi = hiWord(plohi) + hiWord(philo) + hiWord(r1) + phihi;
+}
+
+// 128 bit funnel shift right from a 64 bit hi and lo.
+static void funnelShiftRight(uint64_t hi, uint64_t lo, uint32_t shift,
+                             uint64_t *res_hi, uint64_t *res_lo) {
+  if (shift >= 128) {
+    *res_hi = 0;
+    *res_lo = 0;
+  } else if (shift > 64) {
+    *res_hi = 0;
+    *res_lo = hi >> (shift - 64);
+  } else if (shift == 64) {
+    *res_hi = 0;
+    *res_lo = hi;
+  } else if (shift > 0) {
+    *res_hi = hi >> shift;
+    *res_lo = (hi << (64 - shift)) | (lo >> shift);
+  } else {
+    *res_hi = hi;
+    *res_lo = lo;
+  }
+}
+
+#define NUM_BITS 64
+
+void div10(uint64_t hi, uint64_t lo, uint64_t *res_hi, uint64_t *res_lo) {
+  uint64_t quot = 0, rem = 0;
+
+  for (unsigned i = 0; i < NUM_BITS; ++i) {
+    quot <<= 1;
+    rem <<= 1;
+    rem |= (hi >> ((NUM_BITS - 1) - i)) & 1;
+
+    if (rem >= BASE) {
+      rem -= BASE;
+      quot |= 1;
+    }
+  }
+
+  *res_hi = quot;
+
+  for (unsigned i = 0; i < NUM_BITS; ++i) {
+    quot <<= 1;
+    rem <<= 1;
+    rem |= (lo >> ((NUM_BITS - 1) - i)) & 1;
+
+    if (rem >= BASE) {
+      rem -= BASE;
+      quot |= 1;
+    }
+  }
+
+  *res_lo = quot;
+}
+
+#endif
+
+// Returns: log10(x), rounded down towards negative infinity in the scale
+//          provided in the 2nd argument. If x is zero, INT64_MIN is returned.
+//          If the true result cannot fit inside a signed 64 bit integer,
+//          INT64_MIN is returned. This only affects large, negative results.
+//          If the scale is greater than 60, INT64_MAX is returned to respresent
+//          an invalid scale.
+//
+// Assumptions: x is represents an unsigned scaled integer. When the scale is 0,
+//              this becomes a regular integer log10() function. The scale must
+//              be at most 60 since we will need to respresent a scaled 10 (10
+//              << scale) in an unsigned 64 bit integer.
+
+int64_t __log10scaleddi3(uint64_t x, uint32_t scale) {
+  // An input of zero will return the smallest value possible to represent
+  // negative infinity.
+  if (x == 0)
+    return INT64_MIN;
+
+  // We must be able to represent a scaled 10 in the bit width we have now.
+  // In order to fit a scaled value of 10 in 64 bits, we need 4 bits for the
+  // integer, leaving 60 bits for the scale.
+  if (scale > 60)
+    return INT64_MAX;
+
+  uint64_t oneval = UINT64_C(1) << scale;
+  uint64_t baseval = (uint64_t)BASE << scale;
+  int64_t y = 0;
+
+  // Normalize x up to at least 1.
+  while (x < oneval) {
+    x *= BASE;
+    --y;
+  }
+
+  // Normalize x to at most 10.
+  while (x >= baseval) {
+    x /= BASE;
+    ++y;
+  }
+
+  // Return early since there is not fractional part.
+  if (scale == 0)
+    return y;
+
+  // We overflowed past lowest value we can represent.
+  if (y < (INT64_MIN >> scale))
+    return INT64_MIN;
+
+  y <<= scale;
+
+  // At this point, the value of res is between a scaled 1 and 10.
+  // Here we calculate the fractional part of the result.
+  int64_t b = INT64_C(1) << (scale - 1);
+
+#if defined(__SIZEOF_INT128__)
+  __uint128_t res = x;
+  for (size_t i = 0; i < scale; i++) {
+    assert(res <= UINT64_MAX);
+    res = res * res >> (scale - 1);
+
+    // In order to remain as close as possible to the correct solution, we
+    // attempt to round using the bit before the LSB.
+    int roundbit = res & 1;
+    res >>= 1;
+    if (roundbit)
+      ++res;
+
+    if (res >= baseval) {
+      res /= BASE;
+      y += b;
+    }
+    b >>= 1;
+  }
+#else
+  uint64_t res_hi = 0, res_lo = x;
+  for (size_t i = 0; i < scale; i++) {
+    assert(res_hi == 0);
+    wideMultiply(res_lo, res_lo, &res_hi, &res_lo);
+    funnelShiftRight(res_hi, res_lo, scale - 1, &res_hi, &res_lo);
+
+    // In order to remain as close as possible to the correct solution, we
+    // attempt to round using the bit before the LSB.
+    int roundbit = res_lo & 1;
+    funnelShiftRight(res_hi, res_lo, 1, &res_hi, &res_lo);
+    if (roundbit) {
+      ++res_lo;
+      if (res_lo == 0) // overflow
+        ++res_hi;
+    }
+
+    if (res_hi > 0 || res_lo >= baseval) {
+      div10(res_hi, res_lo, &res_hi, &res_lo);
+      y += b;
+    }
+    b >>= 1;
+  }
+#endif
+
+  return y;
+}

compiler-rt/lib/builtins/log10scaledsi3.c

+//===------- lib/log10scaledsi3.c - Scaled Integer log10() --------*- 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
+//
+//===----------------------------------------------------------------------===//
+//
+// This file implements log10() for 32 bit scaled integers, with a scale of at
+// most 28. This implementation is based on Clay. S. Turner's fast binary
+// logarithm algorithm. If the true result of log10() cannot be precisely stored
+// in the result type, the value is rounded down towards negative infinity.
+// Only up to 28 bits are allowed for the scale since 4 bits are required to
+// represent 10.
+//
+//===----------------------------------------------------------------------===//
+
+#include <stddef.h>
+
+#include "int_lib.h"
+
+#define BASE 10
+
+// Returns: log10(x), rounded down towards negative infinity in the scale
+//          provided in the 2nd argument. If x is zero, INT32_MIN is returned.
+//          If the true result cannot fit inside a signed 32 bit integer,
+//          INT32_MIN is returned. This only affects large, negative results.
+//          If the scale is greater than 28, INT32_MAX is returned to respresent
+//          an invalid scale.
+//
+// Assumptions: x is represents an unsigned scaled integer. When the scale is 0,
+//              this becomes a regular integer log10() function. The scale must
+//              be at most 28 since we will need to respresent a scaled 10 (10
+//              << scale) in an unsigned 28 bit integer.
+
+int32_t __log10scaledsi3(uint32_t x, uint32_t scale) {
+  // An input of zero will return the smallest value possible to represent
+  // negative infinity.
+  if (x == 0)
+    return INT32_MIN;
+
+  // We must be able to represent a scaled 10 in the bit width we have now.
+  // In order to fit a scaled value of 10 in 64 bits, we need 4 bits for the
+  // integer, leaving 60 bits for the scale.
+  if (scale > 28)
+    return INT32_MAX;
+
+  uint32_t oneval = UINT32_C(1) << scale;
+  uint32_t baseval = (uint32_t)BASE << scale;
+  int32_t y = 0;
+
+  // Normalize x up to at least 1.
+  while (x < oneval) {
+    x *= BASE;
+    --y;
+  }
+
+  // Normalize x to at most 10.
+  while (x >= baseval) {
+    x /= BASE;
+    ++y;
+  }
+
+  // Return early since there is not fractional part.
+  if (scale == 0)
+    return y;
+
+  // We overflowed past lowest value we can represent.
+  if (y < (INT32_MIN >> scale))
+    return INT32_MIN;
+
+  y <<= scale;
+
+  // At this point, the value of res is between a scaled 1 and 10.
+  // Here we calculate the fractional part of the result.
+  int32_t b = INT32_C(1) << (scale - 1);
+
+  uint64_t res = x;
+  for (size_t i = 0; i < scale; i++) {
+    res = res * res >> (scale - 1);
+
+    // In order to remain as close as possible to the correct solution, we
+    // attempt to round using the bit before the LSB.
+    int roundbit = res & 1;
+    res >>= 1;
+    if (roundbit)
+      ++res;
+
+    if (res >= baseval) {
+      res /= BASE;
+      y += b;
+    }
+    b >>= 1;
+  }
+
+  return y;
+}

compiler-rt/test/builtins/Unit/log10scaledint.c

+// RUN: %clang_builtins %s %librt -o %t && %run %t
+
+#include <stdint.h>
+#include <stdio.h>
+
+#define CHECK_SUCCESS(x) \
+  if (x)                 \
+    return 1;
+
+#define SCALE 31
+#define SCALED_ONE (UINT64_C(1) << SCALE)
+#define __log10_64(x) __log10scaleddi3(x, SCALE)
+
+#define CHECK_LOG10_64_SCALE(VAL, SCALE, EXPECTED)                                      \
+  if (__log10scaleddi3(VAL, SCALE) != EXPECTED) {                                       \
+    printf(                                                                             \
+        "error in __log10scaleddi3(%llu, %u) on line %d. Expected %lld, found %lld.\n", \
+        VAL, SCALE, __LINE__, EXPECTED, __log10scaleddi3(VAL, SCALE));                  \
+    return 1;                                                                           \
+  }
+
+#define CHECK_LOG10_64(VAL, EXPECTED)                                                   \
+  if (__log10_64(VAL) != EXPECTED) {                                                    \
+    printf(                                                                             \
+        "error in __log10scaleddi3(%llu, %u) on line %d. Expected %lld, found %lld.\n", \
+        VAL, SCALE, __LINE__, EXPECTED, __log10_64(VAL));                               \
+    return 1;                                                                           \
+  }
+
+int64_t __log10scaleddi3(uint64_t x, uint32_t scale);
+int32_t __log10scaledsi3(uint32_t x, uint32_t scale);
+
+int check_log10_64() {
+  // Powers of 10.
+  // These values will not have the precision necessary to represent the true
+  // return values of 1og10() since reciprocal powers of 10 cannot be evenly
+  // represented with sums of recirpocal powers of 2 (1/10^n cannot be
+  // represented exactly with sums of 1/2^m). This is a binary log function so
+  // we can only attempt to represent powers of 10 as close as possible to sums
+  // of powers of 2.
+  // For example, we cannot get -3 exactly from log10(0.001), but we can
+  // precisely get log10((1 << SCALE) / 1000) and get a value close to -3 <<
+  // SCALE. The error between this close result and the true value of -3 <<
+  // SCALE will always be at most 1 if it is rounded.
+  CHECK_LOG10_64(SCALED_ONE / 10,
+                 -2147483652); // -2147483652 / 2^SCALE => -1.000...
+  CHECK_LOG10_64(SCALED_ONE / 100, -4294967317);
+  CHECK_LOG10_64(SCALED_ONE / 1000, -6442451226);
+  CHECK_LOG10_64(SCALED_ONE / 10000, -8589936177);
+
+  // These should produce precise whole integer results since the scaled integer
+  // input is not a fraction.
+  CHECK_LOG10_64(SCALED_ONE, 0);
+  CHECK_LOG10_64(10 * SCALED_ONE, SCALED_ONE);
+  CHECK_LOG10_64(100 * SCALED_ONE, 2 * SCALED_ONE);
+  CHECK_LOG10_64(1000 * SCALED_ONE, 3 * SCALED_ONE);
+  CHECK_LOG10_64(10000 * SCALED_ONE, 4 * SCALED_ONE);
+
+  // Sums of powers of 2
+  CHECK_LOG10_64(SCALED_ONE / 2 + SCALED_ONE / 4, -268303894);
+  CHECK_LOG10_64(SCALED_ONE / 2 + SCALED_ONE / 4 + SCALED_ONE / 8, -124536758);
+  CHECK_LOG10_64(SCALED_ONE / 2 + SCALED_ONE / 4 + SCALED_ONE / 8 +
+                     SCALED_ONE / 16,
+                 -60191226);
+  CHECK_LOG10_64(SCALED_ONE / 8 + SCALED_ONE / 16, -1561217881);
+
+  // Integers
+  CHECK_LOG10_64(3 * SCALED_ONE, 1024610092);
+  CHECK_LOG10_64(5 * SCALED_ONE, 1501026654);
+  CHECK_LOG10_64(7 * SCALED_ONE, 1814834221);
+  CHECK_LOG10_64(11 * SCALED_ONE, 2236373762);
+
+  CHECK_LOG10_64(3 * SCALED_ONE + SCALED_ONE / 2 + SCALED_ONE / 4, 1232722760);
+  CHECK_LOG10_64(5 * SCALED_ONE + SCALED_ONE / 2 + SCALED_ONE / 4 + SCALED_ONE / 8,
+                 1651431827);
+  CHECK_LOG10_64(7 * SCALED_ONE + SCALED_ONE / 2 + SCALED_ONE / 4 +
+                     SCALED_ONE / 8 + SCALED_ONE / 16,
+                 1932056116);
+  CHECK_LOG10_64(11 * SCALED_ONE + SCALED_ONE / 8 + SCALED_ONE / 16, 2252137072);
+
+  // Other numbers that can't evenly be broken into limited sums of powers of 2
+  CHECK_LOG10_64(6746518852, 1067621222); // pi
+  CHECK_LOG10_64(5837465777, 932640298);  // e
+  CHECK_LOG10_64(3037000499, 323228496);  // sqrt(2)
+
+  // Zero returns min val to represent -infinity
+  CHECK_LOG10_64(0, INT64_MIN);
+
+  // Max
+  CHECK_LOG10_64(UINT64_MAX, 21333080777);
+
+  // Min
+  CHECK_LOG10_64(1, -20040166791);
+  CHECK_LOG10_64(2, -19393709798);
+
+  // Scale of zero results in integr log10().
+  CHECK_LOG10_64_SCALE(0, 0, INT64_MIN);
+  CHECK_LOG10_64_SCALE(1, 0, 0);
+  CHECK_LOG10_64_SCALE(9, 0, 0);
+  CHECK_LOG10_64_SCALE(10, 0, 1);
+  CHECK_LOG10_64_SCALE(UINT64_MAX, 0, 19);
+
+  // Small scales. Large values.
+  CHECK_LOG10_64_SCALE(UINT64_MAX, 1, 37);
+  CHECK_LOG10_64_SCALE(UINT64_MAX, 2, 74);
+
+  // Large scales. Large values.
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 49, 50, -338929644074912);
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 54, 55, -10845748610397182);
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 57, 58, -86765988883177456);
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 58, 59, -173531977766354911);
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 59, 60, -347063955532709821);
+
+  // At this point, we are no longer able to fit a scaled 10 in 64 bits.
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 60, 61, INT64_MAX);
+  CHECK_LOG10_64_SCALE(1, 63, INT64_MAX);
+  CHECK_LOG10_64_SCALE(1, 64, INT64_MAX);
+  CHECK_LOG10_64_SCALE(0, 64, INT64_MIN); // Return min bc input is still zero.
+
+  // Large scales. Small values.
+  CHECK_LOG10_64_SCALE(1, 60, INT64_MIN); // Cannot represent 1 / 2^60 in 64 bits.
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 40, 60, -6941279110654196416);
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 34, 60, -9023662843850455340);
+  CHECK_LOG10_64_SCALE(UINT64_C(1) << 33, 60, INT64_MIN);
+
+  return 0;
+}
+
+#undef SCALE
+#undef SCALED_ONE
+
+#define SCALE 15
+#define SCALED_ONE (UINT32_C(1) << SCALE)
+#define __log10_32(x) __log10scaledsi3(x, SCALE)
+
+#define CHECK_LOG10_32_SCALE(VAL, SCALE, EXPECTED)                                \
+  if (__log10scaledsi3(VAL, SCALE) != EXPECTED) {                                 \
+    printf(                                                                       \
+        "error in __log10scaledsi3(%d, %u) on line %d. Expected %d, found %d.\n", \
+        VAL, SCALE, __LINE__, EXPECTED, __log10scaledsi3(VAL, SCALE));            \
+    return 1;                                                                     \
+  }
+
+#define CHECK_LOG10_32(VAL, EXPECTED)                                             \
+  if (__log10_32(VAL) != EXPECTED) {                                              \
+    printf(                                                                       \
+        "error in __log10scaledsi3(%d, %u) on line %d. Expected %d, found %d.\n", \
+        VAL, SCALE, __LINE__, EXPECTED, __log10_32(VAL));                         \
+    return 1;                                                                     \
+  }
+
+// Same as the previous test but for 32 bit ints.
+int check_log10_32() {
+  // Powers of 10.
+  CHECK_LOG10_32(SCALED_ONE / 10, -32772); // -32772 / 2^SCALE => -1.000...
+  CHECK_LOG10_32(SCALED_ONE / 100, -65566);
+  CHECK_LOG10_32(SCALED_ONE / 1000, -98642);
+  CHECK_LOG10_32(SCALED_ONE / 10000, -132328);
+
+  // These should produce precise whole integer results since the scaled integer
+  // input is not a fraction.
+  CHECK_LOG10_32(SCALED_ONE, 0);
+  CHECK_LOG10_32(10 * SCALED_ONE, SCALED_ONE);
+  CHECK_LOG10_32(100 * SCALED_ONE, 2 * SCALED_ONE);
+  CHECK_LOG10_32(1000 * SCALED_ONE, 3 * SCALED_ONE);
+  CHECK_LOG10_32(10000 * SCALED_ONE, 4 * SCALED_ONE);
+
+  // Sums of powers of 2
+  CHECK_LOG10_32(SCALED_ONE / 2 + SCALED_ONE / 4, -4094);
+  CHECK_LOG10_32(SCALED_ONE / 2 + SCALED_ONE / 4 + SCALED_ONE / 8, -1901);
+  CHECK_LOG10_32(SCALED_ONE / 2 + SCALED_ONE / 4 + SCALED_ONE / 8 +
+                     SCALED_ONE / 16,
+                 -919);
+  CHECK_LOG10_32(SCALED_ONE / 8 + SCALED_ONE / 16, -23823);
+
+  // Integers
+  CHECK_LOG10_32(3 * SCALED_ONE, 15634);
+  CHECK_LOG10_32(5 * SCALED_ONE, 22903);
+  CHECK_LOG10_32(7 * SCALED_ONE, 27692);
+  CHECK_LOG10_32(11 * SCALED_ONE, 34124);
+
+  CHECK_LOG10_32(3 * SCALED_ONE + SCALED_ONE / 2 + SCALED_ONE / 4, 18809);
+  CHECK_LOG10_32(5 * SCALED_ONE + SCALED_ONE / 2 + SCALED_ONE / 4 + SCALED_ONE / 8,
+                 25198);
+  CHECK_LOG10_32(7 * SCALED_ONE + SCALED_ONE / 2 + SCALED_ONE / 4 +
+                     SCALED_ONE / 8 + SCALED_ONE / 16,
+                 29480);
+  CHECK_LOG10_32(11 * SCALED_ONE + SCALED_ONE / 8 + SCALED_ONE / 16, 34364);
+
+  // Other numbers that can't evenly be broken into limited sums of powers of 2
+  CHECK_LOG10_32(102943, 16290); // pi
+  CHECK_LOG10_32(89072, 14230);  // e
+  CHECK_LOG10_32(46340, 4931);   // sqrt(2)
+
+  // Zero returns min val to represent -infinity
+  CHECK_LOG10_32(0, INT32_MIN);
+
+  // Max
+  CHECK_LOG10_32(UINT32_MAX, 167690);
+
+  // Min
+  CHECK_LOG10_32(1, -147963);
+  CHECK_LOG10_32(2, -138099);
+
+  // Scale of zero results in integr log10().
+  CHECK_LOG10_32_SCALE(0, 0, INT32_MIN);
+  CHECK_LOG10_32_SCALE(1, 0, 0);
+  CHECK_LOG10_32_SCALE(9, 0, 0);
+  CHECK_LOG10_32_SCALE(10, 0, 1);
+  CHECK_LOG10_32_SCALE(UINT32_MAX, 0, 9);
+
+  // Small scales. Large values.
+  CHECK_LOG10_32_SCALE(UINT32_MAX, 1, 18);
+  CHECK_LOG10_32_SCALE(UINT32_MAX, 2, 36);
+
+  // Large scales. Large values.
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 23, 24, -5050446);
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 24, 25, -10100891);
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 25, 26, -20201782);
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 26, 27, -40403563);
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 27, 28, -80807125);
+
+  // At this point, we are no longer able to fit a scaled 10 in 64 bits.
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 28, 29, INT32_MAX);
+  CHECK_LOG10_32_SCALE(1, 32, INT32_MAX);
+  CHECK_LOG10_32_SCALE(1, 64, INT32_MAX);
+  CHECK_LOG10_32_SCALE(0, 64, INT32_MIN); // Return min bc input is still zero.
+
+  // Large scales. Small values.
+  CHECK_LOG10_32_SCALE(1, 28, INT32_MIN); // Cannot represent 1 / 2^28 in 32 bits.
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 2, 28, -2100985229);
+  CHECK_LOG10_32_SCALE(UINT32_C(1) << 1, 28, INT32_MIN);
+
+  return 0;
+}
+
+int main() {
+  CHECK_SUCCESS(check_log10_64());
+  CHECK_SUCCESS(check_log10_32());
+  return 0;
+}