Monday, July 12, 2010

Curious? Factor to the rescue!

A lot of dynamic language users have strong interest in math (particularly those who use functional programming languages). Some months back I came across a fun blog post about using Clojure to approximate Euler's Number, according to this observation:

"Here is an example of e turning up unexpectedly. Select a random number between 0 and 1. Now select another and add it to the first. Keep doing this, piling on random numbers. How many random numbers, on average, do you need to make the total greater than 1? Answer: 2.71828..."

I wanted to contrast his solution with one using Factor. We'll use these vocabularies:

USING: kernel math random sequences ;

We can use the random vocabulary to create a random float between 0 and 1, counting how many numbers we need to add before the value is greater than 1:

: numbers-added ( -- n )
    0 0 [ dup 1 < ] [
        [ 1 + ] dip 0.0 1.0 uniform-random-float +
    ] while drop ;

Taking the average of a sequence of numbers is pretty easy:

: average ( seq -- n )
    [ sum ] [ length ] bi / ;

Similar to the original blog post, we'll want to estimate the value of e by running numbers-added a certain number of times and then taking the average of the resulting sequence:

: approximate-e ( n -- approx )
    [ numbers-added ] replicate average ;

Running this from one thousand to one million times gives us an increasingly accurate approximation of e:

( scratchpad ) { 1000 10000 100000 1000000 }
               [ approximate-e >float ] map .
{ 2.694 2.7186 2.72149 2.716958 }

( scratchpad ) USING: math.constants ; e .
2.718281828459045

This could be improved by virtualizing the sequences and then performing an incremental average as we generate each of n approximations. That would likely result in both memory and speed improvements (although I haven't tried it yet).

The code for this is on my Github.

2 comments:

Ruben Berenguel said...

Just for fun I wrote it in... Clojure! My first Clojure program aside from the hello world test.

I wanted to do it in Common Lisp (or Forth, but Factor is too similar, better something new), but when I fired slime under emacs, I realised I had it configured for Clojure. I'll post it in my blog tomorrow, probably.

Thanks for sharing this, I enjoyed the read and enjoyed the coding.

Ruben

Ruben Berenguel said...

Now I see the original implementation was already in Clojure. Oh, well. Mine is more lispy :)