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

Problem 49: What 12-digit number do you form by concatenating the three terms in this sequence?

This was a tough problem. Even after dreaming up an algorithm which would perform well enough, it went through several tweaks until I was satisfied. We resort to an array (though no mutation) and sequence expression in order to achieve short-circuiting and overcome F#’s lack of break and continue (again). This algorithm runs in about 1 minute and 500 milliseconds, just over the Project Euler suggested limit. I was able to uglify it a bit to get it to run in 59 seconds and 700 milliseconds, but I won’t even bother showing that here.

let problem49b =
    let fourDigitPrimes =
        {1000..9999}
        |> Seq.filter isPrime
        |> Seq.toArray

    let arithmeticTriples = seq {
        for i in {0..fourDigitPrimes.Length-1} do
            for j in {i+1..fourDigitPrimes.Length-1} do
                for k in {j+1..fourDigitPrimes.Length-1} do
                    let a = fourDigitPrimes.[i]
                    let b = fourDigitPrimes.[j]
                    let c = fourDigitPrimes.[k]
                    //is arithmetic seq other than given in question
                    if b-a = c-b && a <> 1487 then 
                        yield (a,b,c)
    }
    
    let arePerms (a,b,c) = 
        let setA = Set(Digits.fromInt a)
        let setB = Set(Digits.fromInt b)
        setA = setB && setB = Set(Digits.fromInt c)
        
    arithmeticTriples
    |> Seq.find arePerms
    |> (fun (a,b,c) ->
            [Digits.fromInt a; Digits.fromInt b; Digits.fromInt c] 
            |> Seq.concat 
            |> Digits.toInt64)

Problem 48: Find the last ten digits of the series, 1^1 + 2^2 + 3^3 + ... + 1000^1000.

Again, easy problem using functional composition, the “hard” work was done in our previous implementation of the Digits module. The answer was produced surprisingly fast. I would like a better solution to finding the last 10 digits than to convert the entire sequence to an array.

let problem48a =
    let digits =
        {1..1000}
        |> Seq.sumBy (fun n -> BigInteger.Pow(BigInteger(n), n))
        |> Digits.fromBigInt
        |> Seq.toArray

    Array.sub digits (digits.Length - 10) 10
    |> Digits.toInt64 //just makes it easier to copy and paste result

Problem 47: Find the first four consecutive integers to have four distinct primes factors. What is the first of these numbers?

Not the best performing solution (perhaps inherently so, though), but a shining example of the power of functional composition (the hard work was done in implementing cfactorize in a previous problem, everything else is just composing built-in functional operators).

let problem47a =
    (Seq.initInfinite (fun i -> let n = i+3 in (n, cfactorize n)) //starting with 0+3=2
    |> Seq.windowed 4
    |> Seq.find (fun arr -> arr |> Seq.forall (fun (_,lst) -> lst |> List.length = 4))).[0]
    |> fst

Problem 46: What is the smallest odd composite that cannot be written as the sum of a prime and twice a square?

The difficult part of this problem was working around F#'s lack of break, continue, and return in for loops. hasForm would be most easily expressed using nested for loops but would require break, continue, and return to perform acceptably. Another desire is support for infinite range expressions, so that oddComposites could be written as {9..2..} |> Seq.filter (isPrime>>not).

let problem46b = 
    let form p s = p + 2*(s*s)

    //making sure to exclude 1    
    let oddComposites = 
        Seq.initInfinite id 
        |> Seq.filter (fun i -> i <> 1 && i%2=1 && isPrime i |> not)
    
    //cached for performance
    let primes = 
        Seq.initInfinite id 
        |> Seq.filter isPrime 
        |> Seq.cache

    //exclusive
    let primesUpto n = 
        primes 
        |> Seq.takeWhile ((>=) n)

    //think nested loops
    let hasForm n = 
        primesUpto n
        |> Seq.exists 
            (fun p ->
                ({1..n-1} 
                |> Seq.map (fun s -> form p s)
                |> Seq.find (fun r -> (r = n) || (r > n))) = n)

    oddComposites
    |> Seq.find (hasForm>>not)

Problem 45: After 40755, what is the next triangle number that is also pentagonal and hexagonal?

Here we reused the technique from Problem 44 for quickly asserting that a number is in the range of a function by verifying that n = f (f-1(n)).  We could also have tested that f-1(n) is integral.

let problem45a =
    let pent n = n*(3L*n-1L)/2L
    let pentInv n = (1L + sqrtL(24L*n + 1L))/6L
    let isPent n = (n = (pent (pentInv n)))
    
    let hex n = n*(2L*n-1L)
    let hexInv n = (1L + sqrtL(8L*n + 1L))/4L
    let isHex n = (n = (hex (hexInv n)))
    
    let tri n = n*(n+1L)/2L
    let rec triSeq n = seq {yield tri n; yield! triSeq (n+1L)}
    
    triSeq 286L
    |> Seq.find (fun n -> isPent n && isHex n)

Problem 44: Find the pair of pentagonal numbers, Pj and Pk, for which their sum and difference is pentagonal and D = |Pk - Pj| is minimised; what is the value of D?

This problem was a big disappointment.  It turns out that the first pair you find for which their sum and difference is pentagonal is the solution.  The interesting yet difficult minimization of D is a complete waste of time.  More than that, the running time of the solution I came up with is beyond reach and to keep from overflows big numbers would be required.  Looking at the Project Euler forum, most users totally glossed over the requirement for the minimization of D, though a couple others did protest.  What follows is the theory and implementation for the true solution of this problem (which as I’ve pointed out, actually won’t ever produce the result because 1) integer overflow, and 2) unreachable running time even when we do use big numbers).

