Prime Time in Haskell

In a recent blog post, Patrick Vennebush of Math Jokes 4 Mathy Folks noted that 2011 can be expressed as a sum of consecutive prime numbers, and challenged his readers to work out how. He also posed a couple further challenges:

As it turns out, 2011 is extra cool because it can be written as a sum of a prime number of consecutive prime numbers.

When will that happen again?

And finally…

What is the next year that will be a prime number and also be a sum of a prime number of consecutive prime numbers? (Wow, that’s a mouthful, ain’t it?)

Just for fun I’d like to share how I computed solutions to Patrick’s challenges in my favorite programming language, Haskell. It’s great for doing this sort of recreational mathematics (among many other things), as I hope you will see. Don’t worry if you don’t know it, I’ve tried to write this post in such a way that anyone can follow it to some extent. Also, this post is Literate Haskell; you can copy and paste the entire text of the post directly into a file with an .lhs extension and load it in a Haskell interpreter like ghci.

WARNING: spoilers ahead!! If you’d like to work out the challenges yourself, stop reading now.

Okay, first things first. Some smart person has already put together a nice library for efficiently listing prime numbers, so let’s use it:

> import Data.Numbers.Primes (primes, isPrime)

A few other odds and ends we will need:

> import Control.Arrow ((>>>), (&&&))
> import Data.List (tails, find, intercalate)
>
> infixl 0 >$>
> x >$> f = f x   -- backwards function application

First, let’s write a function, primeSums, which takes a target number and figures out if it can be expressed as a sum of consecutive primes.

> -- Find runs of consecutive primes that sum to n.
> primeSums :: Int -> [(Int, Int)]

This says that primeSums takes an Int as input and returns a list of ordered pairs. The first number in each pair is the index of the first prime number in the run (2 is the 0th prime, 3 is the 1st prime, and so on), and the second number is the length of the run. For example, primeSums 120 = [(8,4)], since 120 is equal to the sum of four consecutive primes starting with the ninth prime—that is, 120 = 23 + 29 + 31 + 37.

We’re going to write primeSums by starting with the list of all prime numbers and transforming it incrementally by applying a series of functions to it. Note that it really is the infinite list of all prime numbers—one of the very nice features of Haskell is that it is lazy, which means that we can program with lists and other data structures that are infinite in theory, but only the parts we actually end up using will ever be computed.

Since we are looking for runs of at least two consecutive primes that add up to n, primes larger than half of n are obviously of no use. So we start by keeping only the primes up to half of n [Thanks to Steve for the correction; see comment below!]:

> primeSums n =
>   primes >$> takeWhile (< n)

Next we apply tails which gives us the list of all suffixes, or "tails", of a list. For example, tails [1,2,3] = [[1,2,3], [2,3], [3], []]. This is so we can consider all the possible starting places for runs of consecutive primes. We don’t want the empty list at the end so we use init which throws away the last element of a list.

>          >>> tails
>          >>> init

We now have a list of lists of prime numbers, each list starting at a different place. We want to apply the same processing to each list, which we can do with map:

>          >>> map (

With each list of primes we first compute all the partial sums using scanl (+) 0. For example, scanl (+) 0 [1,2,3,4] = [0,1,3,6,10].

>                scanl (+) 0

The 0 at the beginning is clearly useless so we throw it away with tail.

>            >>> tail

We only care about sums which are less than or equal to our target number n, so we keep only that part of the list. (Note that due to Haskell’s laziness, the bigger sums that would have come later in the list won’t even get computed.)

>            >>> takeWhile (<= n)

So at this point, if we had started with the list [2,3,5,7,11,13,17,19,23] and the target number 28, we would have [2,5,10,17,28] (2 = 2; 5 = 2 + 3; 10 = 2 + 3 + 5; …; 28 = 2 + 3 + 5 + 7 + 11; and the sums including 13 and beyond are too big, so they get thrown away).

Now, all we really care about is the last number in the list we have so far, since it is either equal to the target number, meaning we have found a sequence of consecutive primes summing to it, or it is less than the target number. We also want to remember the length of the list since that tells us how long the run of consecutive primes is.

>            >>> length &&& last
>              )

&&& allows us to apply two functions to a single value and get out a pair of values. For example, in this case (length &&& last) [2,5,10,17,28] = (5,28).

Okay, we are now done processing each list of primes, and in the place of each list we are left with a pair of numbers; those pairs whose second number is equal to the target number mark places where a consecutive run of primes added up to the target number. We want to remember which pair is which (so we know where each run started), so we first use zip [0..] to combine each pair with its index.

>          >>> zip [0..]

Now the elements of the list look like, for example, (0, (5,28)) which means that starting from index 0 there is a sequence of five consecutive primes which add up to 28.

We only want to keep runs that actually add up to the target number; filter keeps only those elements of a list that satisfy some condition:

>          >>> filter ((== n) . snd . snd)

Finally, we can throw away the sum from each element since we know what it is; we keep only the starting index and length of each run.

>          >>> map (fst &&& fst . snd)

And that’s it! We can try it out by loading this file in ghci and evaluating primeSums applied to various arguments at the prompt:

*Main> primeSums 28
[(0,5)]
*Main> primeSums 311
[(4,11),(10,7),(15,5),(25,3)]
*Main> primeSums 2011
[(36,11),(120,3)]

