Three Binary Algorithms
January 15, 2010
We begin with the operations that are built-in to most hardware. R5RS Scheme does not provide bit operations, but most implementations provide an arithmetic-shift function, sometimes called ash as in Lisp. R6RS standardizes the arithmetic-shift function. Simple versions that will work anywhere are given below; for speed, they should be replaced by native operations, and in-lined.
(define (lshift x) (* x 2))
(define (rshift x) (quotient x 2))
(define (add1 x) (+ x 1))
For multiplication, we find it easier to accumulate the product as we go rather than to build the tableau of the ancient Egyptian algorithm:
(define (multiply left right)
(let loop ((left left) (right right) (prod 0))
(if (zero? left) prod
(loop (rshift left) (lshift right)
(if (odd? left) (+ prod right) prod)))))
For division, the outer loop accumulates d · 2k in a temporary variable t, then the inner loop performs the division, stopping when t has been reduced beyond the original denominator, returning both quotient and remainder:
(define (divide n d)
(let loop ((t d))
(if (<= t n) (loop (lshift t))
(let loop ((t (rshift t)) (q 0) (r n))
(cond ((< t d) (values q r))
((<= t r) (loop (rshift t) (add1 (lshift q)) (- r t)))
(else (loop (rshift t) (lshift q) r)))))))
The greatest common divisor function is a literal translation of Stein’s algorithm:
(define (stein-gcd n m)
(let loop ((n n) (m m))
(if (zero? n) m (if (zero? m) n
(if (even? n)
(if (even? m)
(lshift (loop (rshift n) (rshift m)))
(loop (rshift n) m))
(if (even? m)
(loop n (rshift m))
(if (< n m)
(loop n (- m n))
(loop (- n m) m))))))))
Here are some samples:
> (multiply 14 12)
168
> (divide 837 43)
19
20
> (stein-gcd 2322 654)
6
You can run the code at http://programmingpraxis.codepad.org/P9bxsjXo.
[…] 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 ;