Factor[ial] 102

Previously, I wrote an introduction to Factor though Factorial called Factor[ial] 101. If you saw that, you probably thought that's all you'd see of that totally un-compelling example. Well, you're wrong! We can actually implement factorial in a more simple way and more efficiently with large integers.

Let's look at the solution from last time:

: factorial ( n -- n! )
    1 [ 1+ * ] reduce ;

When you saw this, you might have been thinking, "Why not just get a list of the numbers from 1 to n and do their product? Why bother with 1+? In Haskell, it's just product [1..n]." We can actually use this strategy in Factor using the math.ranges library, which has code a word called [1,b] which creates a virtual sequence containing the integers 1 through its argument. We can also use the word product in math.vectors to get the product of a sequence. So the new factorial is just

: factorial ( n -- n! )
    [1,b] product ;

Efficiency and bignums

I previously talked about how to make some simple mathematical functions work well with bignums. The example I used there was a string>number conversion procedure, but it applies equally to getting the product of a list. In short: when multiplying two bignums of size (in bits) n and m, we know there's a lower bound of Ω(nm), since a bignum of nm bits must be constructed. So if we go about finding the product of a sequence by starting with 1, then multiplying that by the first element, then the second, and so on for the entire sequence, then this will take Ω(n²) time where n is the size of the resulting product!

We can do better, and a better strategy is to use binary recursion similar to mergesort: split the sequence in half, find the products of both halves, and multiply them. Then the easy lower bound is lower: Ω(n log n). (Note: the upper bounds for these algorithms is something like the lower bound times log n, with a reasonably efficient multiplication algorithm.)

So here's a better implementation of product, using this strategy:

: halves ( seq -- beginning end )
    dup length 2 /i cut-slice ;

: product ( seq -- num )
    dup length {
        { 0 [ drop 1 ] }
        { 1 [ first ] }
        [ drop halves [ product ] bi@ * ]
    } case ;

Abstraction and combinators

If we want to implement a word to sum an array, we'll be repeating a lot of code. So let's abstract the basic idea of this kind of recursion into a combinator, or higher order function, that we can supply the starting value and combining function to. With this, we should be able to write

: product ( seq -- num ) 1 [ * ] split-reduce ;
: sum ( seq -- num ) 0 [ + ] split-reduce ;

where split-reduce is this new combinator. Now, it's no harder to write code using the binary-recursive strategy than the original naive strategy, if split-reduce is somewhere in the library. Here's how you can implement it:

: split-reduce ( seq start quot -- value )
    pick length {
        { 0 [ drop nip ] }
        { 1 [ 2drop first ] }
        [ drop [ halves ] 2dip [ [ split-reduce ] 2curry bi@ ] keep call ]
    } case ; inline

This looks a little messy, as combinators sometimes get. Let's see how it looks using local variables: (I'm using the locals vocab, which allows the syntax :: word-name ( input variables -- output ) code... ; for lexical scoping. I usually don't use this very often, but in implementing combinators it can make things much cleaner.)

:: split-reduce ( seq start quot -- seq' )
    seq empty? [ start ] [
        seq singleton? [ seq first ]
        [ seq halves [ start quot split-reduce ] bi@ quot call ] if
    ] if ; inline

Which one of these you prefer is purely a matter of taste.

A tangent to mergesort

What if we wanted to use split-reduce to implement mergesort? It might look like you can do this:

: mergesort ( seq -- sorted )
    { } [ merge ] split-reduce ;

However, there's a problem here: in the base case, if we have { 1 }, it'll be changed into 1. But we need the base case to output sequences! (Ignore the fact that 1 is a sequence; it's of the wrong type.) So the cleanest way to do this is to make a new word, split-process, which does the same thing as split-reduce but takes a new parameter specifying what to do in the base case. With this we're able to do

: split-reduce ( seq start quot -- value )
    [ first ] split-process ; inline

To implement this, we just need to modify split-reduce, factoring out the base case code:

:: split-process ( seq start quot base-quot -- seq' )
    seq empty? [ start ] [
        seq singleton? [ seq base-quot call ] [
            seq halves
            [ start quot base-quot split-process ] bi@
            quot call
        ] if
    ] if ; inline

Now mergesort can be implemented as

: mergesort ( seq -- sorted )
    { } [ merge ] [ ] split-process ;

for some suitable implementation of merge.

To binrec and beyond

What if we took this even further: why restrict this to binary recursion on sequences? We can do binary recursion on everything that needs binary recursion! So let's make a combinator out of this, calling it binrec. binrec takes four—four!—quotations. The first one specifies the termination (base case) condition. The second specifies what to do in the base case. The third specifies how to split up the data in the inductive case, and the fourth specifies how to put the two pieces back together after the recursion takes place. Here's how we can implement binrec for a totally general binary recursion combinator:

:: binrec ( data test end split rejoin -- value )
    data test call [ data end call ] [
        data split call
        [ test end split rejoin binrec ] bi@
        rejoin call
    ] if ; inline

In the abstract, this isn't too bad. But how can you read code that uses binrec? You have to remember four quotations, their intended stack effects and their role in calculating this. For me, this is too difficult to do in most cases.

Look at how we can define split-process in terms of binrec:

:: split-process ( seq start rejoin-quot end-quot -- value )
    [ dup singleton? swap empty? or ]
    [ dup singleton? [ end-quot call ] [ drop start ] if ]
    [ halves ] rejoin-quot binrec ; inline

This isn't actually easier than defining split-process directly, and you can argue that it's worse than the original version. Still, it provides an interesting way to avoid explicit recursion.

Pulling it all together

Complicated combinators like binrec can be useful, sometimes, as long as you don't use them directly. One of the great things about Factor is that it's so easy to specialize these things. So why not? Almost every case that you're using binrec follows a particular pattern.

We can tell everyone more loudly about split-reduce, which is much easier to use, and have binrec be hidden in the library for advanced users who want to implement their own similar combinators without repeating the code that's already written in binrec. It's not that recursion is difficult to do, its just that there's no reason to write this code more than once.

So that's how you implement factorial in Factor. Except once all this is in the library, all you have to worry about is [1,b] product.

(BTW If you actually want to use factorial for something practical, where it'll be called multiple times, a memoizing table-based approach might be faster. Or maybe Stirling's approximation is appropriate, depending on the use case. Or maybe one of these algorithms. But that's a topic for another blog post!)

Update: I fixed some typos, found by Slava. Also, Slava added split-reduce into the core as binary-reduce and implemented sum and product with it.

Update 2: Updated the code samples for the current version of Factor, as of late March 2009.