GCC Eliminates Unnecessary Integer Division

Here are free GPU credits. First come first serve :)

A newsletter about applied math and theoretical computer science. Subscribe!

We code the paper on Fast Unsigned Division by Constants (Ammon, 2011) for our number field sieve.

1.0 Introduction

The number field sieve demands us perform repetitive divisibility tests. Observe that GCC runs faster when provided a constant divisor, not a variable, as shown below.

ConstantMod: 0.454 s
VariableMod: 0.562 s

Turns out, GCC finds two magic numbers and replaces division with one multiplication and a comparison. This optimization is Part 8 in our Practical Number Field Sieve for Programmers series:

Part 1: Discrete Logarithms and the Index Calculus Solution.

Part 2: Solving Index Calculus Equations over Integers and Finite Fields.

Part 3: Computation of Discrete Logarithms in Algebraic Number Fields.

Part 4: 2 Dimensional Lattice Basis Reduction.

Part 5: Individual Reduction Phase and Logarithm Collection.

Part 6: Lattice Sieving and Special Q Descent.

Part 7: Continued Fractions and Fast Lattice Sieving.

Part 8: Optimizing Prime Division using Magic Numbers.

1.1 Comparison Code To Observe Assembly

First we write this code for comparison:

#include 
#include 
#include 
#include 
#include 
//Generate Asm: clear && gcc -O2 -S -masm=intel ZeroModPrime.c -o ZeroModPrime.s
//Run code: clear && gcc ZeroModPrime.c -lm -o m.o && ./m.o

#define Constant 17
int ConstantMod(unsigned x)
{
    return x % Constant == 0;
}

int VariableMod(unsigned x, unsigned var)
{
    return x % var == 0;
}

void TimeModulus(void)
{
	const int N = 100000000;
	clock_t start;
	start = clock();
	for(int i = 0; i < N; i++)ConstantMod(i);
	printf("ConstantMod: %.3f s\n",(double)(clock() - start) / CLOCKS_PER_SEC);

	start = clock();
	for(int i = 0; i < N; i++)VariableMod(i, Constant);
	printf("VariableMod: %.3f s\n",(double)(clock() - start) / CLOCKS_PER_SEC);
}


int main()
{
	TimeModulus();
        return(0);
}

Next, we observe the corresponding Assembly:

#Observe the comparison
ConstantMod:
.LFB50:
	.cfi_startproc
	endbr64
	imul	edi, edi, -252645135
	xor	eax, eax
	cmp	edi, 252645135
	setbe	al
	ret
	.cfi_endproc
#Observe the division
VariableMod:
.LFB51:
	.cfi_startproc
	endbr64
	mov	eax, edi
	xor	edx, edx
	div	esi
	xor	eax, eax
	test	edx, edx
	sete	al
	ret
	.cfi_endproc

The variable modulo function explicitly calls div while constant modulo for 17 does x%17==0; // x*0xF0F0F0F1 <= 0x0F0F0F0F.

Our techbro ascendants documented this behavior on StackOverflow (Technosaurus, 2018):

// Source - https://stackoverflow.com/q/53414711
// Posted by technosaurus
// Retrieved 2026-09-14, License - CC BY-SA 4.0

//32bit examples for _Bool mod_n(unsigned x){return x%n==0;};
//note: parameter is unsigned but it becomes a signed multiply
x%3==0;  // x*0xAAAAAAAB <= 0x55555555
x%5==0;  // x*0xCCCCCCCD <= 0x33333333
x%7==0;  // x*0xB6DB6DB7 <= 0x24924924
x%11==0; // x*0xBA2E8BA3 <= 0x1745D174
x%13==0; // x*0xC4EC4EC5 <= 0x13B13B13
x%17==0; // x*0xF0F0F0F1 <= 0x0F0F0F0F
x%19==0; // x*0x286BCA1B <= 0x0D79435E
x%23==0; // x*0xE9BD37A7 <= 0x0B21642C
x%29==0; // x*0x4F72C235 <= 0x08D3DCB0
x%31==0; // x*0xBDEF7BDF <= 0x08421084
x%37==0; // x*0x914C1BAD <= 0x06EB3E45
x%41==0; // x*0xC18F9C19 <= 0x063E7063
x%43==0; // x*0x2FA0BE83 <= 0x05F417D0
x%47==0; // x*0x677D46CF <= 0x0572620A
x%53==0; // x*0x8C13521D <= 0x04D4873E
x%59==0; // x*0xA08AD8F3 <= 0x0456C797
x%61==0; // x*0xC10C9715 <= 0x04325C53
x%67==0; // x*0x07A44C6B <= 0x03D22635
x%71==0; // x*0xE327A977 <= 0x039B0AD1
x%73==0; // x*0xC7E3F1F9 <= 0x0381C0E0
x%79==0; // x*0x613716AF <= 0x033D91D2
x%83==0; // x*0x2B2E43DB <= 0x03159721
x%89==0; // x*0xFA3F47E9 <= 0x02E05C0B
x%97==0; // x*0x5F02A3A1 <= 0x02A3A0FD
///...and even up to 64bit
x%4294967291==0; //x*0x70A3D70A33333333 <= 0x100000005

