Project Euler #27: Quadratic Primes

PROBLEM STATEMENT

To explore a more difficult (and therefore more interesting problem), we will jump ahead to the 27th problem on the site. Unlike our previous problems, the description for this one is a bit more involved:

“Euler discovered the remarkable quadratic formula:

n^2 + n + 41
 
It turns out that the formula will produce 40 primes for the consecutive values n = 0 to 39. However, when n = 40, 402 + 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 was discovered, which produces 80 primes for the consecutive values n = 0 to 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| < 1000 and |n| is the modulus/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.”

SOLUTION STRATEGY

To start, we can observe that if we allow for all the values of a and b then we have an enormous sample space that will likely not compile and solve in less than 1 minute. As such, we should use some mathematical insight to narrow down the problem so that we can find a solution quickly. Let f(n) = n^2 + an + b. In order to have the longest chain of primes it must be the case that, at the least, f(0) should be prime. Thus, we have that f(0) = 0^2 + a*0 + b = b needs to be prime. This cuts down on the number of values of b significantly!

At the same time, we expect that f(1) = 1^2 + a*n + b = 1 + a + b should also be prime. We note that since the prior constraint made it so that b should be prime it follows that (b+1) is either 3 or it is an even number larger than 3. This is because the only even prime is 2, while every other prime is odd. In the first case, in order for (a+b+1) to be odd, a must be even. In the second case, a must be an odd number. Note that these two conditions together imply that a \equiv b(\mod2).

We have imposed some reasonable restrictions on our two values. Now all we need is an algorithm that will find the values for a and b that will yield the largest number of consecutive primes. Our algorithm will be:

  1. Generate the set of primes less than or equal to 1000. This will be the set that we define b from.
  2. For a given fixed value of b, consider all values of a such that -1000 \leq a \leq 1000 and a \equiv b(\mod 2). For the first value of a in this set, test the function where n=2. If the value is prime, then increment the number of found primes by 1 and increment to the next value of n for the given parameters. If the value is not prime, move on to the next set of parameters. That is, if there are more values of a to be tested then do so, otherwise move on to the next value of b.

R IMPLEMENTATION

Time to implement this solution in R! The first thing we need is a function to represent the function value for a given set of parameters. This is easily written:

f <- function(a,b,n) {return(n^2+a*n+b)}

While we can leverage our Sieve function that was mentioned in the post on Primes and Sieves to generate the primes from 1 to 1000, it could be the case that the function evaluates to some prime that is larger than that range. As such, we will want to have a separate function to test whether a given number is prime or not.

IsPrime <- function (N)
{
	if (N < 0) {return(FALSE)}
	if (N == 1)
		return (FALSE)
	
	if (N == 2 || N == 3)
		return (TRUE)

	if (N %% 2 == 0 || N %% 3 == 0)
		return (FALSE)

	Start = 5
	Stop = ceiling(sqrt(N))

	while (Start <= Stop)
	{
		if (N %% Start == 0)
			return (FALSE)

		Start = Start + 2
	}

	return (TRUE)
}

Finally, we need a block of code to manipulate these objects:

ptm <- proc.time()
primes = SieveOfErat(1000)
largestNumPrimes = 2
a_final = 0
b_final = 0

for(i in 1:length(primes))
{
	b = primes[i]
	for(j in (-1000):(1000))
	{
		n = 2
		currentNumPrimes = 2
		if((j + b) %% 2 != 0) {next}
		a = j
		while(IsPrime(f(a,b,n)))
		{
			n = n + 1
			currentNumPrimes = currentNumPrimes + 1
		}
		if (largestNumPrimes < currentNumPrimes)
		{
			a_final = a
			b_final = b
			largestNumPrimes = currentNumPrimes
		}
	}	
}
print(a_final*b_final)
proc.time() - ptm

The result we get is the correct one, -59,231 (obtained with a = -61 and b = 971) and this solution gets it in 10.34 seconds. I am sure there could be room for improvement here on the implementation to make this run faster; consider it a challenge for the interested reader! 🙂

Advertisements

Project Euler #5: Smallest Multiple

Problem 5 deals with the concept of a least common multiple, often abbreviated as LCM. The LCM for a set of numbers can be found conceptually by listing all of the factors of each individual number and then taking the smallest number that is on all such lists. For instance, if I were to ask what the LCM of 3 and 7 is, I could start by considering all of the multiples of each as follows:

  • Multiples of 3 include 3,6,9,12,15,18,21,24,27
  • Multiples of 7 include 7,14,21,28,35,42,49,56

The smallest number that is on both lists is 21, so LCM(3,7) = 21. Our problem, however, is a bit more complicated. It reads:

