theburningmonk.com

Prob­lem

The small­est num­ber express­ible as the sum of a prime square, prime cube, and prime fourth power is 28. In fact, there are exactly four num­bers below fifty that can be expressed in such a way:

28 = 2² + 2³ + 2⁴

33 = 3² + 2³ + 2⁴

49 = 5² + 2³ + 2⁴

47 = 2² + 3³ + 2⁴

How many num­bers below fifty mil­lion can be expressed as the sum of a prime square, prime cube, and prime fourth power?

Solu­tion

// generate all prime numbers under <= this max
let max = int64(sqrt(double(50000000L)))

// initialise the list with 2 which is the only even number in the sequence
let mutable primeNumbers = [2L]

// only check the prime numbers which are <= the square root of the number n
let hasDivisor n =
    primeNumbers
    |> Seq.takeWhile (fun n' -> n' <= int64(sqrt(double(n))))
    |> Seq.exists (fun n' -> n % n' = 0L)

// only check odd numbers <= max
let potentialPrimes = Seq.unfold (fun n -> if n > max then None else Some(n, n+2L)) 3L

// populate the prime numbers list
for n in potentialPrimes do if not(hasDivisor n) then primeNumbers <- primeNumbers @ [n]

// use the same hasDivisor function instead of the prime numbers list as it offers
// far greater coverage as the number n is square rooted so this function can
// provide a valid test up to max*max
let isPrime n = if n = 1L then false else not(hasDivisor(n))

let answer =
    primeNumbers
    |> Seq.collect (fun n ->
        primeNumbers 
        |> Seq.map (fun n' -> pown n 2 + pown n' 3) 
        |> Seq.takeWhile (fun sum -> sum < 50000000L))
    |> Seq.collect (fun sum ->
        primeNumbers 
        |> Seq.map (fun n -> sum + pown n 4) 
        |> Seq.takeWhile (fun sum' -> sum' < 50000000L))
    |> Seq.distinct
    |> Seq.length

The biggest prime that can appear in the equa­tion is equals to the square root of 50000000 – 8 – 16, which, inci­den­tally is just over 7000, so the num­bers of primes involved is rea­son­ably small which bodes well for a fast solution!

The logic in this solu­tion is oth­er­wise sim­ple, using the cached prime num­bers list to first gen­er­ate num­bers (< 50 mil­lion) that can be writ­ten as the sum of a prime square and prime cube; then for each num­ber see if we can add a prime fourth power and get a sum less than 50 million.

One thing that caught me out ini­tially was the need to search for dis­tinct num­bers as some of these num­bers do over­lap, but oth­er­wise this is a solu­tion which runs hap­pily under 3 seconds.

Prob­lem

A num­ber chain is cre­ated by con­tin­u­ously adding the square of the dig­its in a num­ber to form a new num­ber until it has been seen before.

For exam­ple,

There­fore any chain that arrives at 1 or 89 will become stuck in an end­less loop. What is most amaz­ing is that EVERY start­ing num­ber will even­tu­ally arrive at 1 or 89.

How many start­ing num­bers below ten mil­lion will arrive at 89?

Solu­tion

let max = 10000000

// build up a cache for all the numbers from 1 to 10 million
let cache = Array.init (max+1) (fun n ->
    match n with
    | 0 | 1 -> Some(false)
    | 89 -> Some(true)
    | _ -> None)

// define function to add the square of the digits in a number
let addSquaredDigits n =
    n.ToString().ToCharArray()
    |> Array.sumBy (fun c -> pown (int(c.ToString())) 2)

// define function to take an initial number n and generate its number chain until
// it gets to a number whose subsequent chain ends with 1 or 89, which means that
// all previous numbers will also end in the same number
let processChain n =
    let rec processChainRec n (list: int list) =
        if cache.[n] = None then processChainRec (addSquaredDigits n) (list@[n])
        else list |> List.iter (fun n' -> cache.[n'] <- cache.[n])
    processChainRec n []

// go through all the numbers from 2 to 10 million using the above function
[2..10000000] |> List.iter processChain

// how many numbers whose number chain ends with 89?
let answer = cache |> Array.filter (fun n -> n = Some(true)) |> Array.length

This is actu­ally a fairly sim­ple prob­lem, with the main chal­lenge being how to make it run fast enough to obey with the one minute rule which is where the cache comes in.

