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| = 4Find 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.
.75*x^2 + 1.5 * x + 23 = rich in primes and products of two primes when x = even number
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
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
It gave me results in 405ms in a slow language like python..
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.
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 🙂
smartness: b is prime and a is odd !
(more in the spirit of euler)
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 believe this is a pretty interesting problem, and for larger limits on [latex]|a|[/latex] and [latex]|b|[/latex] I’ve noticed this is really time consuming. I’ve done tests for both letting [latex]a\in\mathbb{R}: |a|<1000[/latex] and for [latex]a[/latex] 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 [latex]a[/latex] 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 [latex]a[/latex] 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 [latex]|a|<n[/latex] where [latex]n[/latex] denotes the limit.
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 🙂
So :
– b has to be prime (and positive);
– a has to be negative, with an absolute value smaller than b;
– a has to be odd
I based my code on this. It is in Ruby (I’m a beginner), but it should be readable.
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 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.
I am curious as to how you arrived at the number 87400 as an upper bound for the primes.
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.
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.
@Kartik I am sure I had a reason for it, but today I have no clue why .
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
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.
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.
Yes of course, that is a brilliant insight.
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
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.
You set aodd=-1 for every even index i of bPos, in stead of only for bPos==2 !?
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).
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"));
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:
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.
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
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
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;
}
}
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 […]
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
int[] primes = ESieve(87400);
why did you take 87400 value here…?
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?