The Fastest (long) ilog10 in the West

#include <stdbit.h> #include <stdint.h>

static const uint64_t pow10_lz[65] = { 10000000000000000000ULL, 10000000000000000000ULL, 10000000000000000000ULL, 10000000000000000000ULL, 1000000000000000000, 1000000000000000000, 1000000000000000000, 100000000000000000, 100000000000000000, 100000000000000000, 10000000000000000, 10000000000000000, 10000000000000000, 1000000000000000, 1000000000000000, 1000000000000000, 100000000000000, 100000000000000, 100000000000000, 100000000000000, 10000000000000, 10000000000000, 10000000000000, 1000000000000, 1000000000000, 1000000000000, 100000000000, 100000000000, 100000000000, 10000000000, 10000000000, 10000000000, 1000000000, 1000000000, 1000000000, 1000000000, 100000000, 100000000, 100000000, 10000000, 10000000, 10000000, 1000000, 1000000, 1000000, 100000, 100000, 100000, 10000, 10000, 10000, 10000, 1000, 1000, 1000, 100, 100, 100, 10, 10, 10, 1, 1, 1, 0 };

uint64_t ilog10(uint64_t x) { uint64_t lz = stdc_leading_zeros(x); uint64_t e = (lz * 5) >> 4; return (19 - e) - (x < pow10_lz[lz]); }

But if you’re curious how I got here, and when this is the fastest, read on:

I was attempting to write a silly little C++ program, whose purpose I have since forgotten. As I often do when learning a new language, I imposed a set of restrictions on myself to develop my skills: this time I decided to try and use no branching, loops, recursion, or division.1

At some point in the program I needed to count the (base 10) digits in a number. After some quick research I learned that this is the same as ceil(log10(x + 1)), or more simply expressed ilog10(x) + 1. Quick and easy, no loops or branching, I thought. However, curiosity got the better of me and I decided to take a look at glibc’s implementation of these functions. Let’s look at ceil first:

