I'VE GOT THE BYTE ON MY SIDE

57005 or alive

Nested looping to programmatic depth in F#

Dec 26, 2014 F#

This post is the December 26th entry in the 2014 English F# Blog Advent Calendar. Check out the other great posts from earlier in the month, or head over to the Japanese-language edition!

Ho ho ho! It’s the day after Christmas, and hopefully everyone has some snazzy new toys to play with.  But just in case Santa left you a lump of coal (“… little Billy used mutable state this year …”), don’t pout, don’t cry, I’ve got you covered.  My gift for all the naughty F# developers out there: a mini-framework for imperative-style looping in F#, adding the functionality of “break,” “continue,” and “return,” along with the ability to nest loops down to programmatic depth.

The lack of “break” and “continue” in F# loops is a well-known limitation of the language. (random anecdote: the first and only time I met Miguel de Icaza, at the 2014 Microsoft MVP summit, he mentioned this as one of the things that tripped him up during his first F# dabblings). Tomas Petricek wrote up an excellent pair of posts on this topic way back in 2009, explaining how to implement the behavior of these keywords oneself by using computation expressions.

Unaware of Tomas’s approach until just this week, I re-invented this particular wheel recently. As it turns out, my approach is completely different than his, with some advantages, and some disadvantages.

Prime pairs

Let’s look at a motivating example. Say we want to find the first \(m > 0\) pairs of prime numbers that add up to some number \(n\) . (Goldbach’s Conjecture posits that at least one such pair exists for all even numbers 4 or greater. Such a pair does not exist for all odd numbers.)

An approach in C# might look something like below, assuming I already have an “IsPrime” function. This could be optimized a bit more, but it does take advantage of various shortcuts already.

static void PrimePairs(int n, int m)
{
    var count = 0;

    for (int i = 2; i <= n / 2; i++)
    {
        if (!IsPrime(i))
            continue;

        var j = n - i;
        if (IsPrime(j))
        {
            Console.WriteLine("{0} + {1}", i, j);
            if (++count >= m)
                return;
        }
    }
}

Attempting to transcribe this to F#, we immediately have trouble due to the lack of “continue” or “return”.  We can get something kind of close while still using “for” loops, but it’s fairly ugly and not as well optimized.

let primePairs n m =
    let mutable count = 0
    for i = 2 to n/2 do
        if (count >= m) || not (isPrime i) then
            ()
        else
        let j = n - i
        if isPrime j then
            printfn "%d + %d" i j
            count <- count + 1

We could, if needed, refactor the loop into a recursive function in order to get back control flow that is fully equivalent to the C# example. But that’s not always easy or elegant in practice.

Using the “imperative” computation expression from Tomas’s blog post, we gain back proper “continue” and “return” control flow, and the syntax looks much nicer:

let primePairsImperative n m =
    let count = ref 0
    imperative {
        for i = 2 to n/2 do
            if not (isPrime i) then
                do! continue
            let j = n - i
            if isPrime j then
                printfn "%d + %d" i j
                count := !count + 1
                if !count >= m then
                    return ()
    }

Alright, this is a nice start to build on.

Prime triples

Now we modify the problem so that we are to find the first \(m\) triples of prime numbers that add up to some number \(n\) .

Alright, that’s not too hard, given our existing code.  Just add a nested loop and clean up some of the boundary values.  (We can also refactor so that instead of printing values, our function returns them in a list)

let primeTriplesImperative n m =
    let count = ref 0
    let result = ref []
    imperative {
        for i = 2 to n/3 do
            if not (isPrime i) then
                do! continue
            for j = i to (n - i)/2 do
                if not (isPrime j) then
                    do! continue
                let k = n - i - j
                if isPrime k then
                    result := [i; j; k] :: !result
                    count := !count + 1
                    if !count >= m then
                        return ()
    }
    result

 Prime k-tuples

We can keep going with this, extending the code to find prime 4-tuples, 5-tuples, etc, by adding more and more nested loops. But the code quickly gets gigantic, and each implementation works only for a single, fixed tuple size.

What if we want to find prime \(k\) -tuples that sum to \(n\) , where \(k\) is a runtime parameter?  To maintain this algorithmic pattern, we would need to be nesting to some programmatic depth.

How do we do that?

Solution

I recently ran into this same general problem (though totally unrelated to prime numbers) when working on a project.  My initial solution was a bunch of bespoke code that really only applied to the specific problem I was trying to solve.  Later on, when I needed to refactor and add some new capabilities, it was, predictably, rather difficult.  So I stepped back and thought a bit what is needed to abstract “for” loops in a general way such that they can be nested together dynamically.

The resulting mini-framework I came up with can be found in a gist here.  Even heavily-commented, it’s less than 100 lines of code, but that’s still probably too big to include inline here. Take a look at the gist, and feel free to leave a comment.

A very simple example of using this to do some nested looping is below.  This implements a triply-nested loop (Depth = 3), with each loop running from 0 to 1 (Min = 0; Max = 1), printing the 3 indexes in the innermost loop.

