div

div.c

#include <fcntl.h>
#include <inttypes.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>

// The highest bit in the uintmax_t type to avoid shift overflow.
#define UINTMAX_HIGH (UINTMAX_C(1) << (sizeof(uintmax_t) * 8 - 1))

// Attempt 1:
// Naively divide by counting how many times b can be subtracted from a. This
// algorithm runs in O(a/b) time.
uintmax_t divide1(uintmax_t a, uintmax_t b)
{
	uintmax_t result = 0;
	while ( b <= a )
	{
		a -= b;
		result++;
	}
	return result;
}

// Attempt 2:
// a/b = 2 * (a / 2b) + (a % 2b) / b.
// Recursively reduce the problem to a division with double b and correct the
// result. The remainder of a is less than 2b, and the remaining division can
// be done with a single conditional subtraction. This algorihm runs in time
// O(log2(a/b)) which is O(32) or O(64) on systems with common fixed size
// integers.
uintmax_t divrem2(uintmax_t a, uintmax_t b, uintmax_t* rem)
{
	uintmax_t result = 0;
	if ( !(b & UINTMAX_HIGH) )
	{
		uintmax_t bb = b << 1;
		if ( bb <= a )
			result += 2 * divrem2(a, bb, &a);
	}
	if ( b <= a )
	{
		a -= b;
		result++;
	}
	*rem = a;
	return result;
}

uintmax_t divide2(uintmax_t a, uintmax_t b)
{
	return divrem2(a, b, &a);
}

// Attempt 3:
// Unroll the recursion of attempt 2.
uintmax_t divide3(uintmax_t a, uintmax_t b)
{
	uintmax_t denom = b;
	uintmax_t magnitude = 1;
	while ( !(denom & UINTMAX_HIGH) && denom < a )
	{
		denom = denom << 1;
		magnitude = magnitude << 1;
	}
	uintmax_t result = 0;
	while ( magnitude )
	{
		if ( denom <= a )
		{
			a -= denom;
			result += magnitude;
		}
		magnitude = magnitude >> 1;
		denom = denom >> 1;
	}
	return result;
}

// Attempt 4:
// 2a / 2b = a / b.
// Repeatedly remove common factors in the numerator and denominator, then
// invoke attempt 3.
uintmax_t divide4(uintmax_t a, uintmax_t b)
{
	while ( a && b && !(a & 1) && !(b & 1) )
	{
		a = a >> 1;
		b = b >> 1;
	}
	return divide3(a, b);
}

// Attempt 5:
// When the numerator gets small enough in attempt 3, use 32-bit divison.
uintmax_t divide5(uintmax_t a, uintmax_t b)
{
	if ( a < b )
		return 0;
	if ( a <= UINT32_MAX )
		return (uint32_t) a / (uint32_t) b;
	uintmax_t denom = b;
	uintmax_t magnitude = 1;
	while ( !(denom & UINTMAX_HIGH) && denom < a )
	{
		denom = denom << 1;
		magnitude = magnitude << 1;
	}
	uintmax_t result = 0;
	while ( magnitude )
	{
		if ( denom <= a )
		{
			a -= denom;
			result += magnitude;
			if ( a <= UINT32_MAX )
			{
				if ( a < b )
					return result;
				return result + (uint32_t) a / (uint32_t) b;
			}
		}
		magnitude = magnitude >> 1;
		denom = denom >> 1;
	}
	return result;
}

// Test facility for the chosen attempt.
#define DIVIDE divide5
void test(uintmax_t a, uintmax_t b)
{
	uintmax_t c = DIVIDE(a, b);
	uintmax_t d = a / b;
	if ( c != d )
	{
		printf("WRONG! Got %ju / %ju = %ju but correct is %ju\n", a, b, c, d);
		exit(1);
	}
	printf("%ju / %ju = %ju\n", a, b, c);
}

// Command line program that divides two unsigned integers and prints the
// result, or if no parameters are given, starts doing random divisions.
int main(int argc, char* argv[])
{
	if ( argc <= 1 )
	{
		int fd = open("/dev/urandom", O_RDONLY);
		while ( 1 )
		{
			uintmax_t a = 0;
			uintmax_t b = 0;
			read(fd, &a, sizeof(a));
			while ( b == 0 )
			{
				read(fd, &b, sizeof(b));
				// Uncomment and modify to skew the distribution:
				//b &= 0xFF;
			}
			test(a, b);
		}
		close(fd);
		return 0;
	}
	if ( argc != 3 )
	{
		fprintf(stderr, "Usage: %s numerator denominator\n", argv[0]);
		return 1;
	}
	uintmax_t a = strtoumax(argv[1], NULL, 0);
	uintmax_t b = strtoumax(argv[2], NULL, 0);
	if ( b == 0 )
	{
		printf("Undefined behavior dividing by zero\n");
		return 1;
	}
	test(a, b);
	return 0;
}