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.

About these ads

Pages: 1 2

8 Responses to “Three Binary Algorithms”

  1. […] Praxis – Three Binary Algorithms By Remco Niemeijer In today’s Programming Praxis we have to implement binary algorithms for multiplying, dividing, and finding […]

  2. Remco Niemeijer said

    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)
    
  3. Jebb said

    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;
    }
  4. Jebb said

    …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;
    }
  5. Jebb said

    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…

  6. Mike said

    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) << twos
    
    
  7. Graham said

    Working 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))
    
  8. […] 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 […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 637 other followers

%d bloggers like this: