Project Euler 27: Find a quadratic formula that produces the maximum number of primes for consecutive values of n

We are back to saying something about primes in Problem 27 of Project Euler. So you may already now want to go back to the solution to Problem 10 for an efficient implementation of Sieve of Eratosthenes since we are likely to need that. The problem reads

Euler published the remarkable quadratic formula:

n² + 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² + 41 + 41 is clearly divisible by 41.

Using computers, the incredible formula  n² – 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² + an + b, where |a| < 1000 and |b| < 1000

where |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.

For me this one reeks of brute force, since it is obvious that we can run through all possible values of a and b. Yes that is 4.000.000 solutions we need to check, but it should be possible.

Brute Force

Before presenting you with the main algorithm, I need a small utility function to check if a number is a prime. In this case I have assumed that I have an ordered list of primes larger than the number I want to test, and thus we can write a small function which looks like this.

private bool isPrime(int testNumber) {
int i = 0;
while (primes[i] <= testNumber) {
if (primes[i] == testNumber) {
return true;
}
i++;
}
return false;
}

Now you might ask, why not use the function API list function Contains. I tried it, but the performance completely degraded. So by writing my own function I won a factor 50-100 in execution speed.   With that in place, it is pretty easy to write a brute force function to test all possibilities.

int aMax = 0, bMax = 0, nMax = 0;
int[] primes = ESieve(87400);

for (int a = -1000; a <= 1000; a++) {
for(int b = -1000; b <= 1000; b++){
int n = 0;
while (isPrime(Math.Abs(n * n + a * n + b))) {
n++;
}
if (n > nMax) {
aMax = a;
bMax = b;
nMax = n;
}
}
}

It runs fairly fast, and I get the result

A sequence of length 71, is generated by a=-61, b=971, the product is -59231
Solution took 660 ms

We are well within the one minute limit that the project sets. However, I want to give it another spin, to see if we can decrease the solution space a bit. However, I want to adhere to the brute force method.

Shrinking solution space

The quadratic formula has to provide us with primes all the way from n= 0. What that means is that at n = 0

n2 + an + b  = b

And thus b has to be a prime. That limits b from 2001 possibilities to 336. Or a factor 6 reduction in size.

For n = 1 we have that

n2 + an + b  = 1 + a + b

We know that all primes except for 2 are odd. That means if 1 + a + b has to be odd, a has to be odd as well so the formula to provide a prime for n = 1.  If b = 2, then a has to be even.

That cuts the possibilities of a by 50%. So we have limited the total solution space by around a factor 12.

int aMax = 0, bMax = 0, nMax = 0;
primes = ESieve(87400);
int[] bPos = ESieve(1000);

for (int a = -999; a < 1001; a +=2) {
for (int i = 1; i < bPos.Length; i++) {
for (int j = 0; j < 2; j++) {
int n = 0;
int sign = (j == 0) ? 1 : -1; int aodd = (i % 2 == 0) ? -1 : 0; // Making a even if b is even
while (isPrime(Math.Abs(n * n + (a + aodd) * n + sign*bPos[i]))) {
n++;
}
if (n > nMax) {
aMax = a;
bMax = bPos[i];
nMax = n;
}
}
}
}

The execution of this gives

A sequence of length 71, is generated by a=-61, b=971, the product is -59231
Solution took 115 ms

Which is a factor 6 faster than the original solution.

Wrapping up

There isn’t too much in this problem I think. I have found no curious formulas or anything. Except for the few insights I have provided here, I have found no short cuts to the solution, so if you see something I didn’t don’t hesitate to let me know.

As usual you can find the source code  here to play around with. Posted by Kristian Jerry Iuliano

.75*x^2 + 1.5 * x + 23 = rich in primes and products of two primes when x = even number Kristian

Hi Jerry

Thanks for that quadratic formula. Though it doesn’t really fit the problem description it is a rather interesting set of parameters for a quadratic formula that does indeed produce lots of primes. I quickly tested it up to 20 where all the even numbers result in primes. Pretty interesting.

/Kristian Ronnie

I used both a and be to be prime.. and got passed.. Since I know any number can be a sum of primes and one is prime, there is much possibility that n*(n+a) to be prime as well.. (I have done much work, you can check it on yourself)

But everytime I run it in different ranges, answer produced is correct (same in both slow and fast brute force)

my code goes follows in python

