Hi, I tried doing Project Euler #27 today, which asks the following:

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, 40^(2) + 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.

So, I did my thinking and came up with the following code in a few minutes:

#include <iostream>

using namespace std;

int* sieve(int i)
{
	bool* primes = new bool[i];
	int* primelist = new int[i];
	
	primes[0] = false;
	primes[1] = false;
	
	for (int j = 2; j < i; j++)
		primes[j] = true;
	
	for (int j = 2; j * j < i; j++)
	{
		if (primes[j])
		{
			for (int k = j; k*j < i; k++)
				primes[k*j] = false;
		}
	}
	
	for (int k = 2; k < i; k++)
	{
		if (primes[k])
			primelist[k] = k;
	}
	
	return primelist;
}

bool isPrime(int check)
{
	int* primes = sieve(1000);
	
	for (int i = 0; i < 1000; i++)
	{
		if (check == primes[i])
			return true;
	}
	
	return false;
}

int main()
{
	int* primes = sieve(1000);
	int max_co = 0;
	int max_co_sec = 0;
	int index = 2;
	int max = 0;
	int x;
	
	for (int a = -999; a < 1000; a++)
	{		
		for (int b = -999; b < 1000; b = primes[index++])
		{	
			x = 0;
			
			while (isPrime(x*x+a*x+b))
				x++;
				
			if (x > max)
			{
				max = x;
				max_co = a;
				max_co_sec = b;
			}
		}
	}
	
	cout << max_co * max_co_sec << endl;
	
	return 0;
}

What I did:

Function sieve() returns an array of prime numbers using the sieve of eratosthenes. This part is OK. I validated the function by generating random ranges of prime numbers and using brute-force to check them. It works.

Next, in the main(), there are two loops for checking combinations of n*n+a*n+b. isPrime() goes over the generated primes array to see if the number exists anywhere in it. If it does, then it is a prime.

This is basically what my program does.. Problem is, that I'm getting the wrong output. What I'm getting is something like -75913. The correct answer is something like -51318.

What am I doing wrong here? And, yes, I know there are many design flaws. Please point out as much as you can, as I want to improve my habits.

Thanks for reading :D

Am just a newb so this may be a stupid question, but i don't understand line 58

for (int b = -999; b < 1000; b = primes[index++])

why not

for (int b = -999; b < 1000; b++)

