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 Ω(n2) 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 beyondWhat 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.
No comments:
Post a Comment