Three Binary Algorithms
January 15, 2010
Today’s exercise examines three binary algorithms, for basic arithmetic: multiplication, division, and greatest common divisor.
Four millenia ago, when the ancient Egyptians were building the pyramids, they invented a method for multiplying positive integers that works like this: They start by writing the two numbers to be multiplied at the head of two columns. Then they repeatedly divide the number in the left column by two (right-shift), dropping any remainder, and double the number in the right column (left-shift), writing the two new numbers immediately below their predecessors, until the number in the left column is one. Then they cross out all rows where the number in the left column is even, and add the remaining numbers in the right column, which is the desired product. For instance, the product fourteen times twelve is found like this:
14 12
7 24
3 48
1 96
168
It is easy to see why this method works if you use the grade-school method of multiplication, but with binary numbers instead of decimal numbers:
1 1 0 0 1 2
1 1 1 0 1 4
- - - - - -
0 0 0 0 0 x 1 2 0
1 1 0 0 2 x 1 2 2 4
1 1 0 0 4 x 1 2 4 8
1 1 0 0 8 x 1 2 9 6
- - - - - - - - - - - - - -
1 0 1 0 1 0 0 0 1 4 x 1 2 1 6 8
In binary, 14 is 1 · 23 + 1 · 22 + 1 · 21 + 0 · 20. The odd numbers in the left column correspond to the 1 bits in the binary representation of the multiplicand.
Binary division is the inverse operation to binary multiplication, and is also done by a series of doubling and halving operations. We’ll begin with the grade-school algorithm, in binary:
1 0 0 1 1
- - - - - - - - - - -
1 0 1 0 1 1 ) 1 1 0 1 0 0 0 1 0 1 8 3 7
1 0 1 0 1 1 6 8 8 / 4 3 = 1 6
- - - - - - - - -
1 0 0 1 0 1 0 1 4 9
1 0 1 0 1 1 8 6 / 4 3 = 2
- - - - - - - - - -
1 1 1 1 1 1 6 3
1 0 1 0 1 1 4 3 / 4 3 = 1
- - - - - - - -
1 0 1 0 0 2 0
The quotient has a 1 wherever the divisor is less than the current portion of the dividend, and a 0 where the next digit is “brought down” from the dividend. The algorithmic version of n ÷ d is based on a quantity d · 2k, which we call t, which is initialized with the largest k that leaves t less than n; this can be done by repeated left shifts. Then the quotient q is initialized to zero, the remainder r is initialized to n, and an iteration continues halving t until t is less than the original d; at each stage of the iteration, the quotient is doubled (left-shift) if the current remainder is less than t, otherwise the quotient is doubled (left-shift), then one is added, and the current remainder is reduced by the current value of t. An example is shown below:
t q r
43
86
172
344
688
1376
688 0 837
344 1 149
172 2 149
86 4 149
43 9 63
43/2 19 20
Josef Stein described a binary algorithm for greatest common divisor in 1967 that was devised by Chinese mathematicians in the first century. The greatest common divisor of two numbers is calculated by a recursive algorithm with four rules:
- If either number is zero, their greatest common divisor is the other number.
- If both numbers are even, the two numbers share a factor of two, so their greatest common divisor is double (left-shift) the greatest common divisor of the halves (right-shift) of the two numbers.
- If one number is even and the other is odd, they have no factor of two in common, so their greatest common divisor is the greatest common divisor of the odd number and half (right-shift) the even number.
- Otherwise, both numbers are odd, and a normal Euclidean step is performed, replacing the larger of the two numbers with their difference.
Binary algorithms such as these are quite common; see for instance the ipow and expm functions in the Standard Prelude. They are especially useful for hardware implementations in the microcode of the processor, since bit-shifting operations are extremely fast.
Your task is to write functions that perform the three operations given above, simulating how they are performed in hardware. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] Praxis – Three Binary Algorithms By Remco Niemeijer In today’s Programming Praxis we have to implement binary algorithms for multiplying, dividing, and finding […]
My Haskell solution (see http://bonsaicode.wordpress.com/2010/01/15/programming-praxis-three-binary-algorithms/ for a version with comments):
import Data.Bits left, right :: Bits a => a -> a left = flip shiftL 1 right = flip shiftR 1 binmult :: (Bits a, Integral a) => a -> a -> a binmult 1 b = b binmult a b = binmult (right a) (left b) + if odd a then b else 0 bindiv :: (Bits a, Ord a) => a -> a -> (a, a) bindiv n d = f (right $ until (> n) left d) 0 n where f t q r | t < d = (q, r) | t <= r = f (right t) (left q + 1) (r - t) | otherwise = f (right t) (left q) r bingcd :: (Bits a, Integral a) => a -> a -> a bingcd a 0 = a bingcd 0 b = b bingcd a b | even a && even b = 2 * bingcd (right a) (right b) | even a = bingcd (right a) b | even b = bingcd a (right b) | a > b = bingcd (a - b) b | otherwise = bingcd a (b - a)I’ll break this down in three replies, C is just too verbose. I’m fascinated by how compact your Scheme and Haskell solutions always are… I really wanted to have a serious look at Python as a next step at learning to program, but this could make me change my mind.
Anyway, multiplication:
#define EVEN(A) (!(A & 1)) int mult(unsigned int a, unsigned int b) { int prod = 0; do { if EVEN(a) prod += b; a >>= 1; b <<= 1; } while (a >= 1); printf("Le produit est: %d\n", prod); return prod; }…sorry about the bit in french in the previous one. Division:
int div(unsigned int n, unsigned int d) { unsigned int q = 0, t = d, r = n; while (t <= n) t <<= 1; while ((t >>= 1) >= d) { if (t <= r) { q = 1 + (q << 1); r -= t; } else q <<= 1; } printf("quotient %d, remainder %d\n", q, r); return 0; }And the GCD:
#define EVEN(A) (!(A & 1)) #define MAX(A, B) (A > B ? A : B) #define MIN(A, B) (A > B ? B : A) int gcd(unsigned int m, unsigned int n) { if ((m == 0) || (n == 0)) return MAX(m, n); if (EVEN(m) && EVEN(n)) return (gcd((m >> 1), (n >> 1)) << 1); if (EVEN(m) || EVEN(n)) { //we've just tested for &&, this is now effectively an XOR unsigned int odd, even; EVEN(m) ? (even = m, odd = n) : (even = n, odd = m); return gcd(odd, (even >> 1)); } return gcd(MAX(m, n) - MIN(m, n), MIN(m, n)); }I managed to get this one spectacularly wrong initially, by underestimating the difference between preprocessor macros and functions: I’d initially left out the brackets around the expression of the MAX and MIN macros, and obviously the operators precedence was playing a nasty trick on me in the gcd(MAX(m,n) – MIN(m,n), MIN(m,n)). I’d like to think it’s a lesson I won’t forget…
Python:
def binmul(p, s): prod = 0 while p>=1: if p&1: prod += s p >>= 1 s <<= 1 return prod def bindiv(n, d): t = d while t < n: t <<= 1 t >>= 1 q, r = 0, n while t >= d: q <<= 1 if t < r: q += 1 r -= t t >>= 1 return q,r def bingcd(a, b): twos = 0 while a and b: if a&1: if b&1: a,b = (a-b,b) if a>b else (a,b-a) else: b >>= 1 elif b&1: a >>= 1 else: twos += 1 a >>= 1 b >>= 1 return (a or b) << twosWorking my way through old exercises; here’s my Python solution (note: the previously posted solution doesn’t look like it quite conquered division):
def mul(n, m): prod = 0 while n != 1: if n & 1: prod += m m <<= 1 n >>= 1 return prod + m def div(n, m): t = m while t <= n: t <<= 1 t, q, r = t >> 1, 0, n while t >= m: q <<= 1 if t <= r: q, r = q + 1, r - t t >>= 1 return q, r def gcd(n, m): if n == 0 or m == 0: return m + n elif not n & 1 and not m & 1: return gcd(n >> 1, m >> 1) << 1 elif not n & 1: return gcd(n >> 1, m) elif not m & 1: return gcd(n, m >> 1) else: return gcd(min(n, m), max(n, m) - min(n, m))[…] rozruszania szarych komórek i nabrania praktyki w programowaniu. Na początek jedno z zadań z Programming Praxis, czyli prosta implementacja algorytmu mnożenia dwóch liczb całkowitych przy pomocy przesuwania o […]
My solution in Forth:
: imul 0 >r 2dup < IF swap THEN BEGIN dup 1 and IF over r> + >r THEN swap 2* swap 2/ dup 0= UNTIL 2drop r> ; : idiv 0 >r swap >r ( r: q r) dup BEGIN 2* dup r@ > UNTIL ( d t) BEGIN 2/ 2dup <= WHILE r> r> 2* >r ( d t r) 2dup <= IF r> 1+ >r over - THEN >r REPEAT 2drop r> r> ; : 2sort 2dup > IF swap THEN ; : max 2sort nip ; : even? 1 and 0= ; : 2even? even? swap even? and ; : igcd 0 >r BEGIN 2dup * 0= IF max r> lshift EXIT THEN 2dup 2even? IF r> 1+ >r 2/ swap 2/ swap ELSE dup even? IF 2/ ELSE over even? IF swap 2/ swap ELSE 2sort over - THEN THEN THEN AGAIN ;[…] suggested it would make a good exercise. Indeed it would; in fact, we have already done it twice ([1] […]
[…] X. I describe (with code) the original Euclidean algorithm, the modern Euclidean algorithm, the binary algorithm, and the extended Euclidean algorithm at my […]