60 Days of Euler in F# - Problem 46

The Problem

It was proposed by Christian Goldbach that every odd composite number can be written as the sum of a prime and twice a square.

  • 9 = 7 + 2×12
  • 15 = 7 + 2×22
  • 21 = 3 + 2×32
  • 25 = 7 + 2×32
  • 27 = 19 + 2×22
  • 33 = 31 + 2×12

It turns out that the conjecture was false.

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

The Solution

As we do with some many of these problems, we need some primes. The problem definition asks us to find odd composite numbers that meet the criteria. In this case "odd composite" is simply a way of saying "the multiple of two primes other than 2".

We'll start by stealing the getSomePrimes function from Problem 7. This function generates a specified number of primes.

let intsqrt i = int(sqrt(float i))

let isPrime i =
    if i <= 1 then false
    elif i = 2 then true
    elif (i &&& 1) = 0 then false
    else
        let sr = intsqrt i
        seq { 3..2..sr } |> Seq.forall (fun f -> i%f<>0)

let rec nextPrime i =
    match i with
    | i when isPrime i -> i
    | i -> nextPrime (i + 1)

let getSomePrimes num =
    let primeGenerator (candidate,count) =
        if count > num then
            None
        elif candidate = 2 then
            Some(2,(3,count+1))
        else
            let next = nextPrime candidate
            Some(next,(next+2,count+1))

    Seq.unfold primeGenerator (2,1)

The plan of attack is to multiply primes together to generate the composites. We also need a list of primes to add squares to to test for solutions. How many primes do we need? I have no idea and lack the mathematical analysis capabilities to figure it out. So I think 1000 is a good place to start. We can adjust later if we need to. Because computing primes is expensive, I'm storing them in a List so they only get generated once.

let primes = 
    getSomePrimes 1000 
    |> Seq.skip 1
    |> List.ofSeq

primes will contain a list of 999 primes. It skips over 2 (the first element in the list) so that we don't generate any even numbers.

Now we need some squares to add to our primes.

let generateSquares state =
    Some(state * state, state + 1)

let squares =
    Seq.unfold generateSquares 1
    |> Seq.take 1000
    |> List.ofSeq

We use Seq.unfold to generate squares for the first 1000 integers. Again, 1000 is just a number I picked to get started. We can adjust later if we need to.

Now we need a way to test possible solutions:

let inline isMatch p s n = p + s*2 = n

let findSquare p n = 
    squares 
    |> Seq.filter (fun s -> isMatch p s n)
    |> Seq.tryHead

let solution n =
    primes
    |> Seq.filter (fun p -> p < n)
    |> Seq.map (fun p -> (p,findSquare p n))
    |> Seq.filter (fun (_,square) -> square.IsSome)
    |> Seq.tryHead

The isMatch function takes a prime p, a square s, and a composite number n as input. If these inputs satisfy the conjecture, it returns true.

findSquare takes a prime p and a composite number n as input. It returns the first square that satisfies the conjecture; otherwise it returns None.

The solution function takes a composite number n as input. It takes our list of primes, filters out any below n, and tries to findSquare that satisfies the conjecture with that prime. If a solution is found, it returns the solution; otherwise it returns None.

Now we can solve the problem:

let crossJoin xs ys = 
    xs
    |> Seq.collect (fun x -> ys |> Seq.map (fun y -> x, y))

#time
crossJoin primes primes
|> Seq.map (fun (a,b) -> a*b)
|> Seq.sort
|> Seq.map (fun n -> (n, solution n))
|> Seq.filter (fun (n,solution) -> solution.IsNone)
|> Seq.head
|> fst
#time

We start by cross joining all our primes with each other prime, multiplying them together to get composite numbers, and then sorting the resulting sequence. The sort step is necessary because we are looking for the smallest composite that fails to satisfy the conjecture.

Next, we use Seq.map to generate a solution for each composite. Then we filter out composites where there is no solution. The head of the resulting sequence is the answer.

Note: Above, we arbitrarily chose 1000 primes and 1000 squares. It turns out this generates enough data to find the answer. If it didn't, we'd have to experiment with increasing these numbers until we finally found the answer.

Other Posts in This Series