\fbox{\parbox{\textwidth}{2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder. What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?}}

 
To solve this problem, we note that by the Fundamental Theorem of Arithmetic every natural number has a unique factorization into a product of primes. Another way we can compute the LCM for a set of numbers N_1,N_2, \dots N_k is by the following algorithm:

  1. Determine the unique factorization for each number in the list in terms of primes. That is, for the kth number in our set, let N = p_1^{i_1}*p_2^{i_2} \cdots p_m^{i_m}
  2. , where each of the p‘s is a prime factor.

  3. Once we have this, we go through all the primes one at a time and determine what the largest exponent is in the factorizations of our numbers. That is, for the prime p_i suppose the exponents are e_1,e_2, \dots e_k. Then set e = \max \lbrace e_1,e_2, \dots e_k \rbrace
  4. as the exponent for that prime in our final result.

  5. The LCM will simply be the product of each prime raised to its corresponding exponent calculated in (2).

For the numbers 1 to 20, this is actually quite easy and we will not need to use a computer program to find the result! Consider that the primes that are less than or equal to 20 are 2,3,5,7,11,13,17, and 19. As a result, we know that the minimum exponent on each of these numbers is 1. On the other hand, some numbers less than 20 will have more than 1 multiple of some of these primes. Since 5^2 = 25 > 20, we know that no number less than 20 can have more than 1 factor of 5. This same reasoning applies to all the primes larger than 5, so we can keep all such primes with exponents of 1. For the primes smaller than 5 (2 and 3), we can see that the largest power of 2 is 4 (since 2^4 = 16 and 2^5 = 32) and the largest power of 3 is 2 (since 3^2 = 9 and 3^3 = 27). Thus, our result can be calculated as follows:

2^4*3^2*5*7*11*13*17*19 = 232,792,560
 
It is rare when we can work out one of these without the use of a computer program!

Project Euler #4: Largest Palindrome Product

The fourth problem reads:

\fbox{\parbox{\textwidth}{A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91*99. Find the largest palindrome made from the product of two three digit numbers.}}

 

Our strategy will be to start at the largest number that is a product of two three digit numbers and then determine whether or not it satisfies our requirements. To start, we need a function that will determine whether or not a given number is palindromic. To do this, we can convert the number to a string and compare the characters at corresponding locations to ensure that it is a palindrome. Let N = length of the string. Then a number is a palindrome if the ith character is equal to the N-i+1 character (note that R is a language where the string indices begin at 1 not 0). We can get away with looping from the first character to the middle character; this last character is given by \lfloor\frac{N}{2}\rfloor. The function in R is:

isPalindrome <- function(N)
{
	x = toString(N)
	upper = floor(nchar(x)/2)
	for(i in 1:upper)
	{
		if(substr(x,i,i) != substr(x,nchar(x)-i+1,nchar(x)-i+1)) {return(FALSE)}
	}
	return(TRUE)
}

Another thing that we need to check is that the number is divisible into a product of two 3 digit numbers. To do this, we can write the following function:

is3DigitDivisible <- function(N)
{
	if(N < 100^2 || N > 999^2) {return(FALSE)}
	for(i in 100:999)
	{
		if(N %% i == 0)
		{
			if(nchar(toString(N / i)) == 3)
			{
				return(TRUE)
			}
		}
	}
	return(FALSE)
}

These two functions alone will find the solution if we write some code to use them. The following code uses these functions to locate the largest palindrome that can be factored into a product of 2 three digit numbers.

ptm <- proc.time()
found = 0
for(i in 998001:10000)
{
	if(isPalindrome(i))
	{
		if(is3DigitDivisible(i))
		{
			found = i
			break
		}
	}
}
found
proc.time() - ptm

The result is 906609 and this was found in 1.84 seconds. Let us take a moment to note a few aspects of this solution. First, we can see that the loop starts with the largest number that is a product of 3 digit numbers; this is 999*999 = 998,001. Further, the loop ends at the smallest product of three digit numbers ; this is 100*100 = 10,000. When we reach a number, we first check that this is a palindrome. Supposing it is, we then check that it is divisible by 3. Note that if we reverse the order that this is done (i.e. first check that it is divisible by 3 and then check that it is a palindrome) we get a run time of 39.23 seconds. This is much larger than our previous time and this is of course because the set of numbers that are divisible into two three digit numbers is much larger than the set of palindromic numbers, so we eliminate more unnecessary work by first checking that a number is palindromic.

Now one can further improve this solution by noting that mathematically, if we assume that the number of digits of the largest palindrome is 6 digits, then it must be that it is written ABCCBA. We can represent our palindrome as follows:

N = 100000A + 10000B + 1000C + 100C + 10B + 1A
 
Then collecting like terms, we have:
 
N = 100001A + 10010B + 1100C = 11(9091A + 910B + 100C)
 
What we have just showed is that the palindrome itself must be divisible by ll, regardless of the values of its digits. Thus, we can change our solution to use this fact by starting at the largest product of 3 digit numbers that is also divisible by 11; this is 997997, and we can proceed downwards by multiples of 11.

ptm <- proc.time()
found = 0
for(i in seq(from=997997,to=10000, by =-11))
{
	if(isPalindrome(i))
	{
		if(is3DigitDivisible(i))
		{
			found = i
			break
		}
	}
}
found
proc.time() - ptm

When we do this, we achieve the same result, but it is much faster at 0.21 seconds.