In Problem 146 of Project Euler we are working with primes again, and some quite big ones even. The problem reads

The smallest positive integernfor which the numbersn^{2}+1,n^{2}+3,n^{2}+7,n^{2}+9,n^{2}+13, andn^{2}+27 are consecutive primes is 10. The sum of all such integersnbelow one-million is 1242490.

What is the sum of all such integersnbelow 150 million?

At first I thought I could just make a sieve up to 150 million and then check if the numbers were contained in that. However, rereading the problem I realized I was completely wrong. So in a pure brute force solution we would need to check 150 million values of n and up to 13 numbers for each, since we both need to check that the given numbers are prime. But also that the odd numbers inbetween are not prime. So potentially we have to check 1950 million numbers for primality, which is a moderately expensive operation.

For the primality test we will use a Miller-Rabin (Which I implemented for Problem 58) test which can be made deterministic up to pretty large values. However, there is a still a risk that some of the numbers in the high end of the range that we are going to check here might be false positives.

However, before going heads on into this problem let us see if we can removed some of the numbers.

## Removing some numbers

I decided to optimize the first part by reducing the amount of numbers we need to check. To do that, we need to take a better look at the prime formulas, starting with the first, which is . In order for this to be a prime, it must not be divisible by any prime below it. In particular, if , then must not hold. Let’s see what happens if . Then , which is ok. However, if , then , but then is divisible by , and hence not a prime. So we can see that, in order for to be a prime, must hold.

Which means that we have now reduced the number of checks to half the original level.

In the above section we used the fact that

(a*b) % c = ((a % c) * (b % c) )% c

I know it is written in a slight different notation, since it is copy paste from Problem 48.

Doing the same analysis for 5 we can see that

So either we get a remainder of 0, 1 or 4. If we get a remainder of 1 then the number which means that this is not a series of primes, since the is not prime. If we get a remainder of 4 then cannot be prime. Which means that every number we check must be divisible by 5.

So now we know that every number must be divisible by 5 and 2. Which means they must be divisible by 10. So we have already reduced the numbers by a factor 10.

I have done the same thing for 3 and 7 and concluded that and or .

## Removing some of the checks for in-between numbers

With the knowledge we have now for the values of modulo different numbers, I started looking at the fact that we must check to see if are primes as well.

We know that which means that the +5,+15 and +25 must also be divisible by 5. Thus we don’t need to check them. Reducing the list to

We know that which means that the +11,+17 and +23 must be divisible by 3. Thus we don’t need to check them. Reducing the list to

## Coding a first try

With these piece of analysis I implemented the following piece of code

long limit = 150000000; long result = 0; for (long i = 10; i <= limit; i+=10) { long squared = i * i; if (squared % 3 != 1) continue; if (squared % 7 != 2 && squared % 7 != 3) continue; if (squared % 9 == 0 || squared % 13 == 0 || squared % 27 == 0) continue; if (IsProbablePrime(squared + 1) && IsProbablePrime(squared + 3) && IsProbablePrime(squared + 7) && IsProbablePrime(squared + 9) && IsProbablePrime(squared + 13) && IsProbablePrime(squared + 27) && !IsProbablePrime(squared + 19) && !IsProbablePrime(squared + 21)) result += i; }

This code will run in about 90 seconds on my machine. I did some tweaks to it which complicated the code and got it to around 70 seconds. However, those tweaks did do well enough. So I gave up and cried for help to Bjarki. Which suggested the next step.

## Automating the modulo check

Let us automate the modulo check since that is a far less expensive check to perform on the numbers.

The following function returns a boolean array for a given input number n where each entry *i* gives us false is for any .

public bool[] Mods(int m) { bool[] res = new bool[m]; for (int i = 0; i < m; i++) { res[i] = true; for (int j = 0; j < Add.Length; j++) { if ((i * i + Add[j]) % m == 0) { res[i] = false; break; } } } return res; }

So now we don’t have to perform this by hand and thus limit us to the numbers, 2,3,5 and 7. Now we can do this for all the primes we want.

Doing it for for the first 10 primes gives us

Prime & Valid remainders

2: 0

3: 1, 2

5: 0

7: 3, 4

11: 0, 1, 4, 5, 6, 7, 10

13: 1, 3, 4, 9, 10, 12

17: 0, 1, 3, 6, 7, 8, 9, 10, 11, 14, 16

19: 0, 1, 2, 3, 6, 8, 9, 10, 11, 13, 16, 17, 18

23: 0, 1, 2, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22

29: 0, 1, 2, 3, 5, 6, 8, 9, 10, 11, 13, 16, 18, 19, 20, 21, 23, 24, 26, 27, 28

Now we will go through each integer from 1 to 150 million, but only inspect further if it has a valid remainder for each of the primes in the table. For example, when , we can see that it has a valid remainder for each of the primes in the table. For the prime , we see that , but is a valid remainder for , and for the prime , we see that , which is a valid remainder for . The same can be seen for each of the other primes in the table, and thus we will go on to check if is a valid number for the problem. Another example, when , we can see that, for the prime , , but is not a valid remainder for , so we can safely skip and go on to the next integer.

This means that we can rewrite the main loop as

int limit = 150000000; long result = 0; int[] primes = Sieve(2, 5000); Tuple<int, bool[]>[] mods = new Tuple<int, bool[]>[primes.Length]; for (int i = 0; i < primes.Length; i++) { mods[i] = new Tuple<int, bool[]>(primes[i], Mods(primes[i])); } for (long i = 10; i < limit; i+=10) { bool ok = true; for (int j = 0; ok && j < primes.Length && i*i > primes[j]; j++) { ok = mods[j].Item2[i % mods[j].Item1]; } for (int j = 0; ok && j < Add.Length; j++) { ok = IsProbablePrime(i * i + Add[j]); } for (int j = 0; ok && j < NotAdd.Length; j++) { ok = !IsProbablePrime(i * i + NotAdd[j]); } if (ok) { result += i; } }

Where the first inner loop (line 13-15) checks if the current number falls for the modulo check. If that is not the case then we will throw the expensive primality test after it.

Testing a few different values it seems that all primes under 5000 gives the best speed. Which means around 700 primes. Something we could not have done by hand. Doing it this way means we can solve the problem in around 900ms.

## Wrapping it all up

For this problem there is quite a bit of code which I have not shown and which is necessary for the problem to run. However, you can find it all here. Both my first try and Bjarki’s automation of the mod checking.

Table 1

Sum of exponents n in the sequence of primes:

(n^2) + 1, (n^2) + 3, (n^2) + 7, (n^2) + 9, (n^2) + 13, (n^2) + 27.

http://img11.hostingpics.net/pics/974713pe146tab1.jpg

___

Sources:

1) Project Euler 146

Kristian’s algorithm

http://www.mathblog.dk/files/euler/Problem146.cs

2) Microsoft Visual C# 2010 Express

(Reference added: System.Numerics)

3) Microsoft Office Excel (2007)

4) MathBlog

Prime Factorization Calculator

http://www.mathblog.dk/tools/prime-factorization/

Precision…

In the above comment …

We should read « Sum of integers n » (Instead of « Sum of exponents n »)

Hi! just have a small question, how do you know that n^2 = 3 mod{7} ?

100 % 7 is 2, but it is one of the n.

I don’t conclude that it is a prime number, I am just saying that it is a necessary condition, and if that is not fulfilled then we don’t have to do further checks.

Prime numbers that end with “57” occur more often than any other 2-digit ending among the first one thousand primes.

https://primes.utm.edu/curios/page.php?curio_id=29564