I've already done this but having looked at my code i'm not sure how(one off program so didn't comment on what i was doing) but both answers above are wrong although the second one is closer.

p.s what are your values for a and b?

Edited 5 Years Ago by frogboy77: n/a

Am just a newb so this may be a stupid question, but i don't understand line 58

for (int b = -999; b < 1000; b = primes[index++])

why not

for (int b = -999; b < 1000; b++)

I've already done this but having looked at my code i'm not sure how(one off program so didn't comment on what i was doing) but both answers above are wrong although the second one is closer.

p.s what are your values for a and b?

Weird... In Windows, it seems that my program just crashes! Anyways, a and b in the last few runs before it crashes:

a -414 b 0
a -414 b 0
a -414 b 0
a -414 b 0
a -414 b 0
a -414 b 0
a -414 b 0
a -414 b 0
a -413 b -999
a -412 b -999
a -412 b 0
a -412 b 0
a -411 b -999
a -410 b -999
a -409 b -999
a -409 b 0

Strange..

I did primes[index++] because after doing some thinking I found out that b always needs to be prime, or I'll be wasting computing time. b is assigned the value of primes[index] and only then index is incremented, which gives me an elegant way of moving through primes.

You're accessing the array out of bounds. In main, index exceeds 999, which is the largest available index in primes. Anything can happen, including wacky output and crashes.

what about something like

for (int a = -999; a < 1000; a++)	
{				
      for (int b = 0; b < put size of array here; b++)		
      {				x = 0; 			
                                    while (isPrime(x*x+a*x+primes[b]))
                                    x++; 			
                                    if (x > max)			
                                    {				
                                    max = x;				
                                    max_co = a;				
                                    max_co_sec = primes[b];			
                                    }		
       }	
} 	
cout << max_co * max_co_sec << endl;

Edited 5 Years Ago by frogboy77: n/a

Alright.. I rewrote a bit of my code according to the problems stated in this thread. I decided to use vectors for a safer approach. It seems that a and b have valid values, but the result is still wrong. Now I'm getting -9943.

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

vector<int> sieve(int i)
{
	bool primes[i];
	vector<int> primelist;
	
	primes[0] = false;
	primes[1] = false;
	
	for (int j = 2; j < i; j++)
		primes[j] = true;
	
	for (int j = 2; j * j < i; j++)
	{
		if (primes[j])
		{
			for (int k = j; k*j < i; k++)
				primes[k*j] = false;
		}
	}
	
	for (int k = 2; k < i; k++)
	{
		if (primes[k])
			primelist.push_back(k);
	}
	
	return primelist;
}

bool isPrime(int check)
{
	vector<int> primes = sieve(1000);
	
	if (find(primes.begin(), primes.end(), check) != primes.end())	
		return true;
		
	return false;
}

int main()
{
	vector<int> primes = sieve(1000);
	int max_co = 0;
	int max_co_sec = 0;
	int max = 0;
	int x;
	
	for (int a = -999; a < 1000; a++)
	{		
		for (int b = 0; b < primes.size(); b++)
		{	
			cout << "a " << a << " b " << b << endl;
			
			x = 0;
			
			while (isPrime(x*x+a*x+primes[b]))
				x++;
				
			if (x > max)
			{
				max = x;
				max_co = a;
				max_co_sec = b;
			}
		}
	}
	
	cout << max_co * max_co_sec << endl;
	
	return 0;
}

Thanks :D

> Considering quadratics of the form:

> n² + an + b, where |a| < 1000 and |b| < 1000

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

>> ... I found out that b always needs to be prime ...

This assumption does not appear to be correct. Counter-examples are easy to find:

1. with n == 7, a == 2, n² + an == 49 + 14 == 63, n² + an + b is prime

for b in [ -58, -56, -52, -50, -46, ..., -10, -4, +4, +8, +10, +16, +20, ... ]


2. with n == 5, a == -4, n² + an == 25 - 20 == 5, n² + an + b is prime

for b in [ -3, -2, 0, +2, +6, +8, +12, +14, +18, +24, ..... ]

Wow! I'm getting the correct answer now! Thanks vijayan! I guess my calculations were a bit off. Problem is, that my code now takes almost 30 seconds to run. Anyone got any useful optimization tips? I'm thinking about making the primes vector global, so I won't need to re-create it every time, but I don't think using global variables is very OOP-like. Any suggestions?

Edited 5 Years Ago by mrcpp: n/a

With 2 nested loops(pure brute force) and a simple prime checker mine's works in about 1s but i'm too new to this to suggest where the slow down is.
Read the forum it's usually full of good tips.
As you have the generation of primes done you should take a look at the most recent problem(315). No maths but a nice programming exercise.

>>I'm thinking about making the primes vector global, so I won't need to re-create it every time, but I don't think using global variables is very OOP-like. Any suggestions?
I don't think I would worry too much about how "OOP-like" it is at this point; you're not really doing any OOP in this program that I can see, it's almost all procedural. That being said though, I still wouldn't make it a global because you can't really effectively define it as a constant (you would have to know all the contents at compile time, which just isn't feasible).

More than likely, the slowdown is the continuous regeneration and destruction of the primes vector. I think I would generate it as early as possible in main() then make it a constant reference argument/parameter and pass it around to where it's needed. Keeping it local to main() increases the "security" and reduces the chance of getting it corrupted, as long as you perform your passing correctly.

Edited 5 Years Ago by Fbody: n/a

>>I'm thinking about making the primes vector global, so I won't need to re-create it every time, but I don't think using global variables is very OOP-like. Any suggestions?
I don't think I would worry too much about how "OOP-like" it is at this point; you're not really doing any OOP in this program that I can see, it's almost all procedural. That being said though, I still wouldn't make it a global because you can't really effectively define it as a constant (you would have to know all the contents at compile time, which just isn't feasible).