Using the prop­erty that if a chain start­ing with the num­ber n, i.e. n, n1, n2, …, ends in 89 then the chain start­ing with any of the num­bers n1, n2, … will also end in 89 we can min­imise the amount of com­pu­ta­tion required by caching pre­vi­ous results. Note I didn’t have to mark the cache as muta­ble in this case because the Array con­struct in F# is a fixed-sized, zero-index, muta­ble col­lec­tions of elements.

Prob­lem

The prime 41, can be writ­ten as the sum of six con­sec­u­tive primes:

41 = 2 + 3 + 5 + 7 + 11 + 13

This is the longest sum of con­sec­u­tive primes that adds to a prime below one-hundred.

The longest sum of con­sec­u­tive primes below one-thousand that adds to a prime, con­tains 21 terms, and is equal to 953.

Which prime, below one-million, can be writ­ten as the sum of the most con­sec­u­tive primes?

Solu­tion

// generate all prime numbers under <= this max
let max = 10000

// initialise the list with 2 which is the only even number in the sequence
let mutable primeNumbers = [2]

// only check the prime numbers which are <= the square root of the number n
let hasDivisor n =
    primeNumbers
    |> Seq.takeWhile (fun n' -> n' <= int(sqrt(double(n))))
    |> Seq.exists (fun n' -> n % n' = 0)

// only check odd numbers <= max
let potentialPrimes = Seq.unfold (fun n -> if n > max then None else Some(n, n+2)) 3

// populate the prime numbers list
for n in potentialPrimes do
    if not(hasDivisor n) then primeNumbers <- primeNumbers @ [n]

// use the same hasDivisor function instead of the prime numbers list as it offers
// far greater coverage as the number n is square rooted so this function can
// provide a valid test up to max*max
let isPrime n = if n = 1 then false else not(hasDivisor(n))

// define function which computes the sum of the all consecutive primes starting
// from n, and returns the longest sequence which sums to a prime
let getPrimeSequence n =
    primeNumbers
    |> Seq.filter (fun n' -> n' > n)
    |> Seq.scan (fun (sum, count) n' -> (sum+n', count+1)) (n, 1)
    |> Seq.takeWhile (fun (sum, count) -> sum < 1000000)
    |> Seq.filter (fun (sum, count) -> isPrime sum)
    |> Seq.maxBy (fun (sum, count) -> count)

// for all numbers, find the longest sequence
let answer = primeNumbers |> Seq.map getPrimeSequence |> Seq.maxBy (fun (sum, count) -> count)

Again, uti­liz­ing the PGsim­ple 3 algo­rithm from here, this solu­tion runs in just over a sec­ond as I’ve restricted it to check sequences start­ing with a prime < 10000 which per­haps unsur­pris­ingly, has pro­vided me with the answer.

The only tricky part in this prob­lem is how to build up a list of sequences for each start­ing prime and find the longest sequence start­ing from that num­ber which results in a prime < 1000000. For­tu­nately you can do this rather eas­ily with Seq.scan in F# which works sim­i­larly to Seq.fold but gives you all the inter­me­di­ate results as well as the final result. Used in con­junc­tion with Seq.takeWhile it allows me to get all the sequences which sums up to a num­ber less than 1000000. The rest is pretty straight for­ward, fil­ter out the sums that aren’t prime num­bers and find the longest sequence, and do this for all the primes we have gen­er­ated to find the answer to the problem.

Prob­lem

It is pos­si­ble to show that the square root of two can be expressed as an infi­nite con­tin­ued fraction.

By expand­ing this for the first four iter­a­tions, we get:

1 + 1/2 = 3/2 = 1.5

1 + 1/(2 + 1/2) = 7/5 = 1.4

1 + 1/(2 + 1/(2 + 1/2)) = 17/12 = 1.41666…

1 + 1/(2 + 1/(2 + 1/(2 + 1/2))) = 41/29 = 1.41379…

The next three expan­sions are 99/70, 239/169, and 577/408, but the eighth expan­sion, 1393/985, is the first exam­ple where the num­ber of dig­its in the numer­a­tor exceeds the num­ber of dig­its in the denominator.

In the first one-thousand expan­sions, how many frac­tions con­tain a numer­a­tor with more dig­its than denominator?

Solu­tion

// define function to return all the numerator-denominator pairs for the first n expand
let expand n =
    Seq.unfold (fun (num, denom) -> Some((num, denom), (denom*2I+num, denom+num))) (3I, 2I)
    |> Seq.take n

