2
\$\begingroup\$

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
\$\endgroup\$
1
\$\begingroup\$

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");
}
\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$
3
  • \$\begingroup\$ 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. \$\endgroup\$ – chux - Reinstate Monica Jan 2 at 3:24
  • \$\begingroup\$ As for portability, uintmax_t may be more than 64-bit, so code should not assume 64. \$\endgroup\$ – chux - Reinstate Monica Jan 2 at 3:26
  • \$\begingroup\$ Of course; 64 is short-hand for "whatever your working integer size may be", and yes, it should not be hard-coded. \$\endgroup\$ – Reinderien Jan 2 at 5:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.