More than likely, the slowdown is the continuous regeneration and destruction of the primes vector. I think I would generate it as early as possible in main() then make it a constant reference argument/parameter and pass it around to where it's needed. Keeping it local to main() increases the "security" and reduces the chance of getting it corrupted, as long as you perform your passing correctly.

Thanks for the reply :D

I changed a few segments of my code according to what you said:

bool isPrime(int check, const vector<int>& primes)

and

const vector<int> primes = sieve(1000);

From there, whenever I called isPrime I just passed the primes vector. However, it still takes my program around 40 seconds to run! Do you happen to have some more useful optimization tips? Surely I am doing something wrong here if it is this slow.

I would have thought that would make a fairly significant difference by itself.

I think that, in order to facilitate better suggestions, it would be advisable for you to post your newest code. I don't know how much I can do in the optimization department, but I can certainly look. What compiler and O/S are you using?

Edited 5 Years Ago by Fbody: n/a

a possible speed up would be to generate the primes in a seperate program and output them to a file, then read them into your program into a vector and check that. it would mean only generating them once.

I would have thought that would make a fairly significant difference by itself.

I think that, in order to facilitate better suggestions, it would be advisable for you to post your newest code. I don't know how much I can do in the optimization department, but I can certainly look. What compiler and O/S are you using?

There has been an improvement. It runs a few seconds faster, now that I'm passing by reference. The sieve part isn't too crucial though, because it's already pretty optimized and runs very fast. I still think it could be faster using a bitset, but I'm not so sure on how to do that. I've never tried using bitsets before.

This is my newest code:

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

vector<int> sieve(int i)
{
	bool primes[i];
	vector<int> primelist;
	
	primes[0] = false;
	primes[1] = false;
	
	for (int j = 2; j < i; j++)
		primes[j] = true;
	
	for (int j = 2; j * j < i; j++)
	{
		if (primes[j])
		{
			for (int k = j; k*j < i; k++)
				primes[k*j] = false;
		}
	}
	
	for (int k = 2; k < i; k++)
	{
		if (primes[k])
			primelist.push_back(k);
	}
	
	return primelist;
}

bool isPrime(int check, const vector<int>& primes)
{
	if (find(primes.begin(), primes.end(), check) != primes.end())	
		return true;
		
	return false;
}

int main()
{
	const vector<int> primes = sieve(1000);
	int max_co = 0;
	int max_co_sec = 0;
	int max = 0;
	int x;
	
	for (int a = -999; a < 1000; a++)
	{		
		for (int b = 0; b < 1000; b++)
		{	
			cout << "a " << a << " b " << b << endl;
			
			x = 0;
			
			while (isPrime(x*x+a*x+b, primes))
				x++;
				
			if (x > max)
			{
				max = x;
				max_co = a;
				max_co_sec = b;
			}
		}
	}
	
	cout << max_co * max_co_sec << endl;
	
	return 0;
}

I'm using Ubuntu 10.10. Compiling with g++.


@frogboy77 I'm passing by reference, so it's generated only once.

Edited 5 Years Ago by mrcpp: n/a

In isprime can you not stop the search once the value in the vector is larger than the check value you send it?

I notice that you are returning the vector from sieve(). Try converting sieve() so that it accepts the primes vector as a reference parameter as well:

void sieve(vector<int> &, int);

