Need help with a frustrating code - Getting the wrong output

[mrcpp]

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

[frogboy77]

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?

[mrcpp]

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.

[Narue]

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.

[frogboy77]

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;

[mrcpp]

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

[vijayan121]

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

[mrcpp]

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?

[Narue]

>Wow! I'm getting the correct answer now!
...

[frogboy77]

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.

[Fbody]

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

[mrcpp]

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

[Fbody]

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?

[frogboy77]

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.

[mrcpp]

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.

[frogboy77]

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

[Fbody]

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?

[mrcpp]

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.

[mitrmkar]

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

[vijayan121]

Anyone got any useful optimization tips?

for (vector<int>::const_iterator it = primes.begin(); *it < check + 1; ++it)

performs a linear search inO(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