As expected, 28 is the sum of five consecutive primes starting with the 0th (that is, 2). 311 can be written as a sum of consecutive primes in four different ways, starting with a run of 11 consecutive primes starting with the prime at index 4 (that is, 11). And sure enough, 2011 can be written as the sum of 11 consecutive primes starting at index 36, or three consecutive primes starting at index 120. Uh… which primes are those? Let’s write a little function to "expand" runs so we can actually see what they are:

> expandRun :: (Int, Int) -> [Int]
> expandRun (i,n) = primes >$> drop i >>> take n

Given the pair (i,n) representing a starting index i and number of primes n, we simply take the list of primes, drop the first i of them, and then take the next n.

*Main> map expandRun (primeSums 2011)
[[157,163,167,173,179,181,191,193,197,199,211],[661,673,677]]

That’s better! Apparently 2011 = 157 + 163 + 167 + 173 + 179 + 181 + 191 + 193 + 197 + 199 + 211 = 661 + 673 + 677. And why not add them up just to make sure?

*Main> map (sum . expandRun) (primeSums 2011)
[2011,2011]

Oh good!

Now on to the next bit of Patrick’s challenge. When is the next time the year will be expressible as the sum of a prime number of consecutive prime numbers? First, we need a function to test whether this is true:

> isSumOfPrimeConsecPrimes = any (isPrime . snd) . primeSums

isSumOfPrimeConsecPrimes simply computes the primeSums of a number and then checks whether any of the second (snd) elements of the pairs returned are prime (remember, the second elements give the lengths of the runs). Let’s make sure it’s true for 2011:

*Main> isSumOfPrimeConsecPrimes 2011
True

Sure enough. And when is the next one?

*Main> find isSumOfPrimeConsecPrimes [2012..]
Just 2015
*Main> primeSums 2015
[(76,5)]

Looks like we won’t have to wait very long; 2015 is a sum of five consecutive primes. But wait a minute, 2011 is a sum of a prime number of consecutive primes in two different ways. When is the next time that will happen? (Patrick didn’t ask that, but I’m just curious!)

iSoPCPiMW (is Sum of Prime Consecutive Primes in Multiple Ways) checks just what you would think it does; we get to specify how many ways k we require. It keeps only runs whose length is prime, counts how many such runs there are, and tells us whether that number is no less than k.

> iSoPCPiMW k = (>=k) . length . filter (isPrime . snd) . primeSums
*Main> find (iSoPCPiMW 2) [2012..]
Just 2155
*Main> find (iSoPCPiMW 3) [2012..]
Just 2303
*Main> find (iSoPCPiMW 4) [2012..]
Just 2867

Well, we have to wait a bit longer for this. The year won’t be expressible as a sum of a prime number of consecutive primes in two ways again until 2155! And as you can see, just for fun I also found the next times when there will be three or four such ways. (Interestingly, I also tried searching for a year with five ways; but my computer is still searching after ten minutes, even though all of the above searches took less than a second! So either there isn’t one, or it won’t happen for a very long time, but I don’t know which. Maybe you can figure it out?)

What about Patrick’s last question? 2015 clearly isn’t prime. Well, we can just search through the list of primes instead of searching the list of all numbers bigger than 2011:

*Main> find isSumOfPrimeConsecPrimes $ dropWhile (<=2011) primes
Just 2081

2081, eh? That’s still a long time from now, but unlike 2155 there’s a small chance I’ll still be alive then (I would be 99)!

Well, I hope this has been fun. Feel free to leave a comment with any questions. If you’re interested in learning more about Haskell, you can download it here, and I recommend starting by reading Learn You a Haskell for Great Good!.

About these ads
This entry was posted in arithmetic, number theory, primes, programming and tagged , , , . Bookmark the permalink.

8 Responses to Prime Time in Haskell

  1. Nicely done! I might have to learn me some Haskell!

  2. Eric says:

    Might be interested in a Mathematica solution (hope I can format in a not too distracting way):

    target=2011;
    candidates=Tuples[Range[PrimePi[target/2]+1],2];
    candidates=Select[candidates,#[[1]]<#[[2]]&];

    PrimeListSum[n_Integer?Positive,m_Integer?Positive]:=Plus@@Prime/@Range[n,m];

    Select[candidates,PrimeListSum[#[[1]],#[[2]]]==target&]//Timing

    Out[5]= {0.700243,{{37,47},{121,123}}}

  3. Steve says:

    Just to be a jerk, I thought I’d mention:

    > Since we are looking for runs of at least two consecutive primes
    > that add up to n, primes larger than half of n are obviously of no use.

    Not necessarily true, depending on n — for instance, if I’m looking for two consecutive primes that add up to 12, then 5 and 7 would work, even though 7 is greater than half of 12.

    Anyway, as I said: I’m a jerk. =)

  4. Basil says:

    WoW. Nice work. I love computing power.

  5. Brent says:

    Steve: you jerk! ;) You’re absolutely right though. I should know better than to use the word “obviously”, since it obviously means I can’t be bothered to explain it and am therefore wrong, more likely than not. Anyway, I’ve corrected it above; thankfully it didn’t invalidate any of the results, although the primeSums function was definitely incorrect in general.

  6. JM says:

    I feel like a holiday should be created (following in the footsteps of pi day) to celebrate the greatness of prime numbers. I’m all for the 23rd day of the 5th month.

  7. Pingback: Walking Randomly » Carnival of Mathematics #73 – Chuck Norris Edition

  8. Pingback: Learn You a Haskell! | The Math Less Traveled

Comments are closed.