Even more fascinating, GCC performs the same optimization for regular integer division, as demonstrated in (Qiubit, 2016). An interesting question arises: how do we find these magic numbers?

1.1 Naively Finding Magic Numbers

Reciprocal multiplication is the succint term for ‘replacing division with magic number multiplications and comparisons’. Our secondary sources for the remainder of this section are (Jones, 1999) and (Ammon, 2010).

These magic numbers are precomputed quotients of a rounded up power of 2. At runtime, one multiplies by the magic number, and then divides by the power of 2 to round down. Mathematically, it resembles:

We follow (Ammon, 2010) and call our magic number, mexact and define it below:

mexact is always a fraction. So we must round it up to an integer, and this rounded up integer becomes our magic number. This is called the round-up algorithm.

(Ammon, 2010) proves why magic numbers work and provides the round-up algorithm alongside division below:

1.2 Finding Magic Numbers in Prod

(Ammon, 2011) introduces the round-down algorithm to address the primary shortcoming of the round-up algorithm from Section 1.1: some magic numbers are 33 bits so don’t fit in typical 32-bit hardware register.

The round-down algorithm is guaranteed to find a 32-bit magic number when the round-up algorithm fails. It is provided below:

We follow the LibDivide reference and write this code to find magic numbers and further test if an integer n modulo a divisor d is 0:

#include 
#include 
#include 
#include 
#include 
#include 

//clear && gcc -O2 -S -masm=intel ZeroModPrime.c -o ZeroModPrime.s
//clear && gcc ZeroModPrime.c -lm -o m.o && ./m.o
struct magic_number_struct
{
	uint32_t magicNumber;
	uint32_t dividendShiftPreMultiply;
	uint32_t dividendShiftPostMultiply;
	uint32_t increment; // 0 or 1; if set then increment the numerator, using one of the two strategies
};

//Function definitions
struct magic_number_struct UnsignedMagicInfo(uint32_t divisor, uint32_t numberOfBits);
struct magic_number_struct UnsignedMagicInfo(uint32_t divisor, uint32_t numberOfBits)
{
	//The numerator must fit in a uint32_t
	assert(numberOfBits > 0 && numberOfBits <= sizeof(uint32_t) * CHAR_BIT);
	//divisor must be larger than zero and not a power of 2
	assert(divisor & (divisor-1));
	struct magic_number_struct result;
	const unsigned UINT_BITS = sizeof(uint32_t) * CHAR_BIT;
	// The extra shift implicit in the difference between UINT_BITS and numberOfBits
	const unsigned extraShift = UINT_BITS - numberOfBits;

	// The initial power of 2 is one less than the first one that can possibly work
	const uint32_t initialPowerOf2 = (uint32_t)1 << (UINT_BITS-1);

	// The remainder and quotient of our power of 2 divided by d
	uint32_t quotient = initialPowerOf2 / divisor, remainder = initialPowerOf2 % divisor;

	// ceil(log_2 divisor)
	unsigned ceilLog2divisor;

	// The magic info for the variant "round down" algorithm
	uint32_t downMultiplier = 0;
	unsigned downExponent = 0;
	int hasMagicDown = 0;

	//Compute ceil(log_2 divisor)
	ceilLog2divisor = 0;uint32_t tmp;for(tmp = divisor; tmp > 0; tmp >>= 1)ceilLog2divisor += 1;
	