int main() {
  vector<int> primes;
  sieve(primes, 1000);

  //etc...

Otherwise, I don't really see anything glaringly bad. I think it's just sheer volume. You're talking about 2-million iterations with each iteration having an unknown number of while loop iterations within that. How new/modern is your hardware?

Edited 5 Years Ago by Fbody: n/a

I notice that you are returning the vector from sieve(). Try converting sieve() so that it accepts the primes vector as a reference parameter as well:

void sieve(vector<int> &, int);

int main() {
  vector<int> primes;
  sieve(primes, 1000);

  //etc...

Otherwise, I don't really see anything glaringly bad. I think it's just sheer volume. You're talking about 2-million iterations with each iteration having an unknown number of while loop iterations within that. How new/modern is your hardware?

I changed the sieve function as you said. I also switched to iterators in the isPrime function so I could stop whenever *iterator exceeds the function's given value, as frogboy77 suggested.

Updated code:

#include <iostream>
#include <vector>

using namespace std;

void sieve(vector<int>& primelist, int i)
{
	bool primes[i];
	
	primes[0] = false;
	primes[1] = false;
	
	for (int j = 2; j < i; j++)
		primes[j] = true;
	
	for (int j = 2; j * j < i; j++)
	{
		if (primes[j])
		{
			for (int k = j; k*j < i; k++)
				primes[k*j] = false;
		}
	}
	
	for (int k = 2; k < i; k++)
	{
		if (primes[k])
			primelist.push_back(k);
	}
}

bool isPrime(int check, const vector<int>& primes)
{	
	for (vector<int>::const_iterator it = primes.begin(); *it < check + 1; ++it)    //check+1 and not <= check because addition is much faster than comparison
	{
		if (*it == check)
			return true;
	}
		
	return false;
}

int main()
{
	vector<int> primes;
	int max_co = 0;
	int max_co_sec = 0;
	int max = 0;
	int x;
	
	sieve(primes, 1000);
	
	for (int a = -999; a < 1000; a++)
	{		
		for (int b = 0; b < 1000; b++)
		{	
			cout << "a " << a << " b " << b << endl;
			
			x = 0;
			
			while (isPrime(x*x+a*x+b, primes))
				x++;
				
			if (x > max)
			{
				max = x;
				max_co = a;
				max_co_sec = b;
			}
		}
	}
	
	cout << max_co * max_co_sec << endl;
	
	return 0;
}

I have an Intel Core 2 Duo 3GHZ processor. I bought my PC about 2 years ago, so all my hardware is pretty new, except for my video card, which I bought last year.

Edited 5 Years Ago by mrcpp: n/a

>> There has been an improvement. It runs a few seconds faster ...

Console output is slow, you might simply delete the line

cout << "a " << a << " b " << b << endl;

altogether.

Anyone got any useful optimization tips?

for (vector<int>::const_iterator it = primes.begin(); *it < check + 1; ++it) performs a linear search in O(N) time where N == primes.size()

As the primes are in sorted order, you can reduce the time to O(log N) by doing a binary search.

Better still, you can do this in O(1) (constant time) by indexing into the vector.

For example:

#include <vector>
#include <cassert>
#include <iostream>

std::vector<bool> make_primes( std::size_t N ) // primes < N
{
   std::vector<bool> sieve( N, true ) ;
   sieve[0] = sieve[1] = false ;

   for( std::size_t i = 2; i*i < N ; ++i ) if( sieve[i] )
      for( std::size_t j = i + i ; j < N ; j += i ) sieve[j] = false ;

   return sieve ;
}

inline bool is_prime( int n, const std::vector<bool>& primes )
{ return n>1 && primes[n] ; }

int main ()
{
    enum { AMAX = 1000, BMAX = 1000, NMAX = 200 } ;
    const std::vector<bool>& primes = make_primes( NMAX*NMAX + AMAX*NMAX + BMAX + 1 ) ;

    int max_axb = 0 ;
    int max_n = 0;

    for( int a = -999 ; a < 1000 ; ++a )
    {
        for( int b = -999 ; b < 1000 ; ++b )
        {
            int n = 0 ;
            while( is_prime( n*n + a*n + b, primes ) ) ++n ;
            if( n > max_n )
            {
                max_n = n ;
                max_axb = a * b ;
            }
        }
    }

    assert( max_n <= NMAX ) ;
    std::cout << max_axb << '\n' ;
}

On my laptop (FreeBSD 8.1 RELEASE, Core2 Duo @1.7 GHz, GCC 4.5.2), I get:

> g++45 -Wall -std=c++98 -pedantic -Werror -O3 euler27.cc && ./a.out
-59231
0.087u 0.000s 0:00.08 100.0% 5+1091k 0+0io 0pf+0w

On different runs, the time varied between 0.021u and 0.108u

Edited 5 Years Ago by vijayan121: n/a

for (vector<int>::const_iterator it = primes.begin(); *it < check + 1; ++it) performs a linear search in O(N) time where N == primes.size()

As the primes are in sorted order, you can reduce the time to O(log N) by doing a binary search.

Better still, you can do this in O(1) (constant time) by indexing into the vector.

For example:

#include <vector>
#include <cassert>
#include <iostream>

std::vector<bool> make_primes( std::size_t N ) // primes < N
{
   std::vector<bool> sieve( N, true ) ;
   sieve[0] = sieve[1] = false ;

   for( std::size_t i = 2; i*i < N ; ++i ) if( sieve[i] )
      for( std::size_t j = i + i ; j < N ; j += i ) sieve[j] = false ;

   return sieve ;
}

inline bool is_prime( int n, const std::vector<bool>& primes )
{ return n>1 && primes[n] ; }

int main ()
{
    enum { AMAX = 1000, BMAX = 1000, NMAX = 200 } ;
    const std::vector<bool>& primes = make_primes( NMAX*NMAX + AMAX*NMAX + BMAX + 1 ) ;

    int max_axb = 0 ;
    int max_n = 0;

    for( int a = -999 ; a < 1000 ; ++a )
    {
        for( int b = -999 ; b < 1000 ; ++b )
        {
            int n = 0 ;
            while( is_prime( n*n + a*n + b, primes ) ) ++n ;
            if( n > max_n )
            {
                max_n = n ;
                max_axb = a * b ;
            }
        }
    }

    assert( max_n <= NMAX ) ;
    std::cout << max_axb << '\n' ;
}

On my laptop (FreeBSD 8.1 RELEASE, Core2 Duo @1.7 GHz, GCC 4.5.2), I get:

On different runs, the time varied between 0.021u and 0.108u

Wow, awesome! Thanks! Great code!

Could you just explain a few things? I'm not sure I quite understood everything:

inline bool is_prime( int n, const std::vector<bool>& primes )
{ return n>1 && primes[n] ; }

Why did you choose inline here, and why do you return n>1 && primes[n] ?

enum { AMAX = 1000, BMAX = 1000, NMAX = 200 } ;
    const std::vector<bool>& primes = make_primes( NMAX*NMAX + AMAX*NMAX + BMAX + 1 ) ;

Why did you enumerate AMAX, BMAX, and NMAX? Why not just do make_primes(241001) ?


And a few general things I want to understand better:

Why do you call make_primes and pass such a large value? And why is size_t used throughout the function? Is there something concerning enumerations and size_t that I'm missing here? Also, what does assert(max_n <= NMAX) do? An explanation would be really great.


Thanks for the reply! :D


edit:

Oh and I almost forgot.. Could you please explain why you used the compiling options you used? I'm not familiar with most of them.

Edited 5 Years Ago by mrcpp: n/a

Why did you choose inline here

is_prime() is a small (code size) function that is called millions of times.
See: http://en.wikipedia.org/wiki/Inline_expansion

and why do you return n>1 && primes[n] ?

primes is a vector of bools where primes[n] == true if and only if n is prime.
For example, primes[13] == true and primes[18] == false n>1 && primes[n] => if n is greater than one (to take care of negative n), and primes[n] is true, then n is prime.

Why did you enumerate AMAX, BMAX, and NMAX? Why not just do make_primes(241001) ?

The idea was to make the code a bit more easy to modify; for example if the problem statement changed from
the current n² + an + b, where |a| < 1000 and |b| < 1000 to, say, n² + an + b, where |a| < 500 and |b| < 20000 .

Of course I messed up by writing

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

instead of

for( int a = 1 - AMAX ; a < AMAX ; ++a )
{
    for( int b = 1 - BMAX ; b < BMAX ; ++b )

Why do you call make_primes and pass such a large value?

It is the largest value (plus one) that the quadratric n² + an + b, where |a| <= AMAX and |b| <= BMAX can generate;
we need to be able to check for primes upto that value.

And why is size_t used throughout the function? Is there something concerning enumerations and size_t that I'm missing here?

No significance, could have just used an int instead. The original idea was to avoid a compiler warning: 'comparison between signed and unsigned integer expressions', and I was too lazy to type std::vector<bool>::size_type .

what does assert(max_n <= NMAX) do?

The whole code is based on the assumption that n will never exceed NMAX.
If it does, we want to know about it so that we can revise the value of NMAX upwards.
See: http://en.wikipedia.org/wiki/Assertion_%28computing%29
http://www.stanford.edu/~pgbovine/programming-with-asserts.htm
http://en.wikipedia.org/wiki/Assert.h

Could you please explain why you used the compiling options you used?

-std=c++98 - conform to the IS (International Standard) for C++ -pedantic - reject programs that do not conform strictly to ISO C++, and generate all the warnings required by the IS. -Wall - enable all the warnings about easily avoidable questionable constructs. ( -Wextra is also a handy option.) -Werror - treat all warnings as errors. -O3 - turn on a bunch of compiler optimizations.
For details, see: http://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#Optimize-Options

For doumentation on all the GCC options, see: http://gcc.gnu.org/onlinedocs/gcc/Invoking-GCC.html

Comments
Agreed
comprehensive answer

is_prime() is a small (code size) function that is called millions of times.
See: http://en.wikipedia.org/wiki/Inline_expansion


primes is a vector of bools where primes[n] == true if and only if n is prime.
For example, primes[13] == true and primes[18] == false n>1 && primes[n] => if n is greater than one (to take care of negative n), and primes[n] is true, then n is prime.


The idea was to make the code a bit more easy to modify; for example if the problem statement changed from
the current n² + an + b, where |a| < 1000 and |b| < 1000 to, say, n² + an + b, where |a| < 500 and |b| < 20000 .

Of course I messed up by writing

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

instead of

for( int a = 1 - AMAX ; a < AMAX ; ++a )
{
    for( int b = 1 - BMAX ; b < BMAX ; ++b )

It is the largest value (plus one) that the quadratric n² + an + b, where |a| <= AMAX and |b| <= BMAX can generate;
we need to be able to check for primes upto that value.


No significance, could have just used an int instead. The original idea was to avoid a compiler warning: 'comparison between signed and unsigned integer expressions', and I was too lazy to type std::vector<bool>::size_type .


The whole code is based on the assumption that n will never exceed NMAX.
If it does, we want to know about it so that we can revise the value of NMAX upwards.
See: http://en.wikipedia.org/wiki/Assertion_%28computing%29
http://www.stanford.edu/~pgbovine/programming-with-asserts.htm
http://en.wikipedia.org/wiki/Assert.h -std=c++98 - conform to the IS (International Standard) for C++ -pedantic - reject programs that do not conform strictly to ISO C++, and generate all the warnings required by the IS. -Wall - enable all the warnings about easily avoidable questionable constructs. ( -Wextra is also a handy option.) -Werror - treat all warnings as errors. -O3 - turn on a bunch of compiler optimizations.
For details, see: http://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html#Optimize-Options

For doumentation on all the GCC options, see: http://gcc.gnu.org/onlinedocs/gcc/Invoking-GCC.html

Thank you very much!

This question has already been answered. Start a new discussion instead.