4 different implementations of modulo with fully defined behavior

Question

I've implemented the 3 variants of the modulo operation described on this Wikipedia page.
The goal is to have fully defined behavior for all inputs.

Code

Implementation with truncated division (result has the same sign as the dividend)

int32_t rem(int32_t x, int32_t y) {
    if (y == 0) exit(-1);
    if (y == -1) return 0;
    return x % y;
}

This one is the simplest since it's already how the % operator is implemented in C.
We have to take care of two edge cases though:

• If y == 0 it's an invalid operation; exit(-1).
• INT32_MIN % -1 is undefined behavior, so I've hard-coded the result for y == -1.

Implementation with floored division (result has the same sign as the divisor)

int32_t mod(int32_t x, int32_t y) {
    if (y == 0) exit(-1);
    if (y == -1) return 0;
    return x - y * (x / y - (x % y && (x ^ y) < 0));
}

This one is a bit trickier.

We still have to handle the two previous edge cases. However, the modulo is computed differently: x - y * floor(x / y). We can't use floor here because we're working with integers, so the formula looks like this:

x - y * (x / y - {1 if x isn't a multiple of y and the result of x / y is negative, 0 else})
x - y * (x / y - {1 if x % y and the result of x / y is negative, 0 else})
x - y * (x / y - {1 if x % y and (x or y, but not both, is negative), 0 else})
x - y * (x / y - {1 if x % y and (x ^ y) < 0, 0 else})
x - y * (x / y - (x % y && (x ^ y) < 0))

Implementation with Euclidean division (result is always positive)

Portable implementation

int32_t euc1(int32_t x, int32_t y) {
    if (y == 0) exit(-1);
    if (y == -1) return 0;
    if (y == INT32_MIN) return x >= 0 ? x : INT32_MAX + x + 1;
    y = abs(y);
    return x - y * (x / y - (x < 0 && x % y));
}

It's basically the same implementation as the previous one, except that we use the absolute value of y. Also, (x ^ y) < 0 was simplified to x < 0 because y is always positive (we use the absolute value). However abs(INT32_MIN) is undefined behavior (because abs(INT32_MIN) = INT32_MAX + 1 which is not representable) so we have to take care of this special case:

• If x >= 0, we have x % (INT32_MAX + 1) = x because x is always inferior to INT32_MAX + 1.
• If x < 0, we have x % (INT32_MAX + 1) = INT32_MAX + 1 + x, then we reorder operations to avoid overflow on signed integers which is undefined behavior: INT32_MAX + x + 1 (remember that x is negative).

GCC implementation

int32_t euc2(int32_t x, int32_t y) {
    if (y == 0) exit(-1);
    if (y == -1) return 0;
    if (y != INT32_MIN) y = abs(y);
    int32_t tmp = x / y - (x < 0 && x % y);
    __builtin_mul_overflow(y, tmp, &tmp);
    __builtin_sub_overflow(x, tmp, &tmp);
    return tmp;
}

This implementation is the same as the portable one, except for the handling of the y = INT32_MIN edge case: there is less conditional logic, making it potentially easier to optimize.

If y == INT32_MIN we don't compute the absolute value because it's undefined behavior, then we continue the execution as nothing happened. A multiplication and a subtraction were replaced by their built-in functions counterparts to avoid undefined behavior on overflow. Because of black magic it works:

• If x >= 0 we have x / INT32_MIN - (x < 0 && x % INT32_MIN) = 0, then INT32_MIN * 0 = 0, then x - 0 = x, so the result is x as expected.
• If INT32_MIN < x < 0 we have x / INT32_MIN - (x < 0 && x % INT32_MIN) = -1, then INT32_MIN * -1 = INT32_MIN (built-in function defined behavior), then x - INT32_MIN = INT32_MAX + x + 1 (built-in function defined behavior), so the result is INT32_MAX + x + 1 as expected.
• If x == INT32_MIN we have INT32_MIN / INT32_MIN - (INT32_MIN < 0 && INT32_MIN % INT32_MIN) = 1, then INT32_MIN * 1 = INT32_MIN, then INT32_MIN - INT32_MIN = 0, so the result is 0 as expected.

Questions

Is there an error somewhere, something that i missed? Like a set of inputs that don't give the right result.

If this code really undefined behavior free?

Can it be simplified?

Answer 1

Is there an error somewhere, something that i missed? Is there an error somewhere, something that i missed?

Certainly have most of main street covered.

if (y == -1) return 0; is a nice trick.

y = abs(y); fails 32-bit math when int is 16-bit. Could use labs() or conditional code or an if() block. (Only real portability bug I see.)

Pedantic: int32_t is an optional type. It is missing and code fails to compile on select (dinosaur) machines that do not have 32-bit unpadded 2's complement.

Pedantic: On a select (unicorn) machine where int is 64-bit and int32_t exist, promotions occur. I do not see a problem though.

Pedantic: Usage of x ^ y obliges 2's complement encoding for correct functionality which is the case with intN_t here , but not all integers on all machines.

---

Adding a test harness would help.

---

Can it be simplified?

A shot at simplifying Euclidean division (result is always non-negative). Similar code may apply to the other 2 cases.

// Substitute `int` for `int32_t`, `long long`, etc.
int modulo_Euclidean(int x, int y) {
  if (y == 0) exit(-1);
  if (y == -1) return 0;
  int m = x % y;
  if (m < 0) {
    m = (y < 0) ? m - y : m + y;
  }
  return m;
}

Does not require intN_t math^*, 2's complement, no padding, range constants, nor extended math.

---

It would be interesting to see if 64-bit math was really slower or a reasonable candidate.

int32_t mod_rem_or_euc(int32_t x, int32_t y) {
  if (y == 0) exit(-1);
  int_fast64_t m = (int_fast64_t) x % y;
  ...
}

---

^* intN_t types are 2's complement, no padding.

Answer 2

The code may be free of undefined behavior, but it relies on implementation-defined behavior because of the exit(-1). Read its definition in the C standard, not in POSIX or even Linux. Of these, the C standard gives the fewest guarantees.

Your code is definitely missing a test suite that demonstrates all the edge cases. Without such a test suite, it's not apparent how much thought you invested into this code.

Answer 3

Is there an error somewhere, ... Like a set of inputs that don't give the right result.

Code lacks basic tests to answer that.

Sample test harness. (It lacks % 0 tests.)

OP's euc1() tested successfully on a 32-bit int machine.

#include <assert.h>
#include <limits.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int32_t euc1(int32_t x, int32_t y) {
  if (y == 0) exit(-1);
  if (y == -1) return 0;
  if (y == INT32_MIN) return x >= 0 ? x : INT32_MAX + x + 1;
  y = abs(y);
  return x - y * (x / y - (x < 0 && x % y));
}

int32_t euc1_ref(int32_t x, int32_t y) {
  if (y == 0) return -1;
  int64_t m = (int64_t) x % y;
  if (m < 0) {
    if (y < 0) m -= y;
    else m += y;
  }
  assert(m >= 0);
  assert(m <= INT32_MAX);
  return (int32_t) m;
}

int test_euc(int32_t x, int32_t y) {
  int32_t m1 = euc1_ref(x, y);
  int32_t m2 = euc1(x, y);
  if (m1 != m2) {
    printf("(%ld, %ld) --> %ld, %ld\n", (long) x, (long) y, (long) m1, (long) m2);
    return -1;
  }
  return 0;
}

// Could make this a little more efficient, but something to get things started.
int32_t rand_int32(void) {
  union {
    uint32_t u;
    int32_t i;
  } x;
  x.u = (unsigned) rand();
  x.u = (x.u << 15) ^ (unsigned) rand();
  x.u = (x.u << 15) ^ (unsigned) rand();
  return x.i;
}
    
int tests_euc(void) {
  int32_t t[] = { INT32_MIN, INT32_MIN + 1, -2, -1, 0, 1, 2, INT32_MAX - 1,
  INT32_MAX };
  size_t n = sizeof t / sizeof t[0];
  for (size_t yi = 0; yi < n; yi++) {
    int32_t y = t[yi];
    if (y == 0)
      continue;
    for (size_t xi = 0; xi < n; xi++) {
      int32_t x = t[xi];
      int retval = test_euc(x, y);
      if (retval) {
        return retval;
      }
    }
  }
  n = 1000u*1000*1000; // Number of random tests
  while (n) {
    int32_t y = rand_int32();
    if (y == 0) {
       continue;
    }
    int32_t x = rand_int32();
    int retval = test_euc(x, y);
    if (retval) {
      return retval;
    }
    n--;
  }
  return 0;
}

int main() {
  tests_euc();
  puts("Done");
}