Unsigned integer power with overflow detection

Question

An unsigned integer version of pow() is fairly easy to code, yet I also wanted:

1. overflow detection
2. portability
3. only modest performance impact due to overflow detection.

The below codes uses a simple >= test in the loop to see if the calculation is nearing overflow. When needed, in a 2nd loop, more expensive / tests are used.

[ Code also uses a nifty IMAX_BITS macro to find the number of value bits in a ..._MAX constant as CHAR_BIT * sizeof(type) fails when rare padding exist or the MAX value is not the MAX of the type, like RAND_MAX. ]

Review goals:

• Performance improvement ideas for uintmax_t pow_umax(uintmax_t base, unsigned expo), especially concerning overflow detection.

• General review about pow.h and pow.c.

• General review about test code functionally, not code.

---

pow.h header

#ifndef POW_H
#define POW_H 1
#include <inttypes.h>

// Return base raised to the expo power.
// On overflow, set errno = ERANGE and return UINTMAX_MAX.
uintmax_t pow_umax(uintmax_t base, unsigned expo);

#endif

pow.c code

#include "pow.h"
#include <errno.h>
#include <inttypes.h>
#include <limits.h>

// https://stackoverflow.com/a/4589384/2410359
/* Number of value bits in inttype_MAX, or in any (1<<k)-1 where 0 <= k < 2040 */
#define IMAX_BITS(m) ((m)/((m)%255+1) / 255%255*8 + 7-86/((m)%255+12))

// Value bit width of UINTMAX_MAX (handles padded types).
#define UINTMAX_BITS IMAX_BITS(UINTMAX_MAX)

// When `base >= cubert(UINTMAX_MAX + 1)` potentially the next
// `base *= base` and `pow *= base` may overflow.
// Course lower bound cube root of (UINTMAX_MAX + 1)
#define UINTMAX_MAXP1_CBRT (((uintmax_t) 1) << (UINTMAX_BITS/3))

// Detect overflow in an integer pow() calculation with minimal overhead.
uintmax_t pow_umax(uintmax_t base, unsigned expo) {
  uintmax_t pow = 1;

  for (;;) {
    if (expo % 2) {
      pow *= base;
    }
    expo /= 2;
    if (expo == 0) {
      return pow;
    }
    // If base is so big, it might overflow `base *= base` or `pow *= base`.
    if (base >= UINTMAX_MAXP1_CBRT) {
      break;
    }
    base *= base;
  }

  // `base` is large here, add precise tests for `base *= base` and `pow *= base`.
  for (;;) {
    if (base > UINTMAX_MAX / base) {
      break;
    }
    base *= base;
    if (expo % 2) {
      if (pow > UINTMAX_MAX / base) {
        break;
      }
      pow *= base;
    }
    expo /= 2;
    if (expo == 0) {
      return pow;
    }
  }

  errno = ERANGE;
  return UINTMAX_MAX;
}

Test code

#include <assert.h>
#include <errno.h>
#include <inttypes.h>
#include <limits.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include "pow.h"

// Bit width of various ..._MAX (handles padded types).
#define IMAX_BITS(m) ((m)/((m)%255+1) / 255%255*8 + 7-86/((m)%255+12))
#define UINTMAX_BITS IMAX_BITS(UINTMAX_MAX)
#define RAND_MAX_BITS IMAX_BITS(RAND_MAX)
#define UNSIGNED_BITS IMAX_BITS(UINT_MAX)

// Random number usage on this.
_Static_assert(((RAND_MAX + 1u) & RAND_MAX) == 0, "RAND_MAX is not a Mersenne number");
// Support for failing implementations is not a review goal.

// Return random number [0...UINTMAX_MAX]
uintmax_t rand_umax(void) {
  uintmax_t u = 0;
  for (unsigned bits = 0; bits < UINTMAX_BITS; bits += RAND_MAX_BITS) {
    u = (u << RAND_MAX_BITS) ^ (rand() & RAND_MAX);
  }
  return u;
}