from sympy import isprime

primes=
for a in range(2,1000):
if isprime(a):
primes.append(a)

b=[]
for a in primes:
b+= [-1*a]

for a in b:
primes.append(a)

def fn(a, b, n):
return n**2+a*n+b

ns={}
aa=[]
bb=[]
for a in primes:
for b in primes:
#if a!=b:
n=1
while n:
pr=fn(a,b,n)
if isprime(pr):
aa.append(n)
ns[n]=a*b
#print a,b,n,pr
n+=1

else:
break

a=max(aa)
print ns[a]

It gave me results in 405ms in a slow language like python.. Kristian

Interesting observation. In general I try to avoid using unjustified observations. However, I do like the observation and I wonder if it could be proven that a needs to be a prime too. Ronnie

I’m trying to prove the same too.. but there is nothing like if m=2n, m is even

Theres nothing like that for primes… And sometimes I feel good if I shoot blindly and it hits 🙂 jerome

smartness: b is prime and a is odd !
(more in the spirit of euler) Son-Huy TRAN

I found that if 1 + a + b is prime, so :
1 + a + b > 1
=> a + b > 0

Because b is prime so b > 0, so if a < 0 : |a| < |b|

So we can iterate like this : for each value of b, we iterate the odd numbers from -b to 1000. Robin

I believe this is a pretty interesting problem, and for larger limits on $|a|$ and $|b|$ I’ve noticed this is really time consuming. I’ve done tests for both letting $a\in\mathbb{R}: |a|<1000$ and for $a$ only being odd. However with my code I got two very quick times (in Java)

[java]
public class ProjectEuler27 {
public static boolean isPrime(int n){
if(n<2 || (n>2 && n%2==0)){
return false;
}else if(n==2){
return true;
}else{
for(int i=3; i<=Math.sqrt(n); i=i+2){
if(n%i==0){
return false;
}
}
return true;
}
}
public static void main(String[] args){
int tmp=0; int A = 0; int B = 0; int limit = 1000;
int nbr = 0;
long start = System.currentTimeMillis();
for(int a=-limit+1; a<limit; a=a+2){ //2 if only odd nbrs
for(int b=1; b<limit; b=b+2){
while(isPrime(nbr*nbr+a*nbr+b)){ //n^2+an+b
nbr++;
}
if(nbr>tmp){ //store the values for later out-print
tmp = nbr; A=a; B=b;
}
nbr=0;
}
}
long time = System.currentTimeMillis() - start;
System.out.println("The product a*b is: " + A*B + "\nwith a=" + A +
" and b=" + B + " with limit " + limit +".");
System.out.println("Formula generated " + tmp + " primes.");
System.out.println("Time: " + time + " ms.");
}
}
[/java]

If $a$ only can be odd, I got

The product a*b is: -59231
with a=-61 and b=971 with limit 1000.
Formula generated 71 primes.
Time: 36 ms.

And if we let $a$ increment through all the numbers

The product a*b is: -59231
with A=-61 and B=971 with limit 1000.
Formula generated 71 primes.
Time: 56 ms.

I haven't been able to prove that a needs to be odd to produce lots of primes, but that assumption have worked for my tested values of $|a|<n$ where $n$ denotes the limit. Kristian

Thanks for the comment. I edited it a bit to make the code tags work again. I have not edited the content though. As you can see in my post a indeed has to be odd, and that is provable, so I am glad your tests show the same 🙂 kpwbo

So :
– b has to be prime (and positive);

when n = 0
0^2 + 0*a + b = prime
b = prime

– a has to be negative, with an absolute value smaller than b;

when n = 1
1^2 + 1*a + b = prime
1 + a + b = prime
as a prime is always bigger than 1 :
1 + a + b > 1
a + b > 0
as b is positive, a must be negative
to make the inequality true, abs(a) must be smaller than b

– a has to be odd

when n = 1
1^2 + 1*a + b = prime
1 + a + b = prime
as a prime is almost always odd
1 + a + odd = odd
a + even = odd
a must be odd to make this true

I based my code on this. It is in Ruby (I’m a beginner), but it should be readable.

require "mathn"

largest,max_a,max_b = 0,0,0

(2..999).each { |b|
if b.prime? == false
next
end
(-b..-1).each { |a|
if a.even?
next
end
count = 0
n = 0
while true
if (n**2 + a*n + b).prime?
count += 1
n += 1
else
break
end
end
if count > largest
largest,max_a,max_b = count,a,b
end
}
}