The derivative of p is linear (p' = (3n^s - n)/2), so p is strickly increasing. Therefore, if p(i-1) - p(i) > min {p(s) - p(t) where i > s > t and p(s) and p(t) are pentagonal numbers for which their sum and difference is pentaganol} then the right-hand side is the minimization D.

let problem44a =
    let p n = n*(3*n-1)/2 //p = (3n^2-n)/2
    let pinv n = (1 + sqrtn(24*n + 1))/6 //inverse for positive range of p
    let isPentagonal t = t = (p (pinv t))
    
    let testPair p1 p2 = 
        (isPentagonal (p1-p2)) && (isPentagonal (p1+p2))
            
    let rec loop n minDiff =
        let pn = p n
        if ((p (n+1)) - pn) > minDiff then minDiff
        else
            let rec find m = 
                let pm = p m
                if m = 0 then None
                elif testPair pn pm then Some(pm)
                else find (m-1)
            match find (n-1) with
            | Some(pm) -> loop (n+1) (min (pn-pm) minDiff) //take the first pair
            | None -> loop (n+1) minDiff
    loop 2 System.Int32.MaxValue

Problem 43: Find the sum of all 0 to 9 pandigital numbers with the given property.

Combining our Digits and permutations implementations with F#’s built-in Seq functions, we really see how functional programming shines in describing complex algorithms simply.

let problem43a =  
    let hasProperty p = 
        let ddd = p |> Seq.skip 1
                    |> Seq.windowed 3
                    |> Seq.map Digits.toInt 
        
        Seq.zip ddd [2;3;5;7;11;13;17]
        |> Seq.forall (fun (a,b) -> a % b = 0)
            
    permutationsAsc [0;1;2;3;4;5;6;7;8;9]
    |> Seq.filter hasProperty
    |> Seq.sumBy Digits.toInt64

Problem 42: How many triangle words does the list of common English words contain?

This problem is designed to make use of several tools.

let problem42a =
    let triangle n = let n' = n |> float in (0.5)*(n')*(n'+1.) |> int
    let tseq = Seq.initInfinite (fun i -> triangle (i+1)) |> Seq.cache
    
    let wordValue (s:string) = s |> Seq.sumBy (fun c -> (c |> int)-64)
    
    let text = System.IO.File.ReadAllText("words.txt")
    let words = text.Split([|',';'"'|], System.StringSplitOptions.RemoveEmptyEntries) |> Seq.readonly
    
    let isTriangleWord w =
        let wv = wordValue w
        tseq
        |> Seq.takeWhile (fun t -> t <= wv)
        |> Seq.exists ((=) wv)
        
    words
    |> Seq.filter isTriangleWord
    |> Seq.length

Problem 41: What is the largest n-digit pandigital prime that exists?

My first approach was to try and find the first odd number descending from 987654321 which was both prime and pandigital, but despite attempts at performance tweaks this would never finish executing.

let problem41a = 
    let isPandigital n =
        let dlist = n |> digits |> List.ofSeq
        (List.sort dlist) = [1..List.length dlist]
        
    let rec loop n =
        if isPrime n && isPandigital n then n
        else loop (n-2)
    loop 987654321

Then I realized I would have much better luck generating pandigital numbers and testing them for primality using the lexical permutation algorithm from Problem 24.  The following is a generic implementation with a immutable sequence wrapper.

First we need a couple operators: comparer converts an F# function into an IComparer so we can use the System.Array.Sort overload which allows in-place sub range sorting.  flip reverses the parameters of two parameter functions.

let comparer f = { new System.Collections.Generic.IComparer<'a> with member self.Compare(x,y) = f x y }

let flip f x y = f y x

Next is our permutations function which is generic and takes a comparison function.  The inner permute function mutates the given perm array returning false when the array is at its last permutation.  Then a sequence expression is used to generate a sequence with immutable copies of the array after each permutation.  Finally, permutationsAsc and permutationsDesc are convenience partial applications of permutations.  Despite the fact that the array is copied for each yield of the sequence, this performs quite well.

let permutations f e =
    ///Advances (mutating) perm to the next lexical permutation.
    let permute (perm:'a[]) (f: 'a->'a->int) (comparer:System.Collections.Generic.IComparer<'a>) : bool =
        try
            //Find the longest "tail" that is ordered in decreasing order ((s+1)..perm.Length-1).
            //will throw an index out of bounds exception if perm is the last permuation,
            //but will not corrupt perm.
            let rec find i =
                if (f perm.[i] perm.[i-1]) >= 0 then i-1
                else find (i-1)
            let s = find (perm.Length-1)
            let s' = perm.[s]
            
            //Change the number just before the tail (s') to the smallest number bigger than it in the tail (perm.[t]).
            let rec find i imin =
                if i = perm.Length then imin
                elif (f perm.[i] s') > 0 && (f perm.[i] perm.[imin]) < 0 then find (i+1) i
                else find (i+1) imin
            let t = find (s+1) (s+1)
                
            perm.[s] <- perm.[t]
            perm.[t] <- s'

            //Sort the tail in increasing order.
            System.Array.Sort(perm, s+1, perm.Length - s - 1, comparer)
            true
        with
        | _ -> false
       
    //permuation sequence expression 
    let c = f |> comparer
    let freeze arr = arr |> Array.copy |> Seq.readonly
    seq { let e' = Seq.toArray e
          yield freeze e'
          while permute e' f c do
              yield freeze e' }
              
let permutationsAsc e = permutations compare e
let permutationsDesc e = permutations (flip compare) e

Now we have the solution to our problem.

let problem41c =    
    let rec loop n =
        let maybe = 
            {n..(-1)..1}
            |> permutationsDesc
            |> Seq.map Digits.toInt
            |> Seq.tryFind isPrime
            
        match maybe with
        | Some(value) -> value
        | None -> loop (n-1)
        
    loop 9

Note: apparently if the sum of the digits of a number is divisible by 3 then the number cannot be prime.  This isn’t obvious to me, but it’s late so I’ll look into to it more tomorrow.  However, if we take this fact we find that the number of digits must then either be 4 or 7.  Applying the 7 digit upper bound, both solution (a) and (b) run instantaneously.

Problem 40: If dn represents the nth digit of the fractional part, find the value of the following expression. d1 * d10 * d100 * d1000 * d10000 * d100000 * d1000000

An imperative solution seems appropriate.

let problem40a =
    let sb = System.Text.StringBuilder()
    let mutable n = 1
    while sb.Length <= 1000000 do
        sb.Append(n) |> ignore
        n <- n+1
        
    let mutable prod = 1
    for i in 0..6 do
        prod <- prod * ((sb.[(pown 10 i)-1] |> int) - 48) 
        
    prod

Problem 39: For which value of p ≤ 1000, is the number of solutions maximised?

Brute force algorithm with no real optimizations.

let problem39c =
    let countTriplets p = 
        let mutable count = 0
        for a in 1..p do
            for b in a..p do
                let c = p - (a + b) |> float
                let c' = a*a + b*b |> float |> sqrtF
                if c = c' then
                    count <- count + 1
        count
        
    {1..1000}
    |> PSeq.maxBy countTriplets

Problem 38: What is the largest 1 to 9 pandigital that can be formed by multiplying a fixed number by 1, 2, 3, ... ?

Wow, interesting - tough!  And I finally found Seq.choose (I was looking for Seq.pickAll).

let problem38a =
    let isPosPd (nstr:string) = //we know our chars are digits
        let len = nstr |> Seq.length
        len = 9 
        && (nstr.Contains("0") |> not) 
        && len = (nstr |> Seq.distinct |> Seq.length)
    
    let tryFind n =
        let rec loop i =
            let prodConcats = (seq {for j in 1..i -> n*j |> string }) |> Seq.fold (+) ""
            if isPosPd prodConcats then Some(prodConcats |> int)
            elif i < 9 then loop (i+1)
            else None
        loop 2 //start with 2 since n > 1
        
    {1..9999} //since n > 1, 9999 is clearly an upper limit
    |> PSeq.choose tryFind
    |> PSeq.max

Problem 37: Find the sum of all eleven primes that are both truncatable from left to right and right to left.

Kind of a variation on Problem 35.  I was surprised by how fast it ran!

let problem37a =
    let isTruncatablePrime n =
        if n |> isPrime |> not then false
        else
            let digs = n |> Digits.fromInt
            let truncations =
                seq { for i in 1..(Seq.length digs)-1 -> 
                          digs |> Seq.take i |> Digits.toInt
                      for i in 1..(Seq.length digs)-1 -> 
                          digs |> Seq.skip i |> Digits.toInt }
            Seq.forall isPrime truncations
            
    let rec odds n = seq {yield n; yield! odds (n+2)}
    odds 11 //skip the single digit primes
    |> Seq.filter isTruncatablePrime
    |> Seq.take 11 //we are given there are only 11
    |> Seq.sum

Problem 36: Find the sum of all numbers less than one million, which are palindromic in base 10 and base 2.

Reusing isPalindrome from Problem 4 and leveraging System.Convert for obtaining the binary representation of a given decimal number, this problem is trivial.

let problem36a =
    {1..999999}
    |> Seq.filter (fun n -> (isPalindrome (n|>string)) && (isPalindrome (System.Convert.ToString(n, 2))))
    |> Seq.sum

Problem 35: How many circular primes are there below one million?

Now we get to use our Digits module.  The most interesting part of this algorithm is performing digit rotations (the assignment of r), which went through several revisions before I got it just right.

let problem35c = 
    let isCircularPrime n =
        if not (isPrime n) then false
        else 
            let digs = n |> Digits.fromInt
            let rec loop i = 
                if i = (Seq.length digs)-1 then true
                else
                    let i = i + 1
                    let r = //rotate by i
                        Seq.append 
                            (digs |> Seq.skip i) 
                            (digs |> Seq.take i)
                        |> Digits.toInt
                    if isPrime r then loop i
                    else
                        false
            loop 0

    //we are given that there are 13 such primes below 100
    ({101..2..999999}
    |> Seq.filter isCircularPrime
    |> Seq.length) + 13

Digits

We’ve had a lot of problems involving digit manipulation and so far have implemented ad-hoc string parsing solutions.  But this approach is generally pretty slow.  For primitive integer types (int32 and int64), we can gain two to three times speed increase using arithmetic methods.  On the other hand, generating a sequence of digits from a bigint is actually faster using the former method.  The following module provides type optimized conversions between digit sequences and numbers.

    module Digits =
        //at least twice as fast as fastParse
        let fromInt64 (n:int64) =
            if n = 0L then Seq.singleton 0
            else
                let mutable n = if n < 0L then (abs n) else n
                let mutable powten = 0
                
                while (n/(pown 10L powten) >= 10L) do
                    powten <- powten+1
                    
                let darr = Array.create (powten+1) 0
                let mutable i = 0
                while (powten >= 0) do
                    let d = n/(pown 10L powten)
                    darr.[i] <- d |> int
                    n <- n - d*(pown 10L powten)
                    i <- i + 1
                    powten <- powten-1
                    
                Seq.readonly darr
                            
        let fromInt (n:int) = fromInt64 (n |> int64)
        
        ///parse a positive or negative number string.  if other characters (besides leading '-'
        ///and integers) are present results will be unpredictable.
        let uncheckedParse (nstr:string) =
            if nstr.[0] = '-' then //fast with least num checks
                nstr
                |> Seq.skip 1
                |> Seq.map (fun c -> (c |> int) - 48)
                |> Seq.cache
            else
                nstr
                |> Seq.map (fun c -> (c |> int) - 48)
                |> Seq.cache
        
        ///filter out non-digits
        let parse nstr = 
            uncheckedParse nstr
            |> Seq.filter (fun n -> n >= 0 && n <= 10)
        
        //not sure if using fromString here is faster than direct map string |> parse
        let fromBigInt (n:bigint) = n |> string |> uncheckedParse
                
        let toInt64 (digs:seq<int>) = 
            digs
            |> Seq.fold 
                (fun (e, sum) d -> let e = e-1 in (e, Checked.(+) sum ((d |> int64)*(pown 10L e))))
                (Seq.length digs, 0L)
            |> snd
                    
        ///note: (System.Int32.MinValue) |> Digits.fromInt |> Digits.toInt results
        ///in overflow exceptions since abs(int.MinValue) > int.MaxValue
        let toInt (digs:seq<int>) = 
            digs
            |> Seq.fold 
                (fun (e, sum) d -> let e = e-1 in (e, Checked.(+) sum (d*(pown 10 e))))
                (Seq.length digs, 0)
            |> snd
        
        open System.Text
        let toBigInt (digs:seq<int>) = 
            bigint.Parse(
                digs 
                |> Seq.fold (fun (sb:StringBuilder) d -> ignore <| sb.Append(d) ; sb) (StringBuilder())
                |> string
            ) //not sure if actually fast for bigint

Problem 34: Find the sum of all numbers which are equal to the sum of the factorial of their digits.

The upper bound is 2540160, since (factorial 9)*7 = 2540160 and (factorial 9)^x grows faster than 10^x for x >= 7.

let problem34b =   
    {3..2540160}
    |> Seq.filter (fun i -> i = (i |> string |> Seq.map (string >> int) |> Seq.sumBy factorial))
    |> Seq.sum

Problem 33: Discover all the fractions with an unorthodox cancelling method.

With F#, there are two ways to solve this: the easy way, and the slightly easier way.  The first way, we’ll work entirely with integers, including an implementation of gcd using the Euclid’s algorithm for reducing our final fraction.

gcd and accompanying abs are a sure candidates for inclusion in our Numerics module.

let inline abs (g:G<'a>) n = 
    if n < g.zero then n*g.negone else n
    
let inline gcd (g:G<'a>) n m =
    if n = g.zero || m = g.zero then 
        raise (System.ArgumentOutOfRangeException("n and m must non-zero"))
    else
        let rec gcd n m = 
            if n = m then n
            else
                if n > m then gcd (n-m) m
                else gcd n (m-n)
        gcd (abs g n) (abs g m)

And then our solution is fairly straightforward:

let problem33a =
    let unorthodoxPairs =
        [for numer = 10 to 99 do
            for denom = (numer+1) to 99 do
                let reduced = (numer |> float)/(denom |> float)
                let digitPair n = //decompose digit list of n into a pair
                    let nDigits = n |> string |> Seq.map (string >> float) |> Seq.toArray
                    (nDigits.[0], nDigits.[1])
                let (a,b) = digitPair numer
                let (c,d) = digitPair denom
                let isUnorthodox w x y z = (w/x) = 1. && (y/z) = reduced
                if isUnorthodox a c b d || 
                   isUnorthodox a d b c || 
                   isUnorthodox b d a c || 
                   isUnorthodox b c a d then yield (numer, denom)]
    
    let product = unorthodoxPairs |> List.fold (fun (w,x) (y,z) -> (w*y,x*z)) (1,1)
    snd product / (gcd (fst product) (snd product)) 

Next is our very easy solution which takes advantage of BigRational included in the F# PowerPack:

let problem33b =
    let unorthodoxPairs =
        [for numer = 10 to 99 do
            for denom = (numer+1) to 99 do
                let reduced = (numer |> float)/(denom |> float)
                let digitPair n = //decompose digit list of n into a pair
                    let nDigits = n |> string |> Seq.map (string >> float) |> Seq.toArray
                    (nDigits.[0], nDigits.[1])
                let (a,b) = digitPair numer
                let (c,d) = digitPair denom
                let isUnorthodox w x y z = (w/x) = 1. && (y/z) = reduced
                if isUnorthodox a c b d || 
                   isUnorthodox a d b c || 
                   isUnorthodox b d a c || 
                   isUnorthodox b c a d then yield BigRational.Parse(sprintf "%d/%d" numer denom)]
    
    let product = unorthodoxPairs |> List.fold (*) 1N
    product.Denominator

Here’s a third version which avoids string parsing and the need for conversions between int and float:

let problem33c =
    let unorthodoxPairs =
        [for numer = 10. to 99. do
            for denom = (numer+1.) to 99. do                
                let digitPair n =
                    let tensPlace = System.Math.Truncate(n/10.)
                    (tensPlace, n-tensPlace*10.)    

                let (a,b) = digitPair numer
                let (c,d) = digitPair denom

                let isUnorthodox w x y z = (w/x) = 1. && (y/z) = numer/denom
                if isUnorthodox a c b d || 
                   isUnorthodox a d b c || 
                   isUnorthodox b d a c || 
                   isUnorthodox b c a d then yield BigRational.Parse(sprintf "%.0f/%.0f" numer denom)]
    
    let product = unorthodoxPairs |> List.fold (*) 1N
    product.Denominator

Problem 32: Find the sum of all products whose multiplicand/multiplier/product identity can be written as a 1 through 9 pandigital.

While having lots of pieces, this problem is mostly straightforward.  The only part that required any analysis is determining sufficient lower and upper bounds for the product, since the worst case range (1 to 999999999) is far too large.  Letting p be the product, mn be the multiplican, and mr be the multiplier, we observe that if mn * mr = p, then (len p) = ((len mn) + (len mr)) + (-1 or 0).  Which implies the length of the product is cannot be less than 4 or greater than 5 or we will have either too few or too many digits in our pandigital set. Thus the range  may be restricted to 1234 to 98765.

let problem32c =
    let hasDistinctPositiveDigits (nstr:string) = 
        not (nstr.Contains("0")) && nstr.Length = (nstr |> Seq.distinct |> Seq.length)
    
    let isPositivePanDigitalSet mnstr mrstr pstr =
        let concatedDigitStr:string = mnstr + mrstr + pstr
        concatedDigitStr.Length = 9 && hasDistinctPositiveDigits concatedDigitStr
        
    let findPositivePanDigitalSet p =
        let pstr = p |> string
        if not (hasDistinctPositiveDigits pstr) then None
        else
            let sr = sqrtn p
            let rec loop mn =
                if mn > sr then None
                elif p % mn = 0 then
                    let mr = p/mn
                    if isPositivePanDigitalSet (mn |> string) (mr |> string) pstr then 
                        Some(mn, mr, p) 
                    else loop (mn + 1)
                else loop (mn + 1)
            loop 2
    
    {1234..98765} //product must be between 4 and 5 digits
    |> Seq.map findPositivePanDigitalSet
    |> Seq.filter Option.isSome
    |> Seq.sumBy (fun (Some(_,_,p)) -> p)

Problem 31: How many different ways can £2 be made using any number of coins?

More fun with recursive sequence expressions! In this solution, we generate all of the combinations instead of only counting them.

let problem31c =
    let combinations amt coins = //produce all ascending combinations starting with start
        let rec combinations combo =
            seq{ let sum = List.sum combo
                 if sum = amt then 
                     yield combo
                 elif sum < amt then 
                     yield! coins 
                            |> Seq.filter ((<=) combo.Head) //to generate combinations instead permutations
                            |> Seq.collect (fun c -> combinations (c::combo)) }
        seq {for start in coins do yield! combinations [start]}  //concat all ascending combinations for each start
    combinations 200 [1;2;5;10;20;50;100;200] |> Seq.length

Problem 30: Find the sum of all the numbers that can be written as the sum of fifth powers of their digits.

We must find an upper limit for numbers to search.  As the length of a number increases, 9**5 is the greatest additional contribution a digit can make to the sum of 5th power digits. Meanwhile, the number itself is increasing by powers of 10, which will eventually overcome the rate of increase given summing 5th powers of digits.  So, when we find that the sum of 5th power digits of 9999999 (of length 7) equals 413343 (of length 6), we can be sure that no higher number can produce 5th power digit sum of sufficient length.  (a non-rigorous proof, I know).

let problem30a =
    let powdigsum n e = n |> string |> Seq.map(string >> int32) |> Seq.sumBy(fun i -> pown i e)
    {2..9999999}
    |> Seq.filter(fun n -> (powdigsum n 5) = n)
    |> Seq.sum

Problem 29: How many distinct terms are in the sequence generated by a^(b) for 2 ≤ a ≤ 100 and 2 ≤ b ≤ 100?

Too easy.

let problem29a =
    seq {for a in 2I..100I do for b in 2..100 do yield bigint.Pow(a,b)}
    |> Seq.distinct 
    |> Seq.length

Problem 28: What is the sum of the numbers on the diagonals in a 1001 by 1001 number spiral?

Another pleasant one where by writing down the first few terms of the sequence (by dimension) of sequences (corners in a dimension), it was easy to spot a recursive pattern based on the last element of the previous sequence.  The neatest part is how apt recursive sequence expressions are for capturing this algorithm.  For interest, I show my first-pass solution followed by a refactored version.

let problem28a =
    let diags upperdim =
        let rec diags dim prev =
            seq { let prev = (dim-1) + prev
                  yield prev
                  let prev = (dim-1) + prev
                  yield prev
                  let prev = (dim-1) + prev
                  yield prev
                  let prev = (dim-1) + prev
                  yield prev
                  if dim <> upperdim then 
                    yield! diags (dim+2) prev }
        diags 3 1
    (diags 1001 |> Seq.sum) + 1 //add center 1
                
let problem28b =
    let diags upperdim =
        let rec diags dim prev =
            seq { let next i = i*(dim-1) + prev
                  yield! {1..4} |> Seq.map next
                  if dim <> upperdim then 
                    yield! diags (dim+2) (next 4) }
        diags 3 1
    (diags 1001 |> Seq.sum) + 1 //add center 1

The next two versions use an infinite sequence which includes the center 1.  The later of the two takes advantage of the fact that it’s easy to count how many terms we need to sum.

let problem28c =
    let diags =
        let rec diags dim prev =
            seq { let next i = (dim, i*(dim-1) + prev)
                  yield! {1..4} |> Seq.map next
                  yield! diags (dim+2) (next 4 |> snd) }
        seq { yield (1,1) ; yield! diags 3 1 }
    diags |> Seq.takeWhile(fun (dim,_)  -> dim <= 1001) |> Seq.sumBy snd
    
let problem28d =
    let diags =
        let rec diags dim prev =
            seq { let next i = i*(dim-1) + prev
                  yield! {1..4} |> Seq.map next
                  yield! diags (dim+2) (next 4) }
        seq { yield 1 ; yield! diags 3 1 }
    diags |> Seq.take((4*500)+1) |> Seq.sum

Problem 27: Find the product of the coefficients, a and b, for the quadratic expression that produces the maximum number of primes for consecutive values of n, starting with n = 0.

This one was straight forward, just had to follow instructions (I didn’t apply any mathematical optimizations, because it looked like if once I got started, I’d be able to solve it all by hand).  I needed to adjust isPrime to return false for values of n less than 2.  I got carried away implementing it SQL style; one long run-on sentence (excluding problem27a, not a single let binding!).  Parallelization improves speed by about 25 percent.

let problem27a =
    seq { for a in -999..999 do for b in -999..999 do yield (a,b) }
    |> PSeq.map(fun (a,b) -> (a*b, Seq.initInfinite id
                                   |> Seq.takeWhile (fun n -> n*n + a*n + b |> isPrime)
                                   |> Seq.length))
    |> PSeq.maxBy snd
    |> fst

Problem 26: Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its decimal fraction part.

After taking a break for a few days, having spent 3 consecutive weeks solving problems 1-25 every free hour and deep into the night, this turned out to be a pleasant return despite it’s intimidating appearance.  I simply worked out the first 7 values of d on paper using long division, spotted a pattern corresponding to recurring cycle lengths, and churned out an implementation of the algorithm in F#.  I am finding that looping via recursion is a very gentle way to implement algorithms: you simply reason out each case one at a time and return when you come to the end of a branch.  No worrying what happens next, just pass it on.

let problem26c =     
    let cycleLength d =
        let rec cycleLength (steps:ResizeArray<int>) step =
            if d > step then 
                cycleLength steps (step*10)
            else 
                if steps.Contains(step) then
                    (d, steps.Count - steps.IndexOf(step))
                else
                    steps.Add(step)
                    let step = step - d*(step/d)
                    if step = 0 then
                        (d, 0)
                    else
                        cycleLength steps step
        cycleLength (ResizeArray<int>()) 1
    //on dual core, runs in 60ms vs. 90ms when parallelized, 
    //and twice as fast with larger ranges
    {2..999}
    |> PSeq.map cycleLength 
    |> PSeq.maxBy snd 
    |> fst

Parallelization with PLINQ

The past couple of days I've been interested in exploring PLINQ with F#.  But Since I’m using F# 2.0 for Visual Studio 2008 (i.e. .NET 2.0 build), I had to 1) Install the 3.5 SP1 Reactive Extensions (Rx) for .NET which include a back-port of PLINQ (http://msdn.microsoft.com/en-us/devlabs/ee794896.aspx), and 2) Download the F# PowerPack 2.0 May 2010 source code (http://fsharppowerpack.codeplex.com/SourceControl/list/changesets) and compile the Parallel.Seq project after adding a reference to System.Threading (which was augmented by Rx), then 3) reference the dll in my ProjectEuler project and open Microsoft.FSharp.Collections.

Problem 10 is a perfect candidate for parallelization, since addition is commutative.  The following parallelized version of Problem 10 literally runs twice as fast on my dual core processor than the (nearly) identical non-parallelized version (b).

let problem10e = 
    ({3..2..1999999}
    |> PSeq.filter isPrime
    |> PSeq.sumBy int64) + 2L

To highlight how simple and effective PLINQ is, consider the following attempt at manually parallelizing using F# async workflows.  Because of the overhead involved in trying to parallelize 2 million tasks in one batch, this version actually runs 3 times slower than version (b).

let problem10c =
    let asyncIsPrime i = async { if i |> isPrime then return i else return 0 }
    let tasks = seq {for i in 3..2..1999999 -> asyncIsPrime i}
    ((Async.RunSynchronously (Async.Parallel tasks)) |> Seq.filter((<) 0) |> Seq.sumBy(int64)) + 2L

Problem 25: What is the first term in the Fibonacci sequence to contain 1000 digits?

The key to this one, unlike the traditional non-tail-recursive definition of the Fibonacci sequence, is to generate the next element in the sequence given the previous two. 

let problem25a =
    let rec find n1 n2 term =
        let n0 = n1 + n2
        if (n0 |> string).Length = 1000 then term
        else find n0 n1 (term+1)
    find 1I 1I 3

Playing around with recursive sequences, I came up with this neat solution which runs on average the same as solution (a).

let problem25b =
    let fibseq =    
        let rec fibseq n1 n2 = 
            seq {
                let n0 = n1 + n2 
                yield n0
                yield! fibseq n0 n1 }
        seq {yield 1I ; yield 1I ; yield! (fibseq 1I 1I)}
                
    (fibseq |> Seq.findIndex (fun i -> (i |> string).Length = 1000)) + 1

Problem 24: What is the millionth lexicographic permutation of the digits 0, 1, 2, 3, 4, 5, 6, 7, 8 and 9?

Completing this one marks my advance to “level 1”, having completed 25 problems (I completed Problem 67 when I completed Problem 18), and placing me in the “above what 80.11% of members have failed to do” category.

The algorithm itself was learned from http://stackoverflow.com/questions/352203/generating-permutations-lazily, which produces permutations lexicographically in sequence.  Of course, the example given there is in C++, so my F# solution is much different.  My first attempt was horrid.  After several iterations I finally got it to perform well (about 90 milliseconds, from about 6 seconds) and to at least look kind of functional.

let problem24g =
    let advance (perm:int[]) =
        //Find the longest "tail" that is ordered in decreasing order ((s+1)..perm.Length-1).
        let rec find i =
            if perm.[i] >= perm.[i-1] then i-1
            else find (i-1)
        let s = find (perm.Length-1)
        let s' = perm.[s]
        
        //Change the number just before the tail (s') to the smallest number bigger than it in the tail (perm.[t]).
        let rec find i imin =
            if i = perm.Length then imin
            elif perm.[i] > s' && perm.[i] < perm.[imin] then find (i+1) i
            else find (i+1) imin
        let t = find (s+1) (s+1)
        
        perm.[s] <- perm.[t]
        perm.[t] <- s'
        
        //Sort the tail in increasing order.
        System.Array.Sort(perm, s+1, perm.Length - s - 1)
    
    let perm = [|0;1;2;3;4;5;6;7;8;9|]
    let rec find i =
        if i = 1000000 then perm
        else advance perm ; find (i+1) 
    find  1

Problem 23: Find the sum of all the positive integers which cannot be written as the sum of two abundant numbers.

Since we’ll be reusing the sigma-divisor function again, we’ll generalize it and move it to Calc.Generalize.  Along with that I have added a performance tuned cfactorize, which returns the canonical factorization of a number (e.g. cfactorizeI 24I = [(2I, 3); (3I, 1)], noting that the second component representing the exponent is always int32, while the first is generic).  Note that the apparent silliness of breaking out the numerator and denomenator in sigma is in order to force a “sane” generic type signature.

let inline cfactorize n (lp:Lp.lp<'a>) =
    if n = lp.one then []
    else
        let flist = factorize n lp
        let rec build (f::flist) count cflist =
            if flist = [] then
                (f,count+1)::cflist
            elif f = flist.Head then
                build flist (count + 1) cflist
            else 
                build flist 0 ((f,count+1)::cflist)
        build flist 0 []

let inline sigma n (lp:Lp.lp<'a>) =
    cfactorize n lp
    |> Seq.map 
        (fun (p,e) -> 
            let numer:'a = (pown p (e+1) - lp.one)
            let denom:'a = (p-lp.one)
            numer / denom)
    |> Seq.fold (*) lp.one

Now for our solutions.  This problem definitely calls for arrays and mutation, but it mixes gracefully with other functional techniques.  My first solution (b) suffered from poor performance due to the linear time Remove operations.

let problem23a =
    let isAbundant n = (psigma n) > n
    let abundants = {1..28122} |> Seq.filter isAbundant |> Seq.toArray
    let cannots = System.Collections.Generic.List({1..28122})    
    for i in 0..(abundants.Length-1) do
        for j in i..(abundants.Length-1) do
            ignore (cannots.Remove(abundants.[i] + abundants.[j]))
    cannots |> Seq.sum

I then attempted optimize by breaking out of inner iteration when the abundant sums began to exceed 28122.  But that didn’t make much a dent.  However, by removing the linear time Remove operations and instead zeroing out  elements by array index, speed was improved by 100-fold (from about a a minute to about two seconds).

let problem23b =
    let isAbundant n = (psigma n) > n
    let abundants = {1..28122} |> Seq.filter isAbundant |> Seq.toArray
    let cannots = {1..28122} |> Seq.toArray
    for i in 0..(abundants.Length-1) do
        let rec removeAbundants j =
            let can = abundants.[i] + abundants.[j]
            if can > 28122 then ()
            else
                cannots.[can-1] <- 0 
                removeAbundants (j+1)
        removeAbundants i
    cannots |> Seq.sum

Problem 22: What is the total of all the name “scores” in the file?

Due to the nature of this problem, we don’t shy away from mutation for performance gains.  For example, sorting the names array in place.  The most novel part of this implementation is calculating uppercase alphabetical position based on character to integer conversions.

let problem22b =
    let text = System.IO.File.ReadAllText("names.txt")
    
    let names = text.Split([|',';'"'|], System.StringSplitOptions.RemoveEmptyEntries) 
    names |> Array.sortInPlace

    seq { for i in 0..names.Length-1 do 
              yield (i+1) * (names.[i] |> Seq.sumBy (fun c -> (int c) - 64)) }
    |> Seq.sum

Problem 21: Evaluate the sum of all the amicable numbers under 10000.

First, using our fast factorize function, we implement the well-known sigma-divisor function and it’s proper companion.

let sigma n =
    factorize n
    |> Seq.groupBy id
    |> Seq.map (fun (f,s) -> (f, Seq.length s))
    |> Seq.map (fun (p,e) -> (pown p (e+1) - 1)/(p-1))
    |> Seq.fold (*) 1
    
let psigma n = sigma n - n

Then the solution is straight forward, solution (a) recursively builds a list of amicable pairs, while solution (b) uses a filter with a sequence expression.  Notice the small optimizations skipping primes.

let problem21a =
    let rec build a alist =
        match (a,psigma a) with
        | (a,b) when a <> 1 && a <> b && (psigma b) = a -> build (a+1) (a::alist)
        | (9999,_) -> alist
        | _ -> build (a+1) alist
    (build 4 []) |> List.sum
    
let problem21b =
    let isAmicable a = 
        let b = psigma a
        a <> 1 && a <> b && (psigma b) = a    
    {4..9999} |> Seq.filter isAmicable |> Seq.sum

LanguagePrimitives and Generic Math

Over the past few days I’ve been exploring ways to organize and reuse common functions which repeat across Euler problems.  Something I’ve been focusing on is leveraging F#’s structural typing feature to generalized implementations as much as possible (e.g. implement isPrime once so that it can be reused for int, long and bigint).  The following design has flowed out of this discussion: http://stackoverflow.com/questions/2840714/f-static-member-type-constraints.

The following is an extension of the LanguagePrimitives module which includes a set of functions, types, and values that help us leverage F#’s ability to statically inferred structural type constraints.

namespace Microsoft.FSharp.Core

module LanguagePrimitives =

    let inline zero_of (target:'a) : 'a = LanguagePrimitives.GenericZero<'a>
    let inline one_of (target:'a) : 'a = LanguagePrimitives.GenericOne<'a>
    let inline two_of (target:'a) : 'a = one_of(target) + one_of(target)
    let inline three_of (target:'a) : 'a = two_of(target) + one_of(target)
    let inline negone_of (target:'a) : 'a = zero_of(target) - one_of(target)
    
    let inline any_of (target:'a) (x:int) : 'a = 
        let one:'a = one_of target
        let zero:'a = zero_of target
        let xu = if x > 0 then 1 else -1
        let gu:'a = if x > 0 then one else zero-one

        let rec get i g = 
            if i = x then g
            else get (i+xu) (g+gu)
        get 0 zero 
        
    type G<'a> = {
        negone:'a
        zero:'a
        one:'a
        two:'a
        three:'a
        any: int -> 'a
    }    

    let inline G_of (target:'a) : (G<'a>) = {
        zero = zero_of target
        one = one_of target
        two = two_of target
        three = three_of target
        negone = negone_of target
        any = any_of target
    }

    let gn = G_of 1   //int32
    let gL = G_of 1L  //int64
    let gI = G_of 1I  //bigint
    let gF = G_of 1.0 //float 
    let gM = G_of 1.0M//decimal

And a new module Numerics, together with nested module Generic with the help of our LanguagePrimitives extensions, allows us to implement generic numeric algorithms, and then expose specific (more efficient) partial applications of them.

module Numerics =
    open Microsoft.FSharp.Core.LanguagePrimitives

    module Generic =
        //prime factorization, from greatest to least
        let inline factorize (g:G<'a>) n = 
            let rec factorize n j flist =  
                if n = g.one then flist 
                elif n % j = g.zero then factorize (n/j) j (j::flist) 
                else factorize n (j + g.one) (flist) 
            factorize n g.two []    
        //...

    open Generic

    let inline factorizeG n = factorize (G_of n) n
    let factorizeL = factorize gL
    let factorizeI = factorize gI
    let factorize = factorize gn

    //... 

Here is an example of the kind of types signatures we get using this strategy.

val inline factorizeG :
   ^a ->  ^a list
    when  ^a : (static member get_Zero : ->  ^a) and
          ^a : (static member get_One : ->  ^a) and
          ^a : (static member ( + ) :  ^a *  ^a ->  ^a) and
          ^a : (static member ( - ) :  ^a *  ^a ->  ^a) and
          ^a : (static member ( / ) :  ^a *  ^a ->  ^a) and
          ^a : (static member ( % ) :  ^a *  ^a ->  ^a) and  ^a : equality

Problem 20: Find the sum of digits in 100!

Reusing the fact function we slipped into Problem 15, this is essentially the same as Problem 16.

let problem20a =
    fact 100I |> string |> Seq.sumBy (fun c -> c |> string |> int)

Problem 19: How many Sundays fell on the first of the month during the twentieth century?

Relieved .NET DateTime was sufficient for this.

let problem19a = 
    let mutable cur = new System.DateTime(1901, 1, 1);
    let last = new System.DateTime(2000, 12, 31);
    let mutable count = 0
    while(cur <= last) do
        if cur.Day = 1 && cur.DayOfWeek = System.DayOfWeek.Sunday then
            count <- count + 1
            
        cur <- cur.AddDays(1.0)
            
    count

And here’s a functional version.

let problem19b =
    let last = System.DateTime(2000, 12, 31)
    let rec sundays cur count = 
        if cur > last then count
        elif cur.Day = 1 && cur.DayOfWeek = System.DayOfWeek.Sunday then sundays (cur.AddDays(1.0)) (count+1)
        else sundays (cur.AddDays(1.0)) count
    sundays (System.DateTime(1901, 1, 1)) 0

Problem 18: Find the maximum total from top to bottom of the triangle below.

Wow.  The algorithm just took a little bit of pacing to realize, but I’ve gone through several refinements and variations of the implementation.

First our triangle as a string for parsing.

let tstr = "75
95 64
17 47 82
18 35 87 10
20 04 82 47 65
19 01 23 75 03 34
88 02 77 73 07 63 67
99 65 04 28 06 16 70 92
41 41 26 56 83 40 80 70 33
41 48 72 33 47 32 37 16 94 29
53 71 44 65 25 43 91 52 97 51 14
70 11 33 28 77 73 17 78 39 68 17 57
91 71 52 38 17 14 91 43 58 50 27 29 48
63 66 04 68 89 53 67 30 73 16 69 87 40 31
04 62 98 27 23 09 70 98 73 93 38 53 60 04 23"

Our first solution is imperative.

let problem18c = 
    //parse
    let tarr = Array2D.create 15 16 0 //make first column filled with zeros so no out of bounds checks required
    tstr.Split('\n') |> Array.iteri (fun i row -> row.Split(' ') |> Array.iteri (fun j cell -> (tarr.[i,j+1] <- (cell |> int))))

    //calc
    for i in 1..14 do //start with second row
        for j in 1..i+1 do //shift orientation to right
            tarr.[i,j] <- (max tarr.[i-1,j-1] tarr.[i-1,j]) + tarr.[i,j]

    //find largest        
    let mutable largest = 0        
    for j in 0..14 do 
        largest <- max largest tarr.[14,j] 
           
    largest

Our second solution uses a nicer recursive functional algorithm for computing the answer, but uses (nearly) the same parsing algorithm and array.  It would be nicer to work with a tree structure, but string.Split is just too easy.

let problem18d = 
    let tarr = Array2D.create 15 15 0
    tstr.Split('\n') |> Array.iteri (fun i row -> row.Split(' ') |> Array.iteri (fun j cell -> (tarr.[i,j] <- (cell |> int))))
    
    let rec find (row,col) = 
        if row = 14 then tarr.[row,col]
        else (max (find (row+1,col)) (find (row+1,col+1))) + tarr.[row,col]
    find (0,0)

Problem 17: If all the numbers from 1 to 1000 (one thousand) inclusive were written out in words, how many letters would be used?

After simplifying, and taking care to spell “forty” correctly.

let problem17d =
    "OneThousand".Length +
    
    "OneHundred".Length +
    "TwoHundred".Length +
    "ThreeHundred".Length +
    "FourHundred".Length +
    "FiveHundred".Length +
    "SixHundred".Length +
    "SevenHundred".Length +
    "EightHundred".Length +
    "NineHundred".Length +
    
    99*("OneHundredAnd".Length +
        "TwoHundredAnd".Length +
        "ThreeHundredAnd".Length +
        "FourHundredAnd".Length +
        "FiveHundredAnd".Length +
        "SixHundredAnd".Length +
        "SevenHundredAnd".Length +
        "EightHundredAnd".Length +
        "NineHundredAnd".Length) +
    
    10*("OneTwoThreeFourFiveSixSevenEightNineTenElevenTwelveThirteenFourteenFifteenSixteenSeventeenEighteenNineteen".Length +
        10*("Twenty".Length + "Thirty".Length + "Forty".Length + "Fifty".Length + "Sixty".Length + "Seventy".Length + "Eighty".Length + "Ninety".Length) +
        8*("OneTwoThreeFourFiveSixSevenEightNine".Length))

Problem 16: What is the sum of the digits of the number 2**1000?

Easy.

let problem16a =
    bigint.Pow(2I, 1000) |> string |> Seq.sumBy (fun c -> c |> string |> int)

Problem 15: Starting in the top left corner, how many routes (without backtracking) are there through a 20 by 20 grid to the bottom right corner?

OK.  I am incredibly pleased with myself for solving this one almost entirely through mathematical proof.  Indeed, I’ve solved it generally for NxN grids.  The only time I relied on any outside resources (http://answers.yahoo.com/question/index?qid=20100403083143AAACNJs) was to find a formula for calculating the number of multiset permutations (knowing at that point that it was exactly what I needed).  So, here goes…

Two things are clear from the statement of our problem.  First, every route has length N+N=m.  Second, it takes equally as many downward as rightward movements along the route to make it to the bottom right corner. Hence, we must find the number of permutations of elements in the multiset of length m containing m/2 downward movements and m/2 rightward movements.

let problem15a =
    let fact m = {1I..m} |> Seq.fold (*) (1I)
    let multiset m = (fact m) / ((fact (m/2I))*(fact (m/2I)))
    multiset 40I

Problem 14: Which starting number, under one million, produces the longest “Collatz” chain?

I coded  this solution in just a couple of minutes, but ran into trouble with F#’s unchecked Int32 operations.  We can do checked operations using Checked.(+) for + for example.  After changing to Int64, this runs in about 4 and a half seconds.

let problem14a = 
    let collatz n =
        let rec collatz n count =
            let count = count + 1L
            if n = 1L then count
            elif n % 2L = 0L then collatz (n/2L) count
            else collatz ((3L*n)+1L) count
        collatz n 0L
        
    {1L..999999L} |> Seq.map (fun i -> (i, collatz i)) |> Seq.maxBy snd |> fst

Since subsequences are frequently shared between inputs of n, this is a good opportunity to apply a memoization technique.

let problem14b = 
    let mem = new System.Collections.Generic.Dictionary<int64,int64>()
    let collatz m =
        let rec collatz n count =
            let count = count + 1L
            if mem.ContainsKey n then
                let count = mem.[n] + count - 1L
                mem.Add(m, count)
                count
            elif n = 1L then 
                mem.Add(m, count)
                count
            elif n % 2L = 0L then collatz (n/2L) count
            else collatz ((3L*n)+1L) count
        collatz m 0L
        
    {1L..999999L} |> Seq.map (fun i -> (i, collatz i)) |> Seq.maxBy snd |> fst

That ran in about 1 and a half seconds.

Problem 13: Work out the first ten digits of the sum of the following one-hundred 50-digit numbers.

With F#’s bigint, this is trivial.

let problem13a =    
    [37107287533902102798797998220837590246510135740250I;
    46376937677490009712648124896970078050417018260538I;
    74324986199524741059474233309513058123726617309629I;
    91942213363574161572522430563301811072406154908250I;
    23067588207539346171171980310421047513778063246676I;
    89261670696623633820136378418383684178734361726757I;
    28112879812849979408065481931592621691275889832738I;
    44274228917432520321923589422876796487670272189318I;
    47451445736001306439091167216856844588711603153276I;
    70386486105843025439939619828917593665686757934951I;
    62176457141856560629502157223196586755079324193331I;
    64906352462741904929101432445813822663347944758178I;
    92575867718337217661963751590579239728245598838407I;
    58203565325359399008402633568948830189458628227828I;
    80181199384826282014278194139940567587151170094390I;
    35398664372827112653829987240784473053190104293586I;
    86515506006295864861532075273371959191420517255829I;
    71693888707715466499115593487603532921714970056938I;
    54370070576826684624621495650076471787294438377604I;
    53282654108756828443191190634694037855217779295145I;
    36123272525000296071075082563815656710885258350721I;
    45876576172410976447339110607218265236877223636045I;
    17423706905851860660448207621209813287860733969412I;
    81142660418086830619328460811191061556940512689692I;
    51934325451728388641918047049293215058642563049483I;
    62467221648435076201727918039944693004732956340691I;
    15732444386908125794514089057706229429197107928209I;
    55037687525678773091862540744969844508330393682126I;
    18336384825330154686196124348767681297534375946515I;
    80386287592878490201521685554828717201219257766954I;
    78182833757993103614740356856449095527097864797581I;
    16726320100436897842553539920931837441497806860984I;
    48403098129077791799088218795327364475675590848030I;
    87086987551392711854517078544161852424320693150332I;
    59959406895756536782107074926966537676326235447210I;
    69793950679652694742597709739166693763042633987085I;
    41052684708299085211399427365734116182760315001271I;
    65378607361501080857009149939512557028198746004375I;
    35829035317434717326932123578154982629742552737307I;
    94953759765105305946966067683156574377167401875275I;
    88902802571733229619176668713819931811048770190271I;
    25267680276078003013678680992525463401061632866526I;
    36270218540497705585629946580636237993140746255962I;
    24074486908231174977792365466257246923322810917141I;
    91430288197103288597806669760892938638285025333403I;
    34413065578016127815921815005561868836468420090470I;
    23053081172816430487623791969842487255036638784583I;
    11487696932154902810424020138335124462181441773470I;
    63783299490636259666498587618221225225512486764533I;
    67720186971698544312419572409913959008952310058822I;
    95548255300263520781532296796249481641953868218774I;
    76085327132285723110424803456124867697064507995236I;
    37774242535411291684276865538926205024910326572967I;
    23701913275725675285653248258265463092207058596522I;
    29798860272258331913126375147341994889534765745501I;
    18495701454879288984856827726077713721403798879715I;
    38298203783031473527721580348144513491373226651381I;
    34829543829199918180278916522431027392251122869539I;
    40957953066405232632538044100059654939159879593635I;
    29746152185502371307642255121183693803580388584903I;
    41698116222072977186158236678424689157993532961922I;
    62467957194401269043877107275048102390895523597457I;
    23189706772547915061505504953922979530901129967519I;
    86188088225875314529584099251203829009407770775672I;
    11306739708304724483816533873502340845647058077308I;
    82959174767140363198008187129011875491310547126581I;
    97623331044818386269515456334926366572897563400500I;
    42846280183517070527831839425882145521227251250327I;
    55121603546981200581762165212827652751691296897789I;
    32238195734329339946437501907836945765883352399886I;
    75506164965184775180738168837861091527357929701337I;
    62177842752192623401942399639168044983993173312731I;
    32924185707147349566916674687634660915035914677504I;
    99518671430235219628894890102423325116913619626622I;
    73267460800591547471830798392868535206946944540724I;
    76841822524674417161514036427982273348055556214818I;
    97142617910342598647204516893989422179826088076852I;
    87783646182799346313767754307809363333018982642090I;
    10848802521674670883215120185883543223812876952786I;
    71329612474782464538636993009049310363619763878039I;
    62184073572399794223406235393808339651327408011116I;
    66627891981488087797941876876144230030984490851411I;
    60661826293682836764744779239180335110989069790714I;
    85786944089552990653640447425576083659976645795096I;
    66024396409905389607120198219976047599490197230297I;
    64913982680032973156037120041377903785566085089252I;
    16730939319872750275468906903707539413042652315011I;
    94809377245048795150954100921645863754710598436791I;
    78639167021187492431995700641917969777599028300699I;
    15368713711936614952811305876380278410754449733078I;
    40789923115535562561142322423255033685442488917353I;
    44889911501440648020369068063960672322193204149535I;
    41503128880339536053299340368006977710650566631954I;
    81234880673210146739058568557934581403627822703280I;
    82616570773948327592232845941706525094512325230608I;
    22918802058777319719839450180888072429661980811197I;
    77158542502016545090413245809786882778948721859617I;
    72107838435069186155435662884062257473692284509516I;
    20849603980134001723930671666823555245252804609722I;
    53503534226472524250874054075591789781264330331690I;]
    |> Seq.sum |> string |> Seq.take 10

Problem 12: What is the value of the first triangle number to have over five hundred divisors?

This has been the most sophisticated problem yet.

My first solution took 2 minutes to write, but I gave up on it after 5 minutes of running.

let problem12a = 
    let divisors n = seq {for i in 1..n/2 do if n % i = 0 then yield i}
    let rec find i last = 
        if last |> divisors |> Seq.length > 500 then last
        else find (i+1) (last+(i+1))
    find 1 1

Then I realized I had a fast function for building a list of prime factors from Problem 3 which I could use to generate a list of all divisors.  I abused that into giving me the answer in about 1 minute.

let problem12b =
    let divisors n =
        let rec divisors n j list = 
            if n = 1I then list
            elif n % j = 0I then divisors (n/j) 2I (list @ j::(list |> List.map (fun x -> x*j)))
            else divisors n (j + 1I) (list)
        1I::(divisors n 2I []) |> List.toSeq |> Seq.distinct
        
    let rec find i last = 
        if last |> divisors |> Seq.length > 500 then last
        else find (i+1I) (last+(i+1I))
    find 1I 1I

But then I remembered my Number Theory days.  Given a prime factorization p1**n1…p2**n2…p3**n3… the number of divisors is equal to (n1 + 1)*(n2 + 1)*(n3 + 1)*…

let problem12c =
    let factorize n =
        let rec factorize n j list = 
            if n = 1 then list
            elif n % j = 0 then factorize (n/j) 2 (j::list)
            else factorize n (j + 1) (list)
        factorize n 2 []
    
    let countDivisors n =
        let rec countDivisors flist count =
            match flist with 
            | [] -> count
            | f::list -> 
                let plist = list |> List.partition (fun i -> i = f)
                countDivisors (snd plist) (((fst plist |> List.length)+2)*count)
        countDivisors (factorize n) 1

    let rec find i last = 
        if last |> countDivisors > 500 then last
        else find (i+1) (last+(i+1))
    find 1 1

That gave me the answer in less than a second.

Here’s another version of (c) using Seq.groupBy instead of recursively apply List.partition.  It runs in about the same amount of time.

let problem12d =
    let factorize n =
        let rec factorize n j list = 
            if n = 1 then list
            elif n % j = 0 then factorize (n/j) 2 (j::list)
            else factorize n (j + 1) (list)
        factorize n 2 []
    
    let countDivisors n =
        factorize n
        |> Seq.groupBy id
        |> Seq.fold (fun acc (f,s) -> ((Seq.length s) + 1)*acc) 1

    let rec find i last = 
        if last |> countDivisors > 500 then last
        else find (i+1) (last+(i+1))
    find 1 1

Problem 11: What is the greatest product of four numbers on the same straight line in the 20 by 20 grid?

Not too hard, but it offered some frustration.  First, I couldn’t get the compiler to parse my multi-line list of lists for the grid.  Second, I initially didn’t consider diagonal lines going from top right to bottom left.  But because the question is posed more generally as “four adjacent numbers” conflicting with the title of the problem which states “four numbers on the same straight line,” I was led astray trying to solve the much more difficult adjacent problem.

Our solution exhaustively calculates all possible products, taking care not to go out of bounds of the grid.

let problem11a =
    let grid = array2D [[08;02;22;97;38;15;00;40;00;75;04;05;07;78;52;12;50;77;91;08];[49;49;99;40;17;81;18;57;60;87;17;40;98;43;69;48;04;56;62;00];[81;49;31;73;55;79;14;29;93;71;40;67;53;88;30;03;49;13;36;65];[52;70;95;23;04;60;11;42;69;24;68;56;01;32;56;71;37;02;36;91];[22;31;16;71;51;67;63;89;41;92;36;54;22;40;40;28;66;33;13;80];[24;47;32;60;99;03;45;02;44;75;33;53;78;36;84;20;35;17;12;50];[32;98;81;28;64;23;67;10;26;38;40;67;59;54;70;66;18;38;64;70];[67;26;20;68;02;62;12;20;95;63;94;39;63;08;40;91;66;49;94;21];[24;55;58;05;66;73;99;26;97;17;78;78;96;83;14;88;34;89;63;72];[21;36;23;09;75;00;76;44;20;45;35;14;00;61;33;97;34;31;33;95];[78;17;53;28;22;75;31;67;15;94;03;80;04;62;16;14;09;53;56;92];[16;39;05;42;96;35;31;47;55;58;88;24;00;17;54;24;36;29;85;57];[86;56;00;48;35;71;89;07;05;44;44;37;44;60;21;58;51;54;17;58];[19;80;81;68;05;94;47;69;28;73;92;13;86;52;17;77;04;89;55;40];[04;52;08;83;97;35;99;16;07;97;57;32;16;26;26;79;33;27;98;66];[88;36;68;87;57;62;20;72;03;46;33;67;46;55;12;32;63;93;53;69];[04;42;16;73;38;25;39;11;24;94;72;18;08;46;29;32;40;62;76;36];[20;69;36;41;72;30;23;88;34;62;99;69;82;67;59;85;74;04;36;16];[20;73;35;29;78;31;90;01;74;31;49;71;48;86;81;16;23;57;05;54];[01;70;54;71;83;51;54;69;16;92;33;48;61;43;52;01;89;19;67;48]]

    let calc (i,j) (iu,ju) =
        {0..3} |> Seq.map (fun k -> grid.[i+(k*iu),j+(k*ju)]) |> Seq.fold (*) 1

    seq {for i in 0..16 do 
             for j in 0..16 do
                 yield calc (i,j) (1,0)
                 yield calc (i,j) (0,1)
                 yield calc (i,j) (1,1)
                 if i-3 >= 0 then yield calc (i,j) (-1,1) }
    |> Seq.max

Problem 10: Find the sum of all the primes below two million.

This problem is a variation on Problem 7 allowing us to reuse isPrime and modify nthPrime to build a list of primes up to a max.

//problem 10
let problem10a =
    let primes max = 
        let rec primes i plist =
            if i > max then plist
            elif i |> isPrime then primes (i+2) (i::plist)
            else primes (i+2) plist
        primes 3 [2] //start at 3, and with a list containing our only even prime
        
    primes 1999999 |> List.sumBy int64  //convert to int64 to avoid overflow

Another solution using a sequence expression is even nicer though,

let problem10b =
    (seq { for i in 3..2..1999999 do if i |> isPrime then yield i} //odd primes
    |> Seq.sumBy int64) + 2L //dont forget 2

Problem 9: There exists exactly one Pythagorean triplet for which a + b + c = 1000. Find the product abc.

Here we generate a sequence of all possible triplets meeting the conditions, and find the first which is Pythagorean.  It took a while to get this to perform well, my original solution required three nested loops and took several seconds to compute.  By playing with the conditions a < b < c and a + b + c = 1000, and knowing the first Pythagorean triplet is (3,4,5), I was able to apply some bounds to a and b (though the upper bound for a will never trip when it’s piped into Seq.find).  This solution takes a fraction of a second.

let problem9a =
    seq {for a in 3..332 do for b in 4..498 do yield (a, b, 1000 - a - b)}
    |> Seq.find (fun (a,b,c) -> (pown a 2) + (pown b 2) = (pown c 2))
    |> (fun (a,b,c) -> a*b*c)