// Return random number [0...UINT_MAX]
unsigned rand_u(void) {
  unsigned u = 0;
  for (unsigned bits = 0; bits < UNSIGNED_BITS; bits += RAND_MAX_BITS) {
    u = (u << RAND_MAX_BITS) ^ (rand() & RAND_MAX);
  }
  return u;
}

_Static_assert(sizeof(long double)*CHAR_BIT >= 80, "Need wider long double for verification");

// Test pow_umax() for errno and return value correctness.
int pow_umax_test(uintmax_t base, unsigned expo) {
  errno = 0;
  uintmax_t p = pow_umax(base, expo);
  if (errno && p != UINTMAX_MAX) {
    printf("Error: %ju, %u", base, expo);
    exit(EXIT_FAILURE);
  }
  // Use powl() as a reference function.
  long double y = powl(base, expo);
  if ((errno == 0 && y != p) || (errno && y <= UINTMAX_MAX)) {
    printf("Error: %ju, %u --> p=%ju, y=%Lf", base, expo, p, y);
    exit(EXIT_FAILURE);
  }
  return 0;
}

// Perform multiple pow_umax() tests.
int pow_umax_tests(void) {
  printf("uintmax_t bit width: %ju\n", UINTMAX_BITS);
  printf("unsigned bit width : %u\n", UNSIGNED_BITS);
  printf("RAND_MAX bit width : %u\n", RAND_MAX_BITS);
  printf("uintmaxp1_cbrt     : 0x%jX\n", UINTMAX_MAXP1_CBRT);
  int err = 0;

  // The usual suspects.
  uintmax_t b[] = {0, 1, 2, UINTMAX_MAXP1_CBRT - 1, UINTMAX_MAXP1_CBRT,
  UINTMAX_MAXP1_CBRT + 1, UINTMAX_MAX - 1, UINTMAX_MAX};
  unsigned e[] = {0, 1, 2, 3, 4, 5, 6, 7, 8, UINT_MAX - 1, UINT_MAX};
  size_t bn = sizeof b / sizeof b[0];
  size_t en = sizeof e / sizeof e[0];
  for (unsigned bi = 0; bi < bn; bi++) {
    for (unsigned ei = 0; ei < en; ei++) {
      err += pow_umax_test(b[bi], e[ei]);
    }
  }

  // Random tests
  for (unsigned long i = 0; i < 10000000u; i++) {
    err += pow_umax_test(rand_umax(), rand_u());
  }
  return err;
}

int main(void) {
  pow_umax_tests();
  puts("Done");
}

Sample result (no errors)

uintmax_t bit width: 64
unsigned bit width : 32
RAND_MAX bit width : 31
uintmaxp1_cbrt     : 0x200000
Done

Answer 1

General review about pow.h and pow.c.

When overflow does not occur, pow_umax() can only return UINTMAX_MAX with pow_umax(UINTMAX_MAX, 1).

There is no other base, expo such that pow_umax(base, expo) returns UINTMAX_MAX, regardless of the bit width of uintmax_t due to Mihailescu's Theorem. [Just found that out today.]

So setting errno = ERANGE is only distinctive in that one case. Otherwise code could ignore errno and use:

// Sufficient test except in one case.
if (pow_umax(base, expo) == UINTMAX_MAX) {
  puts("Overflow");
}

Answer 2

As a discussion point (and not strictly a recommendation), I wonder how your overflow detection algorithm would compare with the following:

Assuming a 64-bit integer (substituting as necessary): We know that if \$y \ge 2^{64} \$, \$y\$ has overflowed. If \$ y = x^a \$,

$$ x^a \lt 2^{64} $$ $$ a \log_2 x < 64 $$

\$ \log_2 \$ is cheap and can be done with bit-twiddling, i.e. without using floating-point math; so this should be quick enough. For details on how to do this see https://stackoverflow.com/questions/11376288/fast-computing-of-log2-for-64-bit-integers .

This can be a quick first check before performing the full check, but cannot be a comprehensive check on its own. Take for instance

$$ 3^{63} >> 2^{64} $$ $$ \left\lfloor \log_2 3 \right\rfloor = 1 $$ $$ 63 \times 1 < 64$$

The rounding error from calculating an integer logarithm can yield false negatives.