puts max_a*max_b Chris Goldfrap

I’m a novice and that may be why I don’t see how n^2 – 79n + 1601 is any more exciting than Euler’s original n^2 + n + 41. I’ve input values for n for your formula and get 1523 for both n = 1 and n = 78. It looks like a parabola (?) in which the same primes are repeated on the up curve as on the down. There are no extra ones to Euler’s. What’s remarkable about the latter is that you can hop along starting at 41 lengthening your stride by 2 each time and land on a prime each time 39 times. Has that been surpassed? I’ve been looking at what happens if you continue with Euler’s formula with n = 40 and beyond. There are some interesting prime runs, and computing might eventually reveal some very long ones. Kristian

I must admit that I haven’t really given much thought to investigate the solution further. But it seems to hold quite a few interesting gifts when you look closer. Kartik

I am curious as to how you arrived at the number 87400 as an upper bound for the primes. frgsd

You have to have in mind that the a < |1000| means that 'a' has to be between [-1001, 999]. If the params from the solution were the limits, it would no be right. Kristian

Nope, it means that a is in the interval [-999, 999]. Though you are right, it should be strictly less than and not less than or equal. Kristian

@Kartik I am sure I had a reason for it, but today I have no clue why . Rosendo Salinas

Hi! I found in Wolfram MathWorld that n² – 79n + 1601 = (n-40)² + (n-40) + 41. I calculated (with this formula):
x²-x+41<1000, x<32
And:
(n-31)² + (n-31) + 41 = n² – 61n + 971 🙂
The problem is that I don't know how demonstrate that this sequence generates the maximum number of primes. Kristian

As you can see from the blog post, my best attempt is to brute force the search space. I haven’t been able to find a way to prove it analytically. Darren

Hi Kristian,

I think there is another opportunity to reduce the search space of the solution listed in the post. It can be shown that given a and b, the maximum number of primes for consecutive values of n cannot exceed b (since b^2+ab+b is divisible by b, thus not prime). Therefore, we could have iterated b in the outer loop, from the largest prime that is smaller than 1000 to 3, and tested possible a in the inner loop. Once b is smaller than the current best (nMax in your code), we exit from the loop (later b’s cannot do better than that) and get the answer. This optimization would save many possible (useless) combinations of a and b and end the search quicker.

PS. This site is great. Every time I tackle a problem in Project Euler, I come here to see if there is any better solution. Kristian

Yes of course, that is a brilliant insight. Shaun Mathew

i got for a = -999, and b = 61, 1011 consecutive primes (which is greater than 71, hence this the max product). I am not sure what I am doing wrong here.

using a simple piece of code as follows in python,

count = 0
i = 0
while is_prime(i**2 – 999*i + 61):
count += 1
i += 1

return count

I seem to be very confused by this, unless my algorithm for is_prime is incorrect?
Any feedback provided would be much appreciated Kristian

I think your is_prime is wrong. Perhaps because you test negative numbers. For i=11 you get 10807, which according to http://www.mathblog.dk/tools/prime-factorization/

10807 has the following properties:
It is composite and has the factorization: 101×107
It has 2 factors, and their sum is 208.
It has 2 distinct factors, and their sum is 208.
It has 4 divisors, and their sum is 11016.
There are 10600 numbers, less than or equal to 10807, that are relatively prime to it. In brute force output line 1, you said length 71 but actually length 72 because we get prime for n=0 as well. Martin Pahlplatz

You set aodd=-1 for every even index i of bPos, in stead of only for bPos==2 !? andrew

Shaun,

your is_prime is correct. The values you’ve posted are the ones you get if you don’t plug in the absolute value of f(x), the quadratic output. Alternatively, apply abs to the parameter inside is_prime.

At least for my version of my_prime it works, but if I omit the abs part my program outputs the same numbers as yours, -999, 67 (with 1011 consecutive primes). Sinan

This solved in 10ms and when the ESieve is excluded, the solution loop takes only 1ms.

Stopwatch sw = new Stopwatch();
sw.Restart();
int[] primes = ESieve(1000000);
bool[] isPrime = listifier(primes);
int bestCount = 0;
int bestPrime = 0;
int bestJ = 0;

