Project Euler 27: Quadratic Primes

Euler discovered the remarkable quadratic formula: \(n^2+n+41\)

It turns out that the formula will produce 40 primes for the consecutive integer values \(0 \leq n \leq 39\). However, when \(n=40\), \(40^2+40+41=40(40+1)+41\) is divisible by 41, and certainly when \(n=41\), \(41^2+41+41\) is clearly divisible by 41.

The incredible formula \(n^2-79n+1601\) produces 80 primes for the consecutive values \(0 \leq n \leq 79\). The product of the coefficients, -79 and 1601, is -126479.

Considering quadratics of the form:

\(n^2 + an + b\), where \(|a| < 1000\) and \(|b| \leq 1000\).

\(|n|\) is the absolute value of \(n\), e.g. \(|11|=11\) and \(|-4|=4\).

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\).

The only way to solve this problem is through brute force and search for all possible combinations for \(a\) and \(b\). Without reducing the search space, that would mean testing four million options. We need to reduce the solution space to achieve a manageable running speed.

Because the outcome of the equation must be prime for \(n = 0\), \(b\) also has to be prime. We can use the prime sieve from Euler Problem 3 to use only prime numbers from 2 to 1000, which reduces it from 2000 to 168 permutations.

When we insert \(n = 1\), it follows that \(a\) has to be an even number. If \(1+a+b\) has to be prime and \(b\) is a prime number, then \(a\) can only be an odd number. However, when \(b = 2\), \(a\), also has to be even.

We can use these deliberations to write some code to solve this problem. We cycle through \(a\) and \(b\) and increment a counter every time we find a prime number. The is.prime() function from Project Euler 7 helps to identify whether a number is prime.

ProjectEuler

  ## Project Euler 27: Quadratic Primes
  source("problem007.R", verbose = FALSE)

  seq.a <- seq(-999, 1001, 2) # a has to be odd
  seq.b <- (esieve(1000)) # b has to be prime
  max.n <- 0

  for (a in seq.a) {
    for (b in seq.b) {
      n <- 0 # Find sequence of primes for a and b
      while (is.prime(n^2 + a * n + b)) {
        n <- n + 1
      }
      ## Store maximum values
      if (n > max.n) {
        max.n <- n
        max.a <- a
        max.b <- b
      }
    }
  }

  answer <- max.a * max.b
  print(answer)
  
  ## Test the answer
  n <- 0:(max.count - 1)
  sapply(n^2 + max.a * n + max.b, is.prime)

As seen on R Bloggers