Fast division, modulus and divisibility test with constant divider

div_const.py

#!/usr/bin/env python

# The algorithms implemented in this module are based on the following article:
# D. Lemire, O. Kaser, and N. Kurz, Faster Remainder by Direct Computation, 2018.

import math


def div(value, reciprocal, fractional_bits):
    return (value * reciprocal) >> fractional_bits


def mod(value, reciprocal, fractional_bits):
    fractional_mask = (1 << fractional_bits) - 1
    return (((value * reciprocal) & fractional_mask) * divisor) >> fractional_bits


def is_divisible(value, reciprocal, fractional_bits):
    fractional_mask = (1 << fractional_bits) - 1
    return (value * reciprocal) & fractional_mask < reciprocal


def generate_div_params(divisor, bits):
    """
    Calculate parameters for fast fixed point division by constant value.

    Returned value is a tuple (reciprocal, fractional_bits).

    The following formula is used to calculate division:
    (x / divisor) mod (1 << bits) => (x * reciprocal) >> fractional_bits

    The following formula is used to calculate remainder:
    fractional_mask = (1 << fractional_bits) - 1
    x mod divisor => (((x * reciprocal) & fractional_mask) * divisor) >> fractional_bits

    The following formula is used to check for divisibility:
    fractional_mask = (1 << fractional_bits) - 1
    x mod divisor == 0 => (x * reciprocal) & fractional_mask < reciprocal
    """

    assert(0 < divisor < (1 << bits))

    if is_power_of_2(divisor):
        reciprocal = 1
        fractional_bits = int(math.log2(divisor))
    else:
        for extra_bits in range(0, bits + 1):
            fractional_bits = bits + extra_bits
            val_range = 1 << fractional_bits
            reciprocal = ((val_range - 1) // divisor) + 1
            if divisor <= val_range % divisor + (1 << extra_bits):
                break
        else:
            raise AssertionError("Failed to calculate number of fractional bits")

    return (reciprocal, fractional_bits)


def is_power_of_2(x):
    return x != 0 and x & (x - 1) == 0


if __name__ == '__main__':
    for bits in range(1, 16 + 1):
        print("{} bits".format(bits))
        value_range = 1 << bits

        for divisor in range(1, value_range):
            reciprocal, fractional_bits = generate_div_params(divisor, bits)

            print("x / {} mod {} => x * {} >> {}".format(
                divisor, value_range, reciprocal, fractional_bits))

            for value in range(value_range):
                assert(div(value, reciprocal, fractional_bits) == value // divisor)
                assert(mod(value, reciprocal, fractional_bits) == value % divisor)
                assert(is_divisible(value, reciprocal, fractional_bits) == (value % divisor == 0))