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)

Problem 8: Find the greatest product of five consecutive digits in the 1000-digit number.

First we get our 1000-digit number into a string suitable for digit traversal.  F# has multi-line strings and the backslash character may be used at the end of each line to prevent newline characters from being inserted.

let digits1000 = "73167176531330624919225119674426574742355349194934\
96983520312774506326239578318016984801869478851843\
85861560789112949495459501737958331952853208805511\
12540698747158523863050715693290963295227443043557\
66896648950445244523161731856403098711121722383113\
62229893423380308135336276614282806444486645238749\
30358907296290491560440772390713810515859307960866\
70172427121883998797908792274921901699720888093776\
65727333001053367881220235421809751254540594752243\
52584907711670556013604839586446706324415722155397\
53697817977846174064955149290862569321978468622482\
83972241375657056057490261407972968652414535100474\
82166370484403199890008895243450658541227588666881\
16427171479924442928230863465674813919123162824586\
17866458359124566529476545682848912883142607690042\
24219022671055626321111109370544217506941658960408\
07198403850962455444362981230987879927244284909188\
84580156166097919133875499200524063689912560717606\
05886116467109405077541002256983155200055935729725\
71636269561882670428252483600823257530420752963450"

While the algorithm required for solving this problem is straight forward, it gives us a nice opportunity to compare an iterative approach (a) with a pure functional approach (b).  Note: make sure to convert each char to a string before converting to an int!

let problem8a =
    let mutable max = 0
    for i in 0..(digits1000.Length-5) do
        let mutable next = 1
        for j in 0..4 do
            next <- next * (digits1000.[i+j] |> string |> int)
        if next > max then max <- next
    max
       
let problem8b =
    {0..digits1000.Length-5}
    |> Seq.map
        (fun i -> 
             {0..4} 
             |> Seq.map (fun j -> digits1000.[i+j] |> string |> int) 
             |> Seq.fold (*) 1)
    |> Seq.max

I personally prefer the descriptiveness of (b) and it’s performance is on par with (a) .

Problem 7: What is the 10001st prime number?

I think I’m getting the hang of this!  Some notable optimizations: isPrime only needs to check the divisibility of n by numbers 2..sqrt(n), and nthPrime only needs to check odd numbers for primality by incrementing i by steps of 2 starting at 1.

let problem7a =
    let isPrime n =
        let nsqrt = n |> float |> sqrt |> int
        let rec isPrime i =
            if i > nsqrt then true
            elif n % i = 0 then false
            else isPrime (i+1)
        isPrime 2
        
    let nthPrime n = 
        let rec nthPrime i p count =
            if count = n then p
            elif i |> isPrime then nthPrime (i+2) i (count+1)
            else nthPrime (i+2) p count
        nthPrime 1 1 0
        
    nthPrime 10001

Problem 6: Find the difference between the sum of the squares of the first one hundred natural numbers and the square of the sum.

This one is simple.

let problem6a =
    let range = {1..100}
    let sumOfSquares = range |> Seq.sumBy (fun x -> x * x)
    let sumSquared = pown (range |> Seq.sum) 2
    sumSquared - sumOfSquares

Problem 5: What is the smallest number that is evenly divisible by all of the numbers from 1 to 20?

This one was pretty easy.  We only need to test i against 11 to 20 since the numbers in that range are themselves divisible by those in 1 to 10.  Also, a lower bound for i is the product of all primes in 1 to 20.

let problem5a =
    let rec find i = 
        if {11..20} |> Seq.forall (fun d -> i % d = 0) then i
        else find (i + 1)
    find (2*3*5*7*11*13*17*19)

Problem 4: Find the largest palidromic number made from the product of two 3-digit numbers.

First we implement a function for determining whether a number is palindromic.  We use an efficient tail recursive algorithm comparing the string converted number by character indices working from either end of the string until we reach the middle.

let isPalindrome n =
    let rec isPalindrome (nstr:string) i j =
        if i > j then true 
        elif  nstr.[i] <> nstr.[j] then false
        else isPalindrome nstr (i+1) (j-1)
    let nstr = n |> string
    isPalindrome nstr 0 (nstr.Length - 1)

