Re: Factor

I love reading John D. Cook’s blog which covers various mathematical concepts, code snippets, and other curious explorations. Today, he celebrated his 5,000th post which is a pretty incredible milestone. In honor of that, let’s implement the Recamán’s sequence which he wrote about yesterday:

I recently ran into Recamán’s sequence. N. J. A. Sloane, the founder of the Online Encyclopedia of Integer Sequences calls Recamán’s sequence one of his favorites. It’s sequence A005132 in the OEIS.

This sequence starts at 0 and the nth number in the sequence is the result of moving forward or backward n steps from the previous number. You are allowed to move backward if the result is positive and a number you haven’t already visited. Otherwise you move forward.

Here’s Python code to generate the first N elements of the sequence.

def recaman(N):
    a = [0]*N
    for n in range(1, N):
        proposal = a[n-1] - n
        if proposal > 0 and proposal not in set(a):
            a[n] = proposal
        else:
            a[n] = a[n-1] + n
    return a

For example, recaman(10) returns

[0, 1, 3, 6, 2, 7, 13, 20, 12, 21]

Direct Implementation

The code for this subtract if you can, add if you can’t sequence is not particularly complex, but it is often fun to see what an algorithm looks like when implemented in a different language. So, we’ll try a direct implementation in Factor first and then talk about some variations.

Using local variables, we can keep our version very similar to the original:

:: recaman ( N -- seq )
    N 0 <array> :> a
    N [1..b) [| n |
        n 1 - a nth n - :> proposal
        proposal 0 > proposal a member? not and [
            proposal n a set-nth
        ] [
            n 1 - a nth n + n a set-nth
        ] if
    ] each a ;

And, show that it works:

IN: scratchpad 10 recaman .
{ 0 1 3 6 2 7 13 20 12 21 }

Short-Circuit Logic

One difference that is subtle is that Python’s and operation is short-circuiting so that if the first expression returns False, the second expression is not evaluated. We can make that change ourselves:

- proposal 0 > proposal a member? not and [
+ proposal 0 > [ proposal a member? not ] [ f ] if [

Or use short-circuit combinators to do that for us:

- proposal 0 > proposal a member? not and [
+ proposal { [ 0 > ] [ a member? not ] } 1&& [

Reduce Repetitions

We could also notice repetitions that occur, for example both branches of the if set the nth value of the array. In addition, we could keep the proposal on the stack instead of naming it as a local variable, and then just replace it when the branch evaluates to f.

Incorporating these changes results in this version:

:: recaman ( N -- seq )
    N 0 <array> :> a
    N [1..b) [| n |
        n 1 - a nth n -
        dup { [ 0 > ] [ a member? not ] } 1&& [
            drop n 1 - a nth n +
        ] unless n a set-nth
    ] each a ;

There is still repetition as each iteration looks up the previous value twice. Instead, we could fix that by storing that value as a local variable:

:: recaman ( N -- seq )
    N 0 <array> :> a
    N [1..b) [| n |
        n 1 - a nth :> prev
        prev n - dup { [ 0 > ] [ a member? not ] } 1&&
        [ drop prev n + ] unless n a set-nth
    ] each a ;

Sometimes it is simpler to extract an inner loop and then see how it looks:

:: next-recaman ( prev n a -- next )
    prev n - dup { [ 0 > ] [ a member? not ] } 1&&
    [ drop prev n + ] unless ;

:: recaman ( N -- seq )
    N 0 <array> :> a
    N [1..b) [
        [ 1 - a nth ] [ a next-recaman ] [ a set-nth ] tri
    ] each a ;

The retrieval of the previous value is un-necessary work on each iteration. We could remove that by keeping the previous value on the stack and use the make vocabulary to accumulate values. We also start from index 0 instead of 1 as a simplification:

:: next-recaman ( prev n -- next )
    prev n - dup { [ 0 > ] [ building get member? not ] } 1&&
    [ drop prev n + ] unless ;

:: recaman ( N -- seq )
    [ 0 N <iota> [ next-recaman dup , ] each drop ] { } make ;

Constant-time Lookups

But, there is a big performance issue. Searching the previous values to see if it had been generated before uses a linear-time algorithm. It would be much faster to use a bit-vector to track the previously seen numbers and a constant-time lookup:

:: next-recaman ( prev n seen -- next )
    prev n - dup { [ 0 > ] [ seen nth not ] } 1&&
    [ drop prev n + ] unless t over seen set-nth ;

:: recaman ( N -- seq )
    0 N dup <bit-vector> '[ _ next-recaman dup ] map-integers nip ;

This allows us to generate, for example, one million values in the sequence in a small amount of time:

IN: scratchpad [ 1,000,000 recaman last . ] time
1057164

Running time: 0.049441708 seconds

Or even calculate the nth value without saving the sequence:

:: nth-recaman ( N -- elt )
    0 N dup <bit-vector> '[ _ next-recaman ] each-integer ;

Using Generators

We could even use the generators vocabulary to make an infinite stream that we can sample from:

:: next-recaman ( prev n seen -- next )
    prev n - dup { [ 0 > ] [ seen nth not ] } 1&&
    [ drop prev n + ] unless t over seen set-nth ;

GEN: recaman ( -- gen )
    0 0 ?V{ } clone '[
        [ _ next-recaman dup yield ] [ 1 + ] bi t
    ] loop 2drop ;

And then take some values from it:

IN: scratchpad recaman 10 take .
{ 0 1 3 6 2 7 13 20 12 21 }

It is often pleasantly surprising to explore seemingly simple topics – for example, looking at this algorithm and thinking about different time complexity, data structures, or programming language libraries that might be useful.

I have added this today to the Factor development version.