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

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

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

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 🙂

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.

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.

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

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

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.

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

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

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

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

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<2): #we will set numbers to 0 so this breaks it. we put this first since most important
			continue
		elif(prime > 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]<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<=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("A is", a)
		for b in b_list:
			n=0
			while(is_prime(abs(n**2+a*n+b))):
				n+=1
			if(n>max_primes_generated):
				max_primes_generated=n
				coefficient_product=a*b
				#print("a, b, number of primes is:", a,b,n)

	print max_primes_generated, coefficient_product

if __name__ == "__main__":
    main()

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

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;

your advance solution is wrong……. at row 9

int sign = (j == 0) ? 1 : -1; int aodd = (i % 2 == 0) ? -1 : 0; // Making a even if b is even


it must be
int sign = (j == 0) ? 1 : -1; int aodd = (primes[i] % 2 == 0) ? -1 : 0; // Making a even if b is even

The code dealing with parity of a and b could be removed, since b is prime and it’s easy to show that for b = 2 (the only even prime) the equation n^2 + an + 2 is not a prime already at n = 2, so we can deal with b starting from 3 to 1000. and odd a’s only.
Regarding a – I think there is a better limitation you can impose on it. If we use the quadratic formula for ax^2 + bx + c then we know it has factors if and only if the part under the square is non-negative, so we can add the condition that a^2 – 4b < 0. This condition gives us that a is between -2sqrt(b) and 2sqrt(b).

Is that a valid argument to take?