## Division via Multiplication

In very special circumstances, integer division can be replaced by (integer) multiplication.

Suppose you wish to divide many B-bit numbers N by a fixed odd number D, where you know in advance that D always divides N exactly (i.e. with a zero remainder). Then you can use the following optimisation trick. Instead of dividing by D, multiply by the ‘magic’ number FD.

Example.

Consider B=64 (i.e. 64-bit numbers) with D=123=0x7B. Then it happens that FD=3449391168254631603=0x2FDE,B2FD,EB2F,DEB3.

Trying this on N=123000 (which is known to be divisible by 123, obviously):

123000 * 3449391168254631603 = 0x59D8,0000,0000,0000,03E8 → 0x3E8 = 1000.

The arrow above is loss of overflow bits

Magic!

So, how is FD ralated to D? Well, FD is just the congruent inverse of D modulo 2B. (Congruent inverses are easily calculated by one of two methods. The first is the “extended Euclidean algorithm”. The second is a reduction method.)

Why does it work? That’s because D * FD = 1 (mod M) and D and M are coprime. In our case, M=2B and D is odd, which is sufficient (but not necessary) to guarantee coprimality.

Here’s some Python to calculate congruent inverses:

```### extended Euclidean algorithm
### return (h, s, t) s.t. h = hcf (a, b), h = s * a + t * b
def congruence_solve (a, b):
if b == 0:
### return (a, 1, 0)
return (abs(a), sgn(a), 0)
else:
q = a // b
r = a % b
(h, u, v) = congruence_solve (b, r)
s = v
t = u - v * q
return (h, s, t)```
```### multiplicative inverse in modular arithmetic
### return y s.t.
###   x * y ~ 1 (mod m)
###   0 < y < m  (or  m < y < 0)
### if such a y exists; 0 otherwise
def congruent_inverse_1 (x, m):
(unit, y, t) = congruence_solve (x, m)
if unit == 1:
y = y % m
if (m > 0 and y < 0) or (m < 0 and y > 0):
y += m
return y
else:
return 0```