let answer =
    expand 1000
    |> Seq.filter (fun (num, denom) -> num.ToString().Length > denom.ToString().Length)
    |> Seq.length

If you look at the pat­terns 3/2, 7/5, 17/12, 41/29, and so on, it’s easy to spot a pat­tern where the numer­a­tor and denom­i­na­tor of iter­a­tion n can be derived from the iter­a­tion n-1:

Numerator(n) = Numerator(n-1) + 2 * Denominator(n-1)

Denominator(n) = Numerator(n-1) + Denominator(n-1)

Prob­lem

Start­ing with 1 and spi­ralling anti­clock­wise in the fol­low­ing way, a square spi­ral with side length 7 is formed.

It is inter­est­ing to note that the odd squares lie along the bot­tom right diag­o­nal, but what is more inter­est­ing is that 8 out of the 13 num­bers lying along both diag­o­nals are prime; that is, a ratio of 8/13 62%.

If one com­plete new layer is wrapped around the spi­ral above, a square spi­ral with side length 9 will be formed. If this process is con­tin­ued, what is the side length of the square spi­ral for which the ratio of primes along both diag­o­nals first falls below 10%?

Solu­tion

let hasDivisor(n) =
   let upperBound = int(sqrt(double(n)))
   [2..upperBound] |> Seq.exists (fun x -> n % x = 0)

let isPrime(n) = if n = 1 then false else not(hasDivisor(n))

// define function that returns the number on the corners of a spiral of length n
let getCornerNumbers n =
    match n with
    | 1 -> [1]
    | _ when n % 2 = 0 -> []
    | _ -> [3..-1..0] |> List.map (fun n' -> n*n - n'*(n-1))

let answer =
    let mutable cornerNumbers, primeNumbers, size = 0, 0, 1
    let mutable continueLoop = true

    while continueLoop do
        // get the numbers that appear at the corners of a spiral of the given size
        let newNumbers = getCornerNumbers size

        // increment the totals
        cornerNumbers <- cornerNumbers + newNumbers.Length
        primeNumbers <- primeNumbers + (newNumbers |> List.filter isPrime |> List.length)

        let ratio = double(primeNumbers) / double(cornerNumbers)
        if ratio < 0.1 && size > 1 then continueLoop <- false else size <- size + 2
    size

UPDATE: Hav­ing stum­bled upon some very good algo­rithms for gen­er­at­ing prime num­ber sequence here, I decided to revisit my solu­tion here and using the PGSimple3 algo­rithm the new code now runs in seconds!

// generate all prime numbers under <= this max
let max = 100000

// initialise the list with 2 which is the only even number in the sequence
let mutable primeNumbers = [2]

// only check the prime numbers which are <= the square root of the number n
let hasDivisor n =
    primeNumbers
    |> Seq.takeWhile (fun n' -> n' <= int(sqrt(double(n))))
    |> Seq.exists (fun n' -> n % n' = 0)

// only check odd numbers <= max
let potentialPrimes = Seq.unfold (fun n -> if n > max then None else Some(n, n+2)) 3

// populate the prime numbers list
for n in potentialPrimes do
    if not(hasDivisor n) then primeNumbers <- primeNumbers @ [n]

// use the same hasDivisor function instead of the prime numbers list as it offers
// far greater coverage as the number n is square rooted so this function can
// provide a valid test up to max*max
let isPrime n = if n = 1 then false else not(hasDivisor(n))

// define function that returns the number on the corners of a spiral of length n
let getCornerNumbers n =
    match n with
    | 1 -> [1]
    | _ when n % 2 = 0 -> []
    | _ -> [3..-1..0] |> List.map (fun n' -> n*n - n'*(n-1))

let answer =
    let mutable cornerNumbers, primeNumbers, size = 0, 0, 1
    let mutable continueLoop = true

    while continueLoop do
        // get the numbers that appear at the corners of a spiral of the given size
        let newNumbers = getCornerNumbers size

        // increment the totals
        cornerNumbers <- cornerNumbers + newNumbers.Length
        primeNumbers <- primeNumbers + (newNumbers |> List.filter isPrime |> List.length)

        let ratio = double(primeNumbers) / double(cornerNumbers)
        if ratio < 0.1 && size > 1 then continueLoop <- false else size <- size + 2
    size