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

Written by on 12 February 2011

Topics: Project Euler

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.

33 Comments For This Post I'd Love to Hear Yours!

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

  2. Kristian says:

    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

  3. Ronnie says:

    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=[1]
    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..

  4. Kristian says:

    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.

  5. Ronnie says:

    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 🙂

  6. jerome says:

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

  7. Son-Huy TRAN says:

    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.

  8. Robin says:

    I believe this is a pretty interesting problem, and for larger limits on and I’ve noticed this is really time consuming. I’ve done tests for both letting and for only being odd. However with my code I got two very quick times (in Java)

    If 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 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 where denotes the limit.

  9. Kristian says:

    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 🙂

  10. kpwbo says:

    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
  11. Chris Goldfrap says:

    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.

  12. Kristian says:

    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.

  13. Kartik says:

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

  14. frgsd says:

    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.

  15. Kristian says:

    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.

  16. Kristian says:

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

  17. Rosendo Salinas says:

    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.

    Link: http://mathworld.wolfram.com/Prime-GeneratingPolynomial.html

  18. Kristian says:

    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.

  19. Darren says:

    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.

  20. Kristian says:

    Yes of course, that is a brilliant insight.

  21. Shaun Mathew says:

    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

  22. Kristian says:

    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.

  23. Madhur Panwar says:

    In brute force output line 1, you said length 71 but actually length 72 because we get prime for n=0 as well.

  24. Martin Pahlplatz says:

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

  25. andrew says:

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

  26. Sinan says:

    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"));

  27. Sinan says:
                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;));
  28. Mihnea says:

    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[1]=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()
    
    
  29. 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[1]=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()
  30. Anon says:

    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

  31. Jean-Marie Hachey says:

    Quadratic equations giving primes.
    Prime numbers being only positive, it is evident that only eqn [1] by Euler will produce specifically sequences of primes without duplication.
    On the other hand, eqn [2] and [3] 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.

    http://img11.hostingpics.net/pics/610000pe27fig123quadf.jpg

  32. Julian says:

    Solution:
    Duration in ms: 29

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

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

    public class QuadraticPrimes{
    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();
    primes.add(2);

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

    return primes;
    }
    }

  33. Shruti says:

    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;

Leave a Comment Here's Your Chance to Be Heard!

You can use these html tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

You can use short tags like: [code language="language"][/code]

The code tag supports the following languages: as3, bash, coldfusion, c, cpp, csharp, css, delphi, diff,erlang, groovy, javascript, java, javafx, perl, php, plain, powershell, python, ruby, scala, sql, tex, vb, xml

You can use [latex]your formula[/latex] if you want to format something using latex

This site uses cookies.