double __ceil (double x) { #if USE_CEIL_BUILTIN return __builtin_ceil (x); #else /* Use generic implementation. / int64_t i0, i; int32_t j0; EXTRACT_WORDS64 (i0, x); j0 = ((i0 >> 52) & 0x7ff) - 0x3ff; if (j0 <= 51) { if (j0 < 0) { / return 0 * sign(x) if |x| < 1 / if (i0 < 0) i0 = INT64_C (0x8000000000000000); else if (i0 != 0) i0 = INT64_C (0x3ff0000000000000); } else { i = INT64_C (0x000fffffffffffff) >> j0; if ((i0 & i) == 0) return x; / x is integral / if (i0 > 0) i0 += UINT64_C (0x0010000000000000) >> j0; i0 &= ~i; } } else { if (j0 == 0x400) return x + x; / inf or NaN / else return x; / x is integral / } INSERT_WORDS64 (x, i0); return x; #endif / ! USE_CEIL_BUILTIN */ }

Copyright (C) 1993 Sun Microsystems, Inc. SunPro license Ignoring the very likely possibility of a native implementation, this is not ideal if we’re avoiding branching. However, in this instance I was not concerned with handling of special floating point cases, so instead of ceil we can use truncation when converting from double to long, by taking log10(x) + 1 rather than log10(x + 1). Now let’s look at log10:

double __ieee754_log10 (double x) { double y, z; int64_t i, hx; int32_t k;

EXTRACT_WORDS64 (hx, x);

k = 0; if (hx < INT64_C(0x0010000000000000)) { /* x < 2**-1022 / if (__glibc_unlikely ((hx & UINT64_C(0x7fffffffffffffff)) == 0)) return -two54 / fabs (x); / log(+-0)=-inf / if (__glibc_unlikely (hx < 0)) return (x - x) / (x - x); / log(-#) = NaN */ k -= 54; x = two54; / subnormal number, scale up x / EXTRACT_WORDS64 (hx, x); } / scale up resulted in a NaN number */ if (__glibc_unlikely (hx >= UINT64_C(0x7ff0000000000000))) return x + x; k += (hx >> 52) - 1023; i = ((uint64_t) k & UINT64_C(0x8000000000000000)) >> 63; hx = (hx & UINT64_C(0x000fffffffffffff)) | ((0x3ff - i) << 52); y = (double) (k + i); if (FIX_INT_FP_CONVERT_ZERO && y == 0.0) y = 0.0; INSERT_WORDS64 (x, hx); z = y * log10_2lo + ivln10 * __ieee754_log (x); return z + y * log10_2hi; }

Copyright (C) 1993 Sun Microsystems, Inc. SunPro license This is even worse: not only is there no native implementation, it’s still branching, and there’s no easy alternative to use in its place.

So I decided to write my own alternative. The first thing to consider is that counting the digits of an integer doesn’t require us to deal with any floating point special cases or any of the other baggage that comes with IEEE 754. From this point on I shifted my focus to the more widely applicable ilog10, as you can easily convert that back to counting digits with a + 1. So I started work on…

, which I came up with at the suggestion of a friend:

long ilog10(long x) { int n = -1; while (x != 0) { x /= 10; n++; } return n; }

step or run

While this still contains looping, it is far less complicated than either of the two functions above, and serves as a clean method of understanding what it is we’re trying to compute. The integer coarse logarithm base 10 is ilog⁡10(n)=y  ⟺  10y≤n<10y+1\operatorname{ilog}_{10}(n) = y \iff 10^{y} \le n < 10^{y+1}ilog10​(n)=y⟺10y≤n<10y+1 or in simpler terms, the exponent of the largest power of ten that is at or below the number.

Inspired by the fast inverse square root I decided to search for bitwise hacks to help me make this implementation more efficient. After much googling, and eventually poring over the manual for GCC, I found __builtin_clzl.

__builtin_clzl, or its 32-bit cousin __builtin_clz, is a GCC builtin function which takes an integer as its parameter and calculates the number of leading zeros in front of it. For example, were there an 8-bit version, the number of leading zeros returned by the function for 00001010 would be 4. The number of leading zeros is bounded by the datatype, e.g., the number of leading zeros for the number 1 would be 31 in a 32-bit datatype, and 63 in a 64-bit datatype, so on and so forth. On x86-64 this function becomes the opcode bsr, or, on supported architectures (Intel Haswell and later, AMD K10 and later), lzcnt2.

But you might be thinking, “How does this relate to calculating the ilog10 of a number?” Well, the astute among you may have already noticed you can take the length of the data type, subtract one, and then subtract the number of leading zeros, and get the logarithm base 2 of a number, i.e., ilog2(x) = TYPE_LENGTH - 1 - __builtin_clz{l}{l}(x).

Subtracting one is necessary because bits are counted from 0: in our 8-bit example 00001010 has 4 leading zeros, and 8 - 4 = 4 is how many bits it takes to write 10, but its highest set bit is bit 3, and 23=82^3 = 823=8 is the largest power of two not above 10, therefore ilog⁡2(10)=3\operatorname{ilog}_2(10) = 3ilog2​(10)=3..

At this point, however, I was stuck. While I had the log base 2 I didn’t know where to go next. This led me down several rabbit holes3 until I remembered the logarithm change of base formula:

Which in our case we would substitute 2 and 10 in for aaa and bbb respectively:

Which we can then simplify, for our purposes, to:

Meaning we can take our fast ilog2(x) and transform it into an ilog10(x) by multiplying it by the fraction 12334096\frac{1233}{4096}40961233​..

This, however, presents another problem as we have been dealing with integer types and the decimal version of 12334096\frac{1233}{4096}40961233​ is 0.301025390625 meaning floating point precision is required to perform this operation. Moving from floating point to integers and back is not ideal and claws back a lot of the gains we would have made doing this. However, there is a somewhat obvious way around this: by scaling up the constant (to the degree of precision we require4), multiplying by it, and then dividing by the same scale factor, we can produce a very similar result but with integers. So something like ilog⁡2(x)⋅301025390625/1012\operatorname{ilog}_2(x) \cdot 301025390625 / 10^{12}ilog2​(x)⋅301025390625/1012..

However, this too presents a problem as, if you recall from earlier, I was trying to avoid doing division. This can be solved by choosing our scale a bit more carefully. Instead of scaling by a power of ten we can scale by a power of two, as dividing by a power of two is just a right shift. As luck would have it our decimal version of the constant reduces down perfectly, as 1012=212⋅51210^{12} = 2^{12} \cdot 5^{12}1012=212⋅512 and 301025390625301025390625301025390625 can be rewritten as 1233⋅5121233 \cdot 5^{12}1233⋅512,, meaning:

So our division by 101210^{12}1012 becomes a division by 2122^{12}212,, i.e., a right shift by 12, leaving us with (ilog2(x) * 1233) >> 12, no division required.

Our code so far would be something like this5:

#include <limits.h> #include <stdint.h>

#define LONGSIZE ((int)(sizeof(int64_t) * CHAR_BIT))

uint64_t ilog10(uint64_t x) { uint64_t ilog2 = (LONGSIZE - 1) - __builtin_clzll(x); return (ilog2 * 1233) >> 12; }

step or run

Fixing the off-by-one

Sadly, this also has a problem, which is that this trick of doing ilog10(x) = (ilog2(x) * 1233) >> 12 is not precise enough for all values of x. For example: if we take x=10x = 10x=10 then ilog2(10) = 3, and (3 * 1233) >> 12 is 3699 >> 12, which is 0. However, ilog10(10) is of course 1.

The issue here is one of precision. The trick we are using does not preserve the same level of precision as the slower log⁡2x⋅12334096\log_2 x \cdot \frac{1233}{4096}log2​x⋅40961233​.. This results in soft edges around values slightly above each power of ten, where the lack of precision prevents the correct amount of accumulation from happening, preventing the output from incrementing at the right time in the series; or more simply this trick sometimes has an off-by-one error.

Unfortunately there are no apparent methods of doing this operation, while maintaining (or improving)6 the performance characteristics of the trick, that achieve the level of accuracy needed. Meaning in order for the faster way of doing the operation to work we will need to apply some kind of correction factor.

1101001000x0123result161281024ilog10(x)estimate

The estimate lags each power of ten until the next power of two. Shaded: wrong answers.

And a correction factor we will add; the simplest way of doing so is to just check if we are above the nearest power of ten, and if so add one. But how will we know at what magnitude above the power of ten to stop applying the correction?7 Obviously once we get close to the next power of ten, but we stop needing the correction well before then, before we even hit the midway point between two powers of ten, and applying the correction past when it is needed would become an in-correction. Furthermore, how will we even find the closest power of ten?

The solution lies in the relationship between the number we are correcting and the nearest power of ten. Namely the number we are trying to correct, plus one, is the exponent of the power of ten that is the only one xxx could be over and need correction. This is because at the exact number we would need to stop applying the correction, meaning our estimate has caught up, our estimate would increment and we would be comparing against the next power of ten which we would then be lower than.

891011121314151617181920x10100checked against 10checked against 100+1 applied[1] estimate 0 => 1[2] check jumps 10 => 100[3] 16 < 100; +1 stopsestimate: 0estimate: 1

Once the estimate catches up the check moves to the next power of ten, which x is under, and we no longer need to apply the correction.

The code for this looks like this:

#define _GNU_SOURCE #include <limits.h> #include <math.h> #include <stdint.h>

#define LONGSIZE ((int)(sizeof(int64_t) * CHAR_BIT))

uint64_t ilog10(uint64_t x) { uint64_t ilog2 = (LONGSIZE - 1) - __builtin_clzll(x); uint64_t estimate = (ilog2 * 1233) >> 12; return estimate + (x >= exp10(estimate + 1)); }

step or run

This is a big improvement and now should produce correct values, but we have introduced a new library call whose implementation we need to check, exp10:

double __exp10 (double x) { uint64_t ix = asuint64 (x); uint32_t abstop = (ix >> 52) & 0x7ff;

if (__glibc_unlikely (abstop - SmallTop >= Thresh)) { if (abstop -