On Hacker News, user jimmyjim comments on the recent GoDaddy-related outrage,
In the words of Elie Wiesel, "Take sides. Neutrality helps the oppressor. Never the victim. Never the tormented." ... Too often people think inaction implies neutrality. It does not. It's an enabling behavior. Passivity is a free permission slip for the status quo to persist onward.
Bei Kontext TV gibt es ein interessantes Interview mit Immanuel Wallerstein (deutsch, englisch) zu den „Grenzen des Kapitalismus“ und sein absehbares Ende.
Auch sehenswert ist der Beitrag zum „Wachstum bis zum Kollaps“ (deutsch, englisch) mit Vandana Shiva und Tim Jackson.
While examining Racket's built-in, continuation-based web application library, I noticed that by default, serialized continuations are sent to the HTTP client both unencrypted and unsigned. After a bit of browsing the manual, I discovered that there is some built-in support for signing serialized data using HMAC-SHA1. However, HMAC-SHA1, while not as insecure as a naïve approach using plain SHA1, isn't the cryptographically strongest algorithm around. In addition, Racket's implementation requires the OpenSSL Crypto library to be present on the target system. This encouraged me to write my own little library of cryptographic algorithms.
In addition to HMAC and Whirlpool, which were originally intended to be used as a drop-in replacement for Racket's HMAC-SHA1 function, this library provides implementations of Dan Bernstein's CubeHash hashing function and Salsa20 stream cipher.
You can fetch the Mercurial repository using the following command:
hg clone http://matthias.benkard.de/code/mulkrypt-for-racket/
(Mind the trailing slash. My web server is overly picky there.)
Using the hashing functions (whirlpool, cubehash-{128,160,224,256,384,512,512x}) is pretty straight-forward. They expect a byte string as input, and output an integer representing the hash value.
There is also a cubehash constructor that can be used to create instances of CubeHash with different parameters. These are, in order: the number of initialization rounds (usually 16), the number of rounds per block (usually 16), the block size in bytes (usually 32 except for CubeHash-512x), the number of finishing rounds (usually 32), and the number of output bits.
Using the Salsa20 stream cipher is also pretty simple: The expected inputs are, in order, the encryption key (which must be either 16 or 32 bytes long), the nonce (8 bytes), and a sequence of input bytes. The output is a sequence of bytes containing the ciphertext. The input sequence may be infinite. In this case, the output sequence will also be infinite. Otherwise, the output sequence has the same length as the input sequence.
HMAC-based message authentication is a bit trickier than either encryption or hashing, the reason being that HMAC needs to know about the block size the hashing function uses. The required inputs are, in order: the hashing function (such as whirlpool), the hashing function's block size in bytes (64 for Whirlpool, usually 32 for CubeHash), the size of the hash value in bytes, the authentication key, and finally, the message to authenticate.
The result is the computed authentication code that needs to be transmitted along with the message.
> (define fox-and-dog #"The quick brown fox jumps over the lazy dog")
> (format "~x" (cubehash-512x fox-and-dog))
"a0e462f1cc35c81f73b6e0aad5f92a20ff5cc25d2bf607ff6db9fa47711db65d
65fd2e1df20857609480968e9189825ddb1f7dd57ac1e0d985e9a27fe5f0a5ec"
> (format "~x" (whirlpool fox-and-dog))
"b97de512e91e3828b40d2b0fdce9ceb3c4a71f9bea8d88e75c4fa854df36725f
d2b52eb6544edcacd6f8beddfea403cb55ae31f03ad62a5ef54e42ee82c3fb35"
> (format "~x" (hmac cubehash-128 32 16 #"key" fox-and-dog))
"f17f630db602ad132361129999399d5"
> (sequence->list (salsa20 #"<<<<<key of 16 or 32 bytes.>>>>>"
#"<nonce.>"
fox-and-dog))
(122 146 89 152 145 152 42 248 229 225 154 77 160 250 22 247 6 191 4 202 154
186 141 181 124 241 166 3 1 133 42 128 25 136 203 216 165 24 195 21 159 96 105)
Some time ago, I wrote an implementation of an optionally functional, bitmapped Patricia Tree as part of my project thesis. The project thesis developed into something quite different from what I had imagined, and I consequently never got around to using any Patricia Trees in the actual code. Since it might be useful in a variety of situations, I hereby release the implementation as a library.
You can fetch the Mercurial repository using the following command:
hg clone http://matthias.benkard.de/code/mulklib/
(Mind the trailing slash. My web server is overly picky there.)
Build with something like the following:
make CC=gcc LIB_PREFIX=lib LIB_SUFFIX=.so
The defaults are CC=clang LIB_PREFIX=lib LIB_SUFFIX=.dylib, which is all right for Mac OS X but not so great on other systems. They are also subject to change, since they depend on my development environment. (Of course, if you feel the urge to make the build system a little smarter, patches are always welcome.)
If building and installing a library consisting of a single object file seems like overkill to you (it sure does to me), you may prefer simply integrating the source files (bitmapped_patricia_tree.c, bitmapped_patricia_tree.h) into your own project directly.
// 1. Patricia trees are very amenable to structure sharing. // // 2. Furthermore, big-endian Patricia trees are especially efficient // when indices are allocated sequentially. // // 3. Finally, bitmapping improves the performance of copying because // copying an array is much cheaper than copying an equivalent branch // in a tree.
A Bitmapped Patricia Tree (BPT) maps integer keys to void* values.
BPTs provide a functional, Lisp-like interface. All modifying operations (bpt_assoc and bpt_dissoc, which add and remove entries, respectively) return a new BPT with the requested changes applied, although by default, they may destructively modify the original BPT for improved performance.
Destructive modification can be prohibited by sealing the original BPT using the bpt_seal procedure. Note that there is practically no performance overhead in sealing (the seal being just a boolean flag), so there is nothing wrong with using bitmapped Patricia Trees in a purely functional manner by calling bpt_seal after every operation.
An empty BPT is represented as NULL.
Memory management is done through a reference counting scheme inspired by OpenStep conventions. Newly allocated objects are returned with a reference count of 1. In order to increment the reference count, call bpt_retain; in order to decrement it, bpt_release. bpt_dealloc should never be called directly, although it is part of the public interface.
Client programs can have themselves be notified when a leaf node is deallocated by setting a deallocation callback with bpt_set_dealloc_hook. The callback is passed the key and value that the deallocated node represents. This can, among other things, be used to free the memory pointed to by the value (assuming it represents a pointer).
//
// This is a very basic example of using BPTs that blissfully neglects memory
// management in favor of didactic simplicity.
//
// For a more complete example, see bpt_test.c in the source distribution.
//
#include "stdio.h"
#include "stdlib.h"
#include "bitmapped_patricia_tree.h"
void print_tree(bpt_t b) {
int i;
for (i = 0; i < 10; i++) {
if (bpt_has_key(b, i)) {
printf(" %d -> %s\n", i, bpt_get(b, i));
}
}
}
int main() {
bpt_t b1, b2;
// Create b1.
b1 = bpt_assoc(NULL, 0, "zero"); //functional (obviously)
b1 = bpt_assoc(b1, 1, "one"); //|
b1 = bpt_assoc(b1, 2, "two"); //|destructive!
b1 = bpt_assoc(b1, 3, "three"); //|
b1 = bpt_assoc(b1, 4, "four"); //|
// Make b1 functional.
bpt_seal(b1);
// Make b2 a structure-sharing copy of b1.
b2 = bpt_assoc(b1, 0, "null"); //functional, as b1 is sealed
b2 = bpt_assoc(b2, 3, "a triple"); //destructive
// Modifying b1 does not affect b2.
b1 = bpt_assoc(b1, 5, "five"); //destructive
printf("Map 1:\n");
print_tree(b1);
printf("\nMap 2:\n");
print_tree(b2);
return EXIT_SUCCESS;
}
Output:
Map 1: 0 -> zero 1 -> one 2 -> two 3 -> three 4 -> four 5 -> five Map 2: 0 -> null 1 -> one 2 -> two 3 -> a triple 4 -> four
typedef int32_t bpt_key_t; // Querying void *bpt_get(bpt_t bpt, bpt_key_t key); bool bpt_has_key(bpt_t bpt, bpt_key_t key); // Adding and Removing Entries bpt_t bpt_assoc(bpt_t bpt, bpt_key_t key, void *item); bpt_t bpt_dissoc(bpt_t bpt, bpt_key_t key); // Managing Memory void bpt_retain(bpt_t bpt); void bpt_release(bpt_t bpt); void bpt_dealloc(bpt_t bpt); void bpt_set_dealloc_hook(bpt_t bpt, bpt_key_t key, void (*hook)(bpt_key_t key, void* value)); // Making Maps Functional void bpt_seal(bpt_t bpt);
Arguably, one of the most elegant ways of representing a recursively enumerable set such as the set of prime numbers in a computer program is as an enumeration, or, as some might call it, a generator. A generator, simply speaking, is an object that repeatedly generates the next value of a set when asked for one by the user.
Generators are a sufficiently useful concept that even mainstream programming languages like Python support them with special syntax:
def squares():
i = 0
while True:
yield i*i
i += 1
sqs = squares() # Instantiate generator
print sqs.next() # => 0
print sqs.next() # => 1
print sqs.next() # => 4
print sqs.next() # => 9
Some less mainstream programming languages, such as Haskell, support the concept of a stream instead, which is basically an infinite list:
squares = [x*x | x <- [0..]] main = do println $ take 4 squares
The main difference between a generator and a stream is that a stream need not be instantiated. Whereas a generator has a state that changes every time the user requests a value, a stream always behaves as the same infinite list. As with all things, this entails both disadvantages and benefits, which I am not going to talk about here. However, one thing that makes streams particularly interesting is that since they are static, an implementation can choose to cache the values that have already been generated, so that subsequent uses of the same stream can simply use the cached values instead of generating them from the start again.
This caching behavior is particularly important when generation is expensive—such as is the case for the set of prime numbers.
Clearly, it is not particularly hard to generate all prime numbers. Simply provide a primality test and filter the natural numbers:
def is_prime(n):
for p in primes():
# Break out of the loop if there is
# no chance left of finding a factor.
if p*p > n:
return True
# If p is a factor of n, n ain't a
# prime number.
if n % p == 0:
return False
def primes():
yield 2
i = 3
while True:
if is_prime(i):
yield i
i += 1
Pretty simple, right? Note that the primes generator is used in the definition of is_prime. This works because primes is defined to yield the first prime number without calling is_prime at all (otherwise, this would lead to an infinite recursion).
Now, the above code is pretty inefficient since is_prime repeatedly generates the primes up to sqrt(n). This could easily be solved by making primes a stream, so let's try making it one. In order to do this, I shall switch to the Scheme-derived Racket programming language for the rest of the article, since it provides a rich set of features dealing with generators and streams:
;; Some helper definitions.
(define (divides? p n)
(zero? (remainder n p)))
(define (nats [n 0] [step 1])
(stream-cons n (nats (+ step n) step)))
;; This is the interesting part.
(define (prime? n)
(for/and ([p (stop-after primes (λ (p) (>= (* p p) n)))])
(not (divides? p n))))
(define primes
(stream-cons 2 (stream-filter prime? (nats 3))))
Again, prime? and primes are mutually recursive, but this time, one and the same stream of primes is reused in every invocation of prime?, making it unnecessary to generate the list of primes up to (sqrt n) every single time.
If the functional style is unfamiliar to you, it would have been quite possible to define primes just as in the Python example above, albeit with a little more syntactic overhead (since we still want to convert the generated sequence into a stream):
(define primes
(sequence->stream
(in-generator
(yield 2)
(for ([n (in-naturals 3)])
(when (prime? n)
(yield n))))))
So the above prime generation technique is nice and elegant and everything, but it's also hopelessly inefficient.
> (time (stream-ref primes 20000)) cpu time: 24822 real time: 25316 gc time: 379 224743
Surely we can do better. There are a couple of possible optimizations here. One is to test for primality only those numbers that are +/- 1 mod 6. (Figuring out the reason every prime number greater than 3 has this form is left as an exercise for the reader.) This gives us the following improved primality test:
(define (prime? n)
(or (= n 2)
(= n 3)
(let ([mod6 (modulo n 6)])
(and (> n 1)
(or (= mod6 1)
(= mod6 5))
(for/and ([p (stop-after primes (λ (p) (>= (* p p) n)))])
(not (divides? p n)))))))
This improves performance by a bit, but it still leaves a lot to be desired.
> (time (stream-ref primes 20000)) cpu time: 23223 real time: 23424 gc time: 554 224743
The reason for the sub-par performance is that we are still using a trial-division algorithm. A much more efficient way of generating primes is by using a prime sieve algorithm.
A simple prime sieve algorithm works like this (this is known as the Sieve of Eratosthenes): Assume some bound n, below which we need to enumerate all prime numbers. Allocate an array of length n. Initialize the array with ones, marking all numbers greater than 1 as potential prime numbers. Now, for every new prime number that you find (by simply looking at the first number in the list that is not yet marked), mark off all multiples by setting their array cells to zero. You may stop when you have reached sqrt(n).
Since the times of Eratosthenes, mathematicians have discovered more elaborate prime sieve algorithms. Of note in particular is the Sieve of Atkin, which is known to be significantly more efficient than the Sieve of Eratosthenes.
Unfortunately, it is not at all clear how to use a traditional prime sieve algorithm to generate an infinite stream of primes. Just like the Sieve of Eratosthenes, modern sieve algorithms assume some upper bound n which primes should be generated below of. The reason is that a sieve works by marking off all multiples of some number, and this marking-off process must end at some point. In addition, some of the efficiency benefits that sieve algorithms provide is the fact that they are able to stop relatively early. For instance, the Sieve of Eratosthenes can terminate as soon as the last prime smaller than sqrt(n) has been processed.
Certainly, you could set some preliminary bound and restart the algorithm whenever the bound is crossed, setting a new bound in the process. There may even be a way of choosing the new bound so as to avoid making the asymptotic runtime significantly worse (there's probably a paper or two in here somewhere); but even so, an approach that recomputes a lot of prime numbers again and again is still dissatisfactory.
You would think that there must be some other way—and indeed there is. Let's take a step back for a moment. What is a sieve algorithm actually about? The algorithm steps through the already filtered part of the array, and for each value that it finds, it starts to check off its multiples (or square multiples, or whatever it is the specific algorithm does) up to the end of the array. But wait—why should such a subprocess have to terminate as soon as the end of the array is reached? It might as well sleep until the array becomes larger, at which point it can continue to check off items—after which we can be sure that the enlarged array is also filtered. So if we can somehow make it do that, we have magically solved all of our problems. (Admittedly, this isn't quite true, since the sieve process itself also needs to be suspended whenever it is okay to do so—but that is a negligible issue, since we can just make it stop and continue as well.)
This is indeed feasible. In fact, there are multiple ways of doing something like this. We will take a rather simple and intuitive route: Make every check-off process a coroutine, collect the resulting coroutines into a list, and run them in order every time we resize the array. Each time, the coroutines will stop at the end of the array, waiting for the next iteration to arrive. In pseudo-code,
define sieve coroutine:
loop:
for i from 1 to infinity:
if sieve[i] is 1:
spawn check-off coroutines for i
if i >= stopping number (depending on size of array):
yield to parent
(and sleep until reentered)
define primes:
set up sieve
set coroutines = []
loop:
run sieve coroutine,
saving spawned coroutines into the coroutines list
run coroutines from coroutines list
enumerate all primes from sieve
resize sieve
We implement the above by first defining a facility for creating coroutines:
(define (make-coroutine fn)
(let ([cont #f]
[return #f])
(λ ()
(let/cc return*
(set! return return*)
(if cont
(cont)
(fn (λ () (let/cc k
(set! cont k)
(return)))))))))
Next, we need to implement the main loop. We do this by writing it directly as a stream—though we could also have used a generator, as noted above.
;; Parameters:
;; sieve-constructior — the function that constructs the sieve coroutine
;; starting-numbers — a list of numbers that the sieve does not yield
;; but that still need to be in the stream
;; index->number — a procedure that maps an array index to the number it represents
;; (in the Sieve of Eratosthenes, this is simply the identity function)
(define (sieve->stream sieve-constructor starting-numbers index->number)
;; Starting array size: 100'000 items.
(let* ([sieve (box (make-vector 100000 #t))]
[coroutines (box (list))]
[sieve-coro (sieve-constructor sieve coroutines)])
(define (iter pos)
(letrec
([len
(vector-length (unbox sieve))]
[init
;; Run sieve and filter coroutines.
(λ ()
(sieve-coro)
(for ([coro (unbox coroutines)])
(coro)))]
[loop
;; Yield the generated primes by iterating over the sieve.
(λ (i)
(if (< i len)
(begin
(if (vector-ref (unbox sieve) i)
(stream-cons (index->number i)
(loop (add1 i)))
(loop (add1 i))))
(next-iteration i)))]
[next-iteration
;; Resize array, repeat.
(λ (newpos)
(let* ([old-sieve (unbox sieve)]
[new-size (inexact->exact
(round
(* (vector-length old-sieve)
1.2)))])
(printf "Resizing prime sieve. New size: ~A\n" new-size)
(set-box! sieve (make-vector new-size))
(vector-fill! (unbox sieve) #t)
(vector-copy! (unbox sieve) 0 old-sieve)
(iter newpos)))])
(init)
(loop pos)))
(define (start-iter starting-numbers)
(match starting-numbers
[(list)
(iter 1)]
[(list* x xs)
(stream-cons x (start-iter xs))]))
(start-iter starting-numbers)))
For the sieve algorithm, we can use basically whichever one we want. I have chosen to adapt an algorithm published on StackOverflow.com by one Robert William Hanks (which he says is not his real name). The original version posted on SO is the following:
import numpy
def primesfrom2to(n):
""" Input n>=6, Returns a array of primes, 2 <= p < n """
sieve = numpy.ones(n/3 + (n%6==2), dtype=numpy.bool)
for i in xrange(1,int(n**0.5)/3+1):
if sieve[i]:
k=3*i+1|1
sieve[ k*k/3 ::2*k] = False
sieve[k*(k-2*(i&1)+4)/3::2*k] = False
return numpy.r_[2,3,((3*numpy.nonzero(sieve)[0][1:]+1)|1)]
We translate this into a coroutine-based implementation as below. Note that the sieve routine destructively modifies the sieve array. In addition, we use explicit boxes in order to communicate changes in the list of coroutines from the sieve algorithm to the main loop as well as for the array, which is reallocated on every major iteration.
;; A little helper definition.
(define-syntax-rule (while cond body ...)
(when cond
(let loop ()
body ...
(when cond
(loop)))))
(define (make-rwh-prime-sieve2 sieve coroutines)
(make-coroutine
(λ (yield)
(let loop ([pos 1])
(let* ([n (add1 (* (vector-length (unbox sieve)) 3))])
(for ([i (in-range pos (add1 (quotient (integer-sqrt n) 3)))])
(set! pos (add1 i))
(when (vector-ref (unbox sieve) i)
(let* ([k (bitwise-ior (+ (* 3 i) 1)
1)]
[make-filter
(λ (starting-pos)
(make-coroutine
(λ (yield)
(let loop ([ii (quotient starting-pos 3)])
(while (>= ii (vector-length (unbox sieve)))
(yield))
(vector-set! (unbox sieve) ii #f)
(loop (+ ii (* 2 k)))))))])
(set-box! coroutines
(list* (make-filter (sqr k))
(make-filter (* k
(+ k
(- (* 2 (if (odd? i) 1 0)))
4)))
(unbox coroutines)))))))
(yield)
(loop pos)))))
Finally, putting it all together, we finally get a stream of primes:
(define primes
(sieve->stream make-rwh-prime-sieve2
'(2 3)
(λ (i)
(bitwise-ior (+ (* 3 i) 1) 1))))
So after all this, have we gained anything with regard to performance? Indeed we have:
> (time (stream-ref primes 20000)) cpu time: 387 real time: 391 gc time: 148 224743
In fact, we can now easily generate one order of magnitude more prime numbers than before (note that this does not merely generate the prime numbers up to 100'000—those would be quite a bit fewer!):
> (time (stream-ref primes 100000)) cpu time: 2120 real time: 2322 gc time: 621 1299721
Not too shabby, I'd say, considering that the program has now stored every single prime up to 1'299'721 in a stream of length 100'001. We can verify that this is true by summing over all of them and measuring the time:
> (time (stream-fold + 0 (stream-take primes 100001))) cpu time: 301 real time: 301 gc time: 0 62261998442
Similarly, generating some more primes reuses the existing stream:
> (time (stream-take (stream-drop primes 100001) 10)) cpu time: 270 real time: 269 gc time: 0 '(1299743 1299763 1299791 1299811 1299817 1299821 1299827 1299833 1299841 1299853)
We can conclude that using prime sieves is an efficient way of enumerating primes. More interestingly, even though they are inherently imperative in nature, heavily leveraging side-effects and mutable arrays, they can be given an elegant functional interface using streams and can, with some work, even be used to enumerate unbounded sequences of prime numbers by making use of filtering coroutines.