Archive for November, 2011

Project Euler Problem #231!

November 3rd, 2011

Problem #231 says:

The binomial coefficient \( ^{10}C_{3} = 120\).
\(120 = 2^{3} \ 3 \ 5 = 2\ 2\ 2\ 3\ 5\), and \(2 + 2 + 2 + 3 + 5 = 14\).
So the sum of the terms in the prime factorisation of \(^{10}C_{3}\) is 14.

Find the sum of the terms in the prime factorisation of \(^{20000000}C_{15000000}\).

There’s a simple formula to calculate this. I don’t remember where I found it I wrote this solution a while ago and never posted it. Uses a bunch of prime numbers and the sum of factors to get the answer without brute force.The other way I suppose one could do this is there is a definition for the binomial coefficient as a product of a bunch of stuff. If we could eliminate all rounding error that approach could lead to a working solution.

Solution provided in F# and uses a Sieve of Atkins that seems to be slightly faster than other prime number generators I’ve been using.

#light

open System
open System.Collections
open System.Collections.Generic
open System.Diagnostics
open System.Threading
open Microsoft.FSharp.Collections
open Microsoft.FSharp.Math

let Atkins limit =

    let is_prime = new Dictionary<uint64,bool>()
    let keys = [5UL..limit]
    for i in keys do
        is_prime.Add(i,false)

    let sqrt_limit = uint64 (Math.Sqrt (float limit))
    let xys = [1UL..sqrt_limit]
    for x in xys do
        for y in xys do
            let mutable n = 4UL*(x*x) + (y*y)
            if (n <= limit && (n % 12UL = 1UL || n % 12UL = 5UL)) then
                is_prime.[n] <- not is_prime.[n]
            n <- 3UL*(x*x) + (y*y)
            if (n <= limit && (n % 12UL = 7UL)) then
                is_prime.[n] <- not is_prime.[n]
            n <- 3UL*(x*x) - (y*y)
            if ((x>y) && n <= limit && (n % 12UL = 11UL)) then
                is_prime.[n] <- not is_prime.[n]

    for i in keys do
        if is_prime.[i] then
            let max = limit / (i*i)
            for k in [1UL..max] do
                is_prime.[k*(i*i)] <- false
    [
    yield 2UL;yield 3UL;
    for i in keys do
        if is_prime.[i] then
            yield i
    ]

let primes = Atkins 20000000UL |> Seq.toArray

let SumOfFactors n p  =
    let nt = ref n
    [
        while (!nt > 0UL) do
            nt := !nt / p
            yield !nt
    ] |> Seq.sum

let Problem231 n k =
    primes
    |> PSeq.map (fun i -> ((SumOfFactors n i) - (SumOfFactors k i) - (SumOfFactors (n-k) i) )* i)
    |> Seq.sum

Console.WriteLine "Primes Finished"
let sw = new Stopwatch()
sw.Start()
Console.WriteLine  (Problem231 20000000UL 15000000UL)
sw.Stop()
Console.WriteLine(sw.Elapsed)
Console.ReadLine() |> ignore

Tags: , , ,
Posted in Project Euler | Comments (0)