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.

• chux: When a*log2(x) is close to 64, the calculation of log2(x) and multiplication by a needs to be done very carefully and prone to slight errors. When that product is 64 + margin there's overflow & 64 - margin there's no overflow. Yet when ~64 the integer calculation remains necessary (with overflow precautions) for the remains corner cases. As for "log2 is cheap and can be done with bit-twiddling", (integer part certainly easy) I'd be interest in seeing that. It may need the answer calculate as not +/- 0.5 ULP, but +0/-1.0 ULP (one sided error). Still an insightful idea.
• chux: As for portability, uintmax_t may be more than 64-bit, so code should not assume 64.
• Reinderien: Of course; 64 is short-hand for "whatever your working integer size may be", and yes, it should not be hard-coded.