Efficient long division via Montgomery multiply

TL;DR

This paper introduces a new efficient long division algorithm using Montgomery multiplication, optimized for multi-word dividends and large moduli, with significant performance improvements on modern hardware.

Contribution

The paper presents a novel right-to-left long division algorithm based on Montgomery multiply, optimized for large multi-word dividends and efficient hardware implementation.

Findings

Achieved 6 cycles per dividend word for remainder calculation on x86_64.

Achieved 12.5 cycles per dividend word for full division (quotient and remainder).

Introduced a bit-doubling modular inversion scheme for large precision inverse computation.

Abstract

We present a novel right-to-left long division algorithm based on the Montgomery modular multiply, consisting of separate highly efficient loops with simply carry structure for computing first the remainder (x mod q) and then the quotient floor(x/q). These loops are ideally suited for the case where x occupies many more machine words than the divide modulus q, and are strictly linear time in the "bitsize ratio" lg(x)/lg(q). For the paradigmatic performance test of multiword dividend and single 64-bit-word divisor, exploitation of the inherent data-parallelism of the algorithm effectively mitigates the long latency of hardware integer MUL operations, as a result of which we are able to achieve respective costs for remainder-only and full-DIV (remainder and quotient) of 6 and 12.5 cycles per dividend word on the Intel Core 2 implementation of the x86_64 architecture, in single-threaded execution mode. We further describe a simple "bit-doubling modular inversion" scheme, which allows the entire iterative computation of the mod-inverse required by the Montgomery multiply at arbitrarily large precision to be performed with cost less than that of a single Newtonian iteration performed at the full precision of the final result. We also show how the Montgomery-multiply-based powering can be efficiently used in Mersenne and Fermat-number trial factorization via direct computation of a modular inverse power of 2, without any need for explicit radix-mod scalings.