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.
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))