Fountain Codes

September 4, 2012

We begin with the calculation of the degree distribution:

(define (degree-distribution n)
  (let* ((len (* n (+ n 1) 1/2))
         (k (randint len)))
    (let loop ((i 0) (s 0) (t 1))
      (if (< k s) (- n i -1)
        (loop (+ i 1) (+ s t) (+ t 1))))))

A fountain is a function that, each time it is called, returns a packet consisting of an xor’ed combination of blocks in its car and the list of xor’ed blocks in its cdr; it calls the degree-distribution function to determine the number of blocks in the list, then calls shuffle to generate the list of blocks without duplicates:

(define (make-fountain str)
  (let* ((cv (list->vector (map char->integer (string->list str))))
         (len (vector-length cv)))
    (lambda ()
      (let* ((d (degree-distribution len))
             (blocks (take d (shuffle (range len)))))
        (let loop ((bs blocks) (x 0))
          (if (null? bs) (cons x blocks)
            (loop (cdr bs) (logxor x (vector-ref cv (car bs))))))))))

For example, here is the function that makes a fountain function from the input string “Programming Praxis:”

(define fountain (make-fountain "Programming Praxis"))

Decoding a fountain isn’t hard if you remember the steps, though it is a bit messy; we employ three auxiliary functions to cope with the mess. The cv vector keeps track of the current state of the blocks; a known block has a value, an unknown block is #f. Store stores a singleton packet (a packet with a value and a single block) in the cv array, then returns #t for a new letter and #f for a letter that is already known. Reduce applies the known blocks to a packet; it returns (0) if the packet provides no new information because all of the blocks in the packet are known, or a packet with smaller degree than the original if some but not all of the blocks in the packet are known, or the original packet unchanged if none of the blocks in the packet are known. Ripple takes a hold, applies reduce to each packet in the hold, and produces a new hold; during a ripple operation, if any new blocks become known, ripple is called recursively, repeating until a ripple comes up empty. Overall, the decode function keeps catching drops from the fountain until there are no unknown blocks; if a packet is a singleton, it calls ripple if a new block has been identified, otherwise it adds the new packet to the hold and gets another packet from the fountain. Here’s decode:

(define (decode fountain len)
  (let* ((cv (make-vector len #f)) (n len))
    (define (store packet)
      (cond ((vector-ref cv (cadr packet)) #f)
            (else (vector-set! cv (cadr packet) (car packet))
                  (set! n (- n 1)) #t)))
    (define (reduce packet)
      (let loop ((x (car packet)) (bs (cdr packet)) (zs (list)))
        (cond ((null? bs) (cons x zs))
              ((vector-ref cv (car bs))
                (loop (logxor x (vector-ref cv (car bs)))
                      (cdr bs) zs))
              (else (loop x (cdr bs) (cons (car bs) zs))))))
    (define (ripple hold)
      (let ((done? #t))
        (let loop ((hold hold) (out-hold (list)))
          (if (null? hold)
              (if done? out-hold (ripple out-hold))
              (let ((p (reduce (car hold))))
                (cond ((equal? p (list 0))
                        (loop (cdr hold) out-hold))
                      ((null? (cddr p))
                        (if (store p) (set! done? #f))
                        (loop (cdr hold) out-hold))
                      ((equal? (cdar hold) (reverse (cdr p)))
                        (loop (cdr hold) (cons p out-hold)))
                      (else (set! done? #f)
                            (loop (cdr hold) (cons p out-hold)))))))))
    (let loop ((hold (list)))
      (if (zero? n)
          (list->string (map integer->char (vector->list cv)))
          (let ((p (fountain)))
            (if (null? (cddr p))
                (if (store p)
                    (loop (ripple hold))
                    (loop hold))
                (loop (cons p hold))))))))

Here’s an example:

> (decode fountain 18)
"Programming Praxis"

We used take, range, logxor, randint and shuffle from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/TvJutmlJ.

About these ads

Pages: 1 2

One Response to “Fountain Codes”

  1. Paul Hofstra said

    A python version. I had a lot of fun with this exercise.

    txt = "Programming Praxis Part rocks!"
    
    import random
    from operator import xor
    
    def get_boxes(d, n, prng):
        L = range(n)
        prng.shuffle(L)
        return sorted(L[:d])
            
    def get_nr_packets(n, r):
        for i in xrange(n):
            r -= n - i
            if r < 0:
                return i + 1
        
    def make_fountain(source):
        """ a generator function that emits the fountain """
        prng = random.Random()
        n = len(source)
        yield n
        source = [ord(c) for c in source]
        den = n * (n + 1) / 2
        r = prng.randint(0, 2 ** 32- 1)
        prng.seed(r)
        yield r # the random seed
        while 1:
            d = get_nr_packets(n, prng.randint(0, den-1))
            boxes = get_boxes(d, n, prng)
            yield reduce(xor, [source[b] for b in boxes])
            
    def decode_fountain(datagen):
        """ """
        prng = random.Random()
        n = datagen.next()
        prng.seed(datagen.next())
        den = n * (n + 1) / 2
        park = {}
        solved = {}
        npackets = 0
        for info in datagen:
            npackets += 1
            d = get_nr_packets(n, prng.randint(0, den-1))
            boxes = get_boxes(d, n, prng)
            s = set(boxes).intersection(solved.keys())
            for si in s:
                boxes.remove(si)
                info ^= solved[si]
            tobeprocessed = []
            if len(boxes) > 1:
                park[tuple(boxes)] = info
            elif len(boxes) == 1:
                tobeprocessed.append((boxes[0], info))
                while tobeprocessed:
                    box, info = tobeprocessed.pop()
                    if box not in solved:
                        solved[box] = info
                        for b2, info2 in park.items():
                            if box in b2:
                                del park[b2]
                                b2 = list(b2)
                                b2.remove(box)
                                info2 ^= info
                                if len(b2) == 1:
                                    tobeprocessed.append((b2[0], info2))
                                else:
                                    park[tuple(b2)] = info2
            if len(solved) == n:
                L = solved.items()
                L.sort()
                #print "number of packets", npackets, n
                return "".join([chr(z) for i, z in L])
            
    print decode_fountain(make_fountain(txt))
    

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 575 other followers

%d bloggers like this: