# 4 different implementations of modulo with fully defined behavior

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?

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.

• Thanks for this review. I like your Euclidean modulo simplification because it gets rid of the abs function that can cause problem when fed with int32_t. Dec 23, 2020 at 12:31

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.

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");
}