for (int i = 1; primes[i] < 1000; i++) {
int currentPrime = primes[i];//b has to be a prime number
for (int j = - 1 * (int) Math.Sqrt(4*currentPrime); j < 999; j++) {
//(j = a) a^2 bestCount) {
bestCount = primesGenerated;
bestPrime = currentPrime;
bestJ = j;
}
}

}

sw.Stop();
Console.Write("result\t{0}\t{1}\n", bestJ, bestPrime);
Console.Write("time:\t{0}\n\n", ((double)(sw.Elapsed.TotalMilliseconds)).ToString("0.000 000 ms")); Sinan
Stopwatch sw = new Stopwatch();
sw.Restart();
int[] primes = ESieve(1000000);
bool[] isPrime = listifier(primes);
int bestCount = 0;
int bestPrime = 0;
int bestJ = 0;

for (int i = 1; primes[i] &lt; 1000; i++) {
int currentPrime = primes[i];//b has to be a prime number
for (int j = - 1 * (int) Math.Sqrt(4*currentPrime); j &lt; 999; j++) {
//(j = a) a^2&lt;4b because otherwise there will be a non-negative solution
//non-negative solution means that the output of the polinomial
//will take negative values
int primesGenerated = 1;
while (true) {
if (!isPrime[primesGenerated * primesGenerated + j * primesGenerated + currentPrime])
//I created an array which would tell me if a number is prime
break;
primesGenerated++;
}
if (primesGenerated &gt; bestCount) {
bestCount = primesGenerated;
bestPrime = currentPrime;
bestJ = j;
}
}

}

sw.Stop();
Console.Write(&quot;result\t{0}\t{1}\n&quot;, bestJ, bestPrime);
Console.Write(&quot;time:\t{0}\n\n&quot;, ((double)(sw.Elapsed.TotalMilliseconds)).ToString(&quot;0.000 000 ms&quot;)); Mihnea

hey there, I’m getting the right result but my run time is insanely slow (203s), and I can’t seem to find where the egregious inefficiency is. Im using python, does somebody mind taking a look at it:

def createPrimes(n):
global primes
numbers_for_primes=range(0,n+1)
# identify all prime numbers
for prime in numbers_for_primes:
if(prime&lt;2): #we will set numbers to 0 so this breaks it. we put this first since most important
continue
elif(prime &gt; n**0.5): # we end loop when greater than sqrt(n) because it wouldve marked prime before that
break
else:
for i in range(prime**2,n,prime): #mark all numbers with 0
numbers_for_primes[i]=0
numbers_for_primes=0

#build a new array of primes based on whats not set to 1
primes = filter(lambda x: x != 0, numbers_for_primes)

def is_prime(number):
for i in range(0,len(primes)):
if(primes[i]&lt;number):
continue
elif(number==primes[i]):
return True
else: #primes now bigger than number so no need
return False

def main():
# find the first 78.5k primes in 1M numbers or 9.5K given 100K
global primes
createPrimes(87400)
max_a = 999 #maximum a coefficient
max_b = 999 #maximum b coefficient
# when n=0, b has to be prime, so b must be in list of primes less than 1K
pos_b_list=filter(lambda x: x&lt;=1000, primes)
neg_b_list=[-1*i for i in pos_b_list]
b_list=pos_b_list + neg_b_list
coefficient_product = 0 # we want a*b to be the number we return
max_primes_generated = 0 # we need to keep this meta data to wait for product
# we will need to look at negative coefficients too
for a in range(-1*max_a,max_a+1):
#print(&quot;A is&quot;, a)
for b in b_list:
n=0
while(is_prime(abs(n**2+a*n+b))):
n+=1
#numPrimes=consecutive_primes(a,b)
if(n&gt;max_primes_generated):
max_primes_generated=n
coefficient_product=a*b
print(&quot;a, b, number of primes is:&quot;, a,b,n)

print max_primes_generated, coefficient_product

if __name__ == &quot;__main__&quot;:
main() John Jacobs

hey guys

my algorithm was able to get the correct answer but it took a crazy long amount of time (200 s). I’ve looked it over time and time again and even compared it to other solutions (and just simple logic of the fact that there are only, at most, 4M numbers to look through) and I just can’t seem to find what I’m doing wrong. Does somebody mind taking a look at it for me. I’d greatly appreciate it.