Now we just need to enumerate all products made from two 3-digit (100 through 999) numbers, taking the max of those which are palindromic.  Solution (a) uses an imperative approach justifiable in nested loop scenarios while (b) and (c) use sequence expressions and typical function applications.

//fastest, reasonable use of iterative loop + mutable
let problem4a =        
    let mutable max = 0
    for i in 100..999 do
        for j in 100..999 do
            let p = i * j
            if p > max && p |> isPalindrome then max <- p    
    max

//notably slower since has to process isPalindrome for every possibility
let problem4b =
    seq { for i in 100..999 do for j in 100..999 do yield i * j }
    |> Seq.filter isPalindrome
    |> Seq.max

//still lower than 4a
let problem4c =
    seq { for i in 100..999 do for j in 100..999 do yield i * j }
    |> Seq.fold (fun acc i -> if i > acc && i |> isPalindrome then i else acc) 0

While (b) is the clearest, I like solution (a) the best since it’s clear enough and notably faster than either (b) or (c).

Problem 3: What is the largest prime factor of the number 600851475143?

This one required a flash of inspiration.  We recursively build a list of prime factors by finding the first even divisor j of n, necessarily prime, and then repeating with n/j until we reach n = 1.  It works out that the head of the list is the largest prime factor.

let problem3a =
    let factorize n =
        let rec factorize n j list = 
            if n = 1I then list
            elif n % j = 0I then factorize (n/j) j (j::list)
            else factorize n (j + 1I) (list)
        factorize n 2I []
    600851475143I |> factorize |> List.head

Now, we only need the largest prime factor so this algorithm can be specialized.

let problem3b =
    let rec largestFactor n j largest = 
        if n = 1I then largest
        elif n % j = 0I then largestFactor (n/j) j j
        else largestFactor n (j + 1I) largest
    largestFactor 600851475143I 2I 1I

Problem 2: Find the sum of all the even-valued terms in the Fibonacci sequence which do not exceed four million.

First we implement a function for computing terms of the Fibonacci sequence.

let rec fib n = if (n = 1 || n = 0) then 1 else fib(n-1) + fib(n-2)

We pipe that into a sequence generator, taking all even terms not exceeding four million, and then sum them.

let problem2a =
    fib
    |> Seq.initInfinite
    |> Seq.takeWhile (fun i -> i <= 4000000) 
    |> Seq.filter (fun i -> i % 2 = 0)
    |> Seq.sum

Problem 1: Add all the natural numbers below one thousand that are multiples of 3 or 5.

Solutions (a)-(c) take the same approach but use different Seq function applications. Every integer from 1 to 999 is tested.

let problem1a =
    {1..999}
    |> Seq.filter (fun i -> i % 3 = 0 || i % 5 = 0)
    |> Seq.sum

let problem1b =
    {1..999}
    |> Seq.fold (fun acc i -> if i % 3 = 0 || i % 5 = 0 then acc + i else acc) 0
    
let problem1c =    
    {1..999}
    |> Seq.sumBy (fun i -> if i % 3 = 0 || i % 5 = 0 then i else 0)

Solution (d) uses inclusion/exclusion instead of an exhaustive search.

let problem1d = 
    let a = Seq.append {0..3..999} {0..5..999} |> Seq.sum
    let b = {0..15..999} |> Seq.sum
    a - b    

I like solution (a) the best since it's the clearest.

The beginning

Back in 2008 I happened upon Expert F# while browsing my local book store.  Having a background in mathematics, I was enticed by the Introduction which pitched F# as a programming language close in expression to what mathematicians are accustomed.  I purchased the book and read it in large part while also following the language on the web.  In this first pass I mostly absorbed the concepts of functional programming but didn't practice much with the language itself.  However, in the intervening years I've been applying a lot of the functional concepts I learned to my day work in C#.  Now, with the release of F# version 2.0, I've found renewed interest and have determined to become experienced with the language itself.  The mathematical nature of the problems in Project Euler make it a particularly well-suited set of challenges for becoming familiar with F#.  While it's been done, here I will document my journey.