Rash thoughts about .NET, C#, F# and Dynamics NAV.


"Every solution will only lead to new problems."

Monday, 8. December 2008


Project Euler in F# – Problem 53 – Dynamic Programming

Filed under: English posts,F#,Theoretische — Steffen Forkmann at 12:23 Uhr

Claudio Cherubino (from fsharp.it/) posted his solution to Euler Project – problem 53. As I dealed a lot with Dynamic Programming in the recent time, I tried to solve the problem with a dynamic program in F#.

Project Euler – Problem 53:

How many values of C(n,k), for 1 ≤ n ≤ 100, exceed one-million?

Remark: C(n,k) are the binomial coefficients.

As it turned out, this is not that complicated if one knows the recursive function for the binomial coefficients (see Wikipedia):

Recursive definition for the Binomial coefficient with Recursive definition for the Binomial coefficient

This is easily transformed into a F# program:

let binomials = Array2.create 101 101 1I
let answer = ref 0
for n in [1..100] do
  for k in [1..n - 1] do
    binomials.[n, k] <- binomials.[n - 1,k] + binomials.[n - 1,k - 1]
    if binomials.[n, k] > 1000000I then
      answer := !answer + 1

!answer |> printf "Answer: %A"

Claudio’s program took 1315ms on my computer. The dynamic program needs only 63ms. But we can still do better if we use the symmetry of Pascal’s triangle.

Symmetry of Pascal's triangle

This leads to an algorithm, which calculates only half of the binomial coefficients.

let binomials = Array2.create 101 101 1I
let answer = ref 0
for n in [1..100] do
  for k in [1..n/2] do
    let b = binomials.[n - 1,k] + binomials.[n - 1,k - 1]
    binomials.[n, k] <- b
    binomials.[n, n - k] <- b
    if b > 1000000I then
      if k = n-k then
        answer := !answer + 1
      else
        answer := !answer + 2
!answer |> printf "Answer: %A"

This version needs only 45ms – but we are not ready. I mentioned Pascal’s triangle and its symmetry. But we can use another property. We don’t have to calculate the complete row, if we exceed 100000. All values behind this threshold have to be greater.

let binomials = Array2.create 101 101 1I
let answer = ref 0
for n in [1..100] do
  let threshold_reached = ref false
  let c = ref 0
  for k in [1..n/2] do
    if not !threshold_reached then
      let b = binomials.[n - 1,k] + binomials.[n - 1,k - 1]
      binomials.[n, k] <- b
      binomials.[n, n - k] <- b
      if b > 1000000I then
        threshold_reached := true
      else
        c := !c + 1

  if !threshold_reached then
    answer := !answer + (n - 1) - (2 * !c)
!answer |> printf "Answer: %A"

This final version took only 29 ms.

In the next posting I will show a version using memoization.

Tags: , , , ,

6 Comments »

  1. Great efficiency improvements! I wrote the following version of your algorithm, using recursion instead of an imperative style and also using (k-1) instead of c. My version seems to be a little more efficient.

    let euler53() =
        let binomial (* n k *) =
            let binomials = Array2.create 101 101 1I
            fun n k -> 
                let b = binomials.[n-1, k] + binomials.[n-1, k-1]
                binomials.[n, k] <- b
                binomials.[n, n-k]  100 then acc 
                                         else kLoop n 1 acc
        and kLoop n k acc =
            if k > n/2 
                then nLoop (n+1) acc
            elif binomial n k > 1000000I 
                then nLoop (n+1) (acc + n-1 - 2*(k-1))
            else kLoop n (k+1) acc                  
                          
        nLoop 1 0

    Comment by GrahameTheHornist — Tuesday, 9. December 2008 um 0:47 Uhr

  2. Shoot, no formatting! Sorry about that.

    Comment by GrahameTheHornist — Tuesday, 9. December 2008 um 0:48 Uhr

  3. Great post!
    Dynamic programming is exactly what I had in mind when I asked for optimization strategies.

    Comment by Claudio — Tuesday, 9. December 2008 um 9:09 Uhr

  4. Thanks Grahame. But I think something went wrong. Where is the definition of the nLoop?

    (i added a <pre> for code formatting)

    Comment by Steffen Forkmann — Tuesday, 9. December 2008 um 9:10 Uhr

  5. let euler53() =
        let binomial (* n k *) =
            let binomials = Array2.create 101 101 1I
            fun n k -> 
                let b = binomials.[n-1, k] + binomials.[n-1, k-1]
                binomials.[n, k] <- b
                binomials.[n, n-k] <- b
                b
            
        let rec nLoop n acc = if n > 100 then acc 
                                         else kLoop n 1 acc
        and kLoop n k acc =
            if k > n/2 
                then nLoop (n+1) acc
            elif binomial n k > 1000000I 
                then nLoop (n+1) (acc + n-1 - 2*(k-1))
            else kLoop n (k+1) acc                  
                          
        nLoop 1 0 
    

    Comment by GrahameTheHornist — Tuesday, 9. December 2008 um 17:56 Uhr

  6. @GrahameTheHornist:

    Thank you. I tried your solution on my computer and with 16ms it is the fastest so far. Well done.

    Comment by Steffen Forkmann — Tuesday, 9. December 2008 um 18:18 Uhr

RSS feed for comments on this post. | TrackBack URI

Leave a comment

XHTML ( You can use these tags): <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> .