	// Begin a loop that increments the exponent, until we find a power of 2 that works.
	unsigned exponent;
	for (exponent = 0; ; exponent++)
	{
		// Quotient and remainder is from previous exponent; compute it for this exponent.
		if(remainder >= divisor - remainder)
		{
			// Doubling remainder will wrap around divisor
			quotient = quotient * 2 + 1;
			remainder = remainder * 2 - divisor;
		}
		else
		{
			// Remainder will not wrap
			quotient = quotient * 2;
			remainder = remainder * 2;
		}

		// We're done if this exponent works for the round_up algorithm.
		// Note that exponent may be larger than the maximum shift supported,
		// so the check for >= ceilLog2divisor is critical.
		if ((exponent + extraShift >= ceilLog2divisor) || (divisor - remainder) <= ((uint32_t)1 << (exponent + extraShift)))
		break;

		// Set magic_down if we have not set it yet and this exponent works for the round_down algorithm
		if (! hasMagicDown && remainder <= ((uint32_t)1 << (exponent + extraShift)))
		{
			hasMagicDown = 1;
			downMultiplier = quotient;
			downExponent = exponent;
		}
	}

	if(exponent < ceilLog2divisor) 
	{
		// magic_up is efficient
		result.magicNumber = quotient + 1;
		result.dividendShiftPreMultiply = 0;
		result.dividendShiftPostMultiply = exponent;
		result.increment = 0;
	}
	else if(divisor & 1)
	{
		// Odd divisor, so use magic_down, which must have been set
		assert(hasMagicDown);
		result.magicNumber = downMultiplier;
		result.dividendShiftPreMultiply = 0;
		result.dividendShiftPostMultiply = downExponent;
		result.increment = 1;
	}
	else
	{
		// Even divisor, so use a prefix-shifted dividend
		unsigned preShift = 0;
		uint32_t shiftedDivisor = divisor;
		while ((shiftedDivisor & 1) == 0)
		{
			shiftedDivisor >>= 1;
			preShift += 1;
		}
		result = UnsignedMagicInfo(shiftedDivisor, numberOfBits - preShift);
		assert(result.increment == 0 && result.dividendShiftPreMultiply == 0); //expect no increment or preShift in this path
		result.dividendShiftPreMultiply = preShift;
	}
	return result;
}
    
static int MagicNumberDivisibilityTest(uint32_t number, uint32_t divisor, uint32_t magic, uint32_t preShift, uint32_t postShift, int increment)
{
	uint32_t original = number;
	if(preShift) number >>= preShift;
	if(increment && number != UINT32_MAX) number++;
	uint32_t q = (uint32_t)(((uint64_t)number * magic) >> 32);
	q >>= postShift;
	return original == q * divisor;
}


void TestMagicNumber()
{
	uint32_t number  = 99;
	uint32_t divisor = 7;
	uint32_t numberOfBits = sizeof(uint32_t) * CHAR_BIT;
	struct magic_number_struct magic = UnsignedMagicInfo(divisor, numberOfBits);
	printf("%u\n",numberOfBits);
	printf("Unsigned division by %u\n", divisor);
	printf("Multiplier : %u\n", magic.magicNumber);
	printf("Pre-shift  : %u\n", magic.dividendShiftPreMultiply);
	printf("Post-shift : %u\n", magic.dividendShiftPostMultiply);
	printf("Increment  : %d\n", magic.increment);
	int result = MagicNumberDivisibilityTest(number, divisor, magic.magicNumber, magic.dividendShiftPreMultiply, magic.dividendShiftPostMultiply, magic.increment);
	printf("result: %s\n", result ? "divisible" : "not divisible");
}
int main()
{
	//TimeModulus();
	TestMagicNumber();
        return(0);
}

An example with 7 is:

Unsigned division by 7
Multiplier : 1227133513
Pre-shift  : 0
Post-shift : 1
Increment  : 1

Try it yourself!

References

Ammon, P. (2011). Labor of Division (Episode 3): Fast Unsigned Division by Constants. PDF.

Technosaurus. (2018). Math Behind gcc9+ Modulus Optimizations. StackOverflow. Link.

Qiubit. (2016). Why Does GCC Use Multiplication by a Strange Number in Implementing Division. StackOverflow. Link.

Jones, D. (1999). Reciprocal Multiplication. University of Iowa Department of Computer Science. Arithmetic Tutorial Collection. Link.

Ammon, P. (2010). Labor of Division (Episode 1). Ridiculous Fish. Link.

Libdivide Authors. (2011). divide_by_constants_codegen_reference. GitHub. Link.