Portable integer log2 using floating point primitives

[ewx]

One of the things I needed for the user management window below was to convert from single bit in a bitmap to its bit number. (Rights are represented internally as a bitmap, but you want the bit number to index the array of GtkCheckButton widgets that represented them.)

Henry Warren suggests (among a number of other things) using the system's ability to perform this operation as part of converting integers to normalized floating point values, and inferring the answer from the exponent it chooses. Unfortunately the code he presents is very unportable, as it explicitly inspects the bytes representing the floating point value. His code is:

int nlz6(unsigned k) {
   union {
      unsigned asInt[2];
      double asDouble;
   };
   int n;

   asDouble = (double)k + 0.5;
   n = 1054 - (asInt[LE] >> 20);
   return n;
}

(nlz stands for “number of leading zeroes”. In a fixed-size word this is essentially the same question as leftmost bit or log₂.)

As it happens C has long had a function called frexp which will portably extract the exponent, eliminating the need to make any assumptions about the representation of double. The code is very simple:

/** @brief Compute index of leftmost 1 bit
 * @param n Integer
 * @return Index of leftmost 1 bit or -1
 *
 * For positive @p n we return the index of the leftmost bit of @p n.  For
 * instance @c leftmost_bit(1) returns 0, @c leftmost_bit(15) returns 3, etc.
 *
 * If @p n is zero then -1 is returned.
 */
int leftmost_bit(uint32_t n) {
  /* See e.g. Hacker's Delight s5-3 (p81) for where the idea comes from.
   * Warren is computing the number of leading zeroes, but that's not quite
   * what I wanted.  Also this version should be more portable than his, which
   * inspects the bytes of the floating point number directly.
   */
  int x;
  frexp((double)n, &x);
  /* This gives: n = m * 2^x, where 0.5 <= m < 1 and x is an integer.
   *
   * If we take log2 of either side then we have:
   *    log2(n) = x + log2 m
   *
   * We know that 0.5 <= m < 1 => -1 <= log2 m < 0.  So we floor either side:
   *
   *    floor(log2(n)) = x - 1
   *
   * What is floor(log2(n))?  Well, consider that:
   *
   *    2^k <= z < 2^(k+1)  =>  floor(log2(z)) = k.
   *
   * But 2^k <= z < 2^(k+1) is the same as saying that the leftmost bit of z is
   * bit k.
   *
   *
   * Warren adds 0.5 first, to deal with the case when n=0.  However frexp()
   * guarantees to return x=0 when n=0, so we get the right answer without that
   * step.
   */
  return x - 1;
}

The code should certainly work up to the number of mantissa bits in a double, but I would worry about rounding causing 2ⁿ-1 to overflow to 2ⁿ for values larger than around 2⁵³. I suppose this could be checked manually, avoiding including knowledge of how the rounding is actually done, but I don't need anything bigger than 32 bits right now anyway.

(comments locked due to excess spam which for some reason targets this article in particular.)

Comments

[gjm11.livejournal.com]

As mentioned above, getting your hands on that particular constant isn't so very bad on an ARM (more precisely, you don't need to get the constant as such, you just need to multiply by it).

If the speed of this thing actually matters, it's probably being done a lot, so the relevant bit of your literal pool is likely to be in the dcache.

(For the avoidance of doubt, I'm not claiming that Seal's approach is best on all architectures. I just don't think it's so narrowly targetted as to justify calling it "ARM-centric".)