Problem 50: Which prime, below one-million, can be written as the sum of the most consecutive primes?

Alright, we completed through problem 50!

First, I thought it was about time I implemented a decent prime sequence generator instead of continuing to combine Seq.initInfinite and isPrime. It was important to me that we create an infinite sequence of primes, rather than the typical finite prime sieve. I spent a considerable amount of effort playing with various algorithms and implementations of those algorithms. The final algorithm I settled on isn’t quite a sieve, but works on the general principle, allowing us to generate the millionth prime in about 4 seconds. Compare that to this infinite, true prime sieve which can generate the millionth prime in just over a second, but I prefer to stick with my algorithm, since I worked the core of it out independently. It may be worth pointing out that the memory used in my implementation is O(n), where n is the count of the prime, while any true sieve is O(p), where p is the prime itself.

To help us out, we finally get around to implementing an extension to the Seq module, infiniteRange (note that by marking the function inline, the compiler will automatically give a structural type signature to this function so that it can be used with int32, int64, bigint, float, etc.):

let inline infiniteRange start skip = 
    seq {
        let n = ref start
        while true do
            yield n.contents
            n.contents <- n.contents + skip
    }

Next now we show our our primes implementation which is structurally generic as usual with integral type specific versions:

///p is "prime", s=p*p, c is "multiplier", m=c*p
type SievePrime<'a> = {mutable c:'a ; p:'a ; mutable m:'a ; s:'a}

///A cached, infinite sequence of primes
let inline primes_of (g:G<'a>) =
    let primeList = ResizeArray<_>()
    primeList.Add({c=g.three ; p=g.three ; m=g.three*g.three ; s=g.three*g.three})

    //test whether n is composite, if not add it to the primeList and return false
    let isComposite n = 
        let rec loop i = 
            let sp = primeList.[i]
            while sp.m < n do
                sp.c <- sp.c+g.one
                sp.m <- sp.c*sp.p

            if sp.m = n then true
            elif i = (primeList.Count-1) || sp.s > n then
                 primeList.Add({c=n ; p=n ; m=n*n ; s=n*n})
                 false
            else loop (i+1)
        loop 0

    seq { 
        yield g.two ; yield g.three

        for i in 1..primeList.Count-1 do
            yield primeList.[i].p

        yield! Seq.infiniteRange (primeList.[primeList.Count-1].p + g.two) g.two 
               |> Seq.filter (isComposite>>not)
    }

let primes = primes_of gn
let primesL = primes_of gL
let primesI = primes_of gI

Now for problem 50. This was a pretty straight forward problem. We use an array in our algorithm but do not perform any mutation:

let problem50b =
    let primesArr = primes |> Seq.takeWhile ((>=)1000000) |> Seq.toArray
    let count p =
        let rec count p cur curMax =
            let start = primesArr.[cur]
            if start >= p/2 then curMax
            else
                let rec consecutiveSum i sum =
                    if sum >= p then (i,sum)
                    else consecutiveSum (i+1) (sum + primesArr.[i])

                let i,sum = consecutiveSum cur 0            
                let curMax = max curMax (if sum = p then i-cur else 0)
                count p (cur+1) curMax
        count p 0 0

    primesArr
    |> Seq.maxBy count