# Uses the sieve of erostenes to find all primes up to a certain n
# Returns an array of prime numbers
def createPrimes(n):
global primes
numbers_for_primes=range(0,n+1)
# identify all prime numbers
for prime in numbers_for_primes:
if(prime&lt;2): #we will set numbers to 0 so this breaks it. we put this first since most important
continue
elif(prime &gt; n**0.5): # we end loop when greater than sqrt(n) because it wouldve marked prime before that
break
else:
for i in range(prime**2,n,prime): #mark all numbers with 0
numbers_for_primes[i]=0
numbers_for_primes=0

#build a new array of primes based on whats not set to 1
primes = filter(lambda x: x != 0, numbers_for_primes)

def is_prime(number):
for i in range(0,len(primes)):
if(primes[i]&lt;number):
continue
elif(number==primes[i]):
return True
else: #primes now bigger than number so no need
return False

def main():
# find the first 78.5k primes in 1M numbers or 9.5K given 100K
global primes
createPrimes(87400)
max_a = 999 #maximum a coefficient
max_b = 999 #maximum b coefficient
# when n=0, b has to be prime, so b must be in list of primes less than 1K
pos_b_list=filter(lambda x: x&lt;=1000, primes)
neg_b_list=[-1*i for i in pos_b_list]
b_list=pos_b_list + neg_b_list
coefficient_product = 0 # we want a*b to be the number we return
max_primes_generated = 0 # we need to keep this meta data to wait for product
# we will need to look at negative coefficients too
for a in range(-1*max_a,max_a+1):
#print(&quot;A is&quot;, a)
for b in b_list:
n=0
while(is_prime(abs(n**2+a*n+b))):
n+=1
if(n&gt;max_primes_generated):
max_primes_generated=n
coefficient_product=a*b
#print(&quot;a, b, number of primes is:&quot;, a,b,n)

print max_primes_generated, coefficient_product

if __name__ == &quot;__main__&quot;:
main() Anon

Have you tried solving this using the heegner number? the discriminant must be one of the heegner numbers to be able to generate primes. i.e. a^2-4b should be one of 1,2,3,7,11,19,43,67, and 163. Also the upper bound for the loop can be set to the heegner number Jean-Marie Hachey

Quadratic equations giving primes.
Prime numbers being only positive, it is evident that only eqn  by Euler will produce specifically sequences of primes without duplication.
On the other hand, eqn  and  will generate two sequences of duplicated primes on the descending and ascending branch of the parabola.

Fig. 1, Fig. 2 and Fig. 3 illustrate the results. Julian

Solution:
Duration in ms: 29

maxN: 71
maxA: -61
maxB: 971
maxA * maxB = -59231
Code:

import java.util.BitSet;
import java.util.HashSet;

public static void main(String[] args){
int n = 0;
int maxN = 0;
int maxA = 0;
int maxB = 0;

long startTime = System.currentTimeMillis();
HashSet primeList = ESieve(10000);

for(int a = -999; a < 1000; a+=2){

for(int b = 2; b maxN){
maxN = n;
maxA = a;
maxB = b;
}
}
}

System.out.println(“Duration in ms: ” + (System.currentTimeMillis()-startTime));

System.out.println(“\nmaxN: ” + maxN);
System.out.println(“maxA: ” + maxA);
System.out.println(“maxB: ” + maxB);
System.out.println(“maxA * maxB = ” + (maxA*maxB));

}

public static HashSet ESieve(int limit){
int border = (((int)Math.sqrt(limit))-1)/2;
int halfSize = (limit-1)/2;
BitSet valStore = new BitSet();
valStore.set(0, halfSize+1, true);

int tmpPrime = 0;
for(int i = 1; i <= border; i++){
tmpPrime = i*2+1;

if(valStore.get(i)){
for(int c = (tmpPrime*tmpPrime-1)/2; c <= halfSize; c += tmpPrime){
valStore.set(c, false);
}
}
}

HashSet primes = new HashSet();

for(int k = 1; k < halfSize; k++){
if(valStore.get(k)){
}
}

return primes;
}
} Shruti

In the line where aodd is assigned:
int aodd = (i % 2 == 0) ? -1 : 0;

Shouldn’t it have been:
int aodd = (bPos[i] % 2 == 0) ? -1 : 0;

[…] this problem is through brute force and reduce the solution space to optimise it for speed (source: mathblog.dk). Because the outcome of the equation must be prime for , also has to be prime. We can use the […]