nestedFor {
    Depth = 3; InitialRange = { Min=0; Max=1 }
    GetNextRange = fun _ -> { Min=0; Max=1 }
    Body = function
        | { Depth=3; Indexes=idxs } ->
            printfn "%A" (List.rev idxs)
            Continue
        | _ -> Descend
}

// result:
// [0; 0; 0]
// [0; 0; 1]
// [0; 1; 0]
// [0; 1; 1]
// [1; 0; 0]
// [1; 0; 1]
// [1; 1; 0]
// [1; 1; 1]

Some things to note about this framework

Using this toolset, we can re-write the prime-triples function like this:

let primeTriplesNested n m =
    let count = ref 0
    let result = ref []
    nestedFor {
        Depth = 2; InitialRange = { Min = 2; Max = n/3 }
        GetNextRange = fun { Prev = (prev :: _ ) } ->
            { Min=prev; Max = (n - prev)/2 }
        Body = function
            | { Indexes = [i] } ->
                if not (isPrime i) then Continue else Descend
            | { Indexes = [j; i] } ->
                if not (isPrime j) then Continue else
                let k = n - i - j
                if isPrime k then
                    result := [i; j; k] :: !result
                    count := !count + 1
                    if !count >= m then
                        Return
                    else Continue
                else Continue
        }
    result

The syntax is definitely clunkier than the “imperative” computation expression, but the performance is significantly better (about 13x on my laptop).

> let _ = primeTriplesImperative 10000 1000 ;;
// Real: 00:00:00.880, CPU: 00:00:00.890, GC gen0: 75, gen1: 1, gen2: 1

> let _ = primeTriplesNested 10000 1000 ;;
// Real: 00:00:00.066, CPU: 00:00:00.078, GC gen0: 21, gen1: 0, gen2: 0

And now we have the tools we need to efficiently generate \(k\) -tuples.

We notice that our nested looping algorithm basically only has 2 cases, no matter how far down we have to nest the loops:

  1. We are not yet at the final depth. All we need to do is check if the new index is a prime. If it is, then descend down to the next loop level. If it isn’t, continue to the next index at this level.
  2. We are at the final level. Check if the latest index is prime, just like before,  continuing if not. If it is, compute the final number and attempt to complete the sum.

This can be implemented like so:

let primeKTuplesNested n m k =
    let count = ref 0
    let result = ref []
    let finalDepth = k - 1
    nestedFor {
        Depth = finalDepth; InitialRange = { Min = 2; Max = n/k }
        GetNextRange = fun { Prev = (prev :: _ ) as p } ->
            { Min = prev; Max = (n - (List.sum p))/2 }
        Body = function
            | { Depth = d; Indexes = i :: _ } when d < finalDepth ->
                if not (isPrime i) then Continue else Descend
            | { Depth = d; Indexes = i :: prev } when d = finalDepth ->
                if not (isPrime i) then Continue else
                let j = n - i - (List.sum prev)
                if isPrime j then
                    result := (j :: i :: prev) :: !result
                    count := !count + 1
                    if !count >= m then
                        Return
                    else Continue
                else Continue
        }
    result

Even this fully-general solution performs pretty well for the 3-tuple problem (though not as well as the dedicated version):

> let _ = primeKTuplesNested 10000 1000 3 ;;
// Real: 00:00:00.119, CPU: 00:00:00.109, GC gen0: 22, gen1: 0, gen2: 0

And now we can create prime sums as large as we like!

> primeKTuplesNested 10000 20 5 ;;
// val it : int list list =
//   [[9221; 773; 2; 2; 2]; [9293; 701; 2; 2; 2]; [9311; 683; 2; 2; 2];
//    [9341; 653; 2; 2; 2]; [9377; 617; 2; 2; 2]; [9431; 563; 2; 2; 2];
//    [9437; 557; 2; 2; 2]; [9473; 521; 2; 2; 2]; [9491; 503; 2; 2; 2];
//    [9533; 461; 2; 2; 2]; [9551; 443; 2; 2; 2]; [9677; 317; 2; 2; 2];
//    [9743; 251; 2; 2; 2]; [9767; 227; 2; 2; 2]; [9803; 191; 2; 2; 2];
//    [9851; 143; 2; 2; 2]; [9857; 137; 2; 2; 2]; [9887; 107; 2; 2; 2];
//    [9923; 71; 2; 2; 2]; [9941; 53; 2; 2; 2]]

> primeKTuplesNested 10000 1000 10 |> List.take 5 ;;
// val it : int list list =
//   [[5381; 4591; 13; 3; 2; 2; 2; 2; 2; 2];
//    [5449; 4523; 13; 3; 2; 2; 2; 2; 2; 2];
//    [5479; 4493; 13; 3; 2; 2; 2; 2; 2; 2];
//    [5483; 4489; 13; 3; 2; 2; 2; 2; 2; 2];
//    [5521; 4451; 13; 3; 2; 2; 2; 2; 2; 2]]

This approach turned out to be very useful for me, and helped me solve my particular problem.  Hopefully it’s useful for others, as well.  Surely, there are other known approaches for attacking these kinds of algorithmic challenges, but it was a lot of fun solving it for myself this way!