Â asked on # Prime number counter

I got really bored one day so I wrote this program that counts up prime numbers. First is C, second is TI-83 language. My challege is for people to find faster ways to find Prime numbers.

// *************************************************

// File Name:prime_1.c Â

// Author: Samuel Dockery

// UserID: dockesa

// Date Created: September 28, 2002 Â

// Purpose: Compute Prime Numbers and Display them to stdout. Â

//

// Notes: Compiled using Microsoft Visual C++ 6.0

// *************************************************

// Include Files

#include <stdio.h> Â

// Global Constants

// None

// Global Types

// None

// Global Variables

// None

// Function Prototypes

// None

//****************************************************

// Function name:main Â

// Author: Samuel Dockery

// Purpose: Compute Prime Numbers and Display them to stdout.

// Â

// Input parameters: none Â

// Output parameters: none Â

// Returns:0 Â

// Error handling: none Â

// Record of changes: none

//****************************************************

int main (void)

{

long halfPrime = 1;//Used to speed up process

unsigned long currentPrime = 2;//Current Prime number computer is calculating.

unsigned long smallNumber = 2;//Number used in checking to see if currentPrime is a prime number.

for (;;)//Loop infinitly.

Â Â Â {

Â Â Â Â Â if ( (currentPrime % smallNumber) != 0)//Checks to see if currentPrime is not evenly divisible.

Â Â Â Â Â Â Â Â { Â Â Â Â Â

Â Â Â Â Â Â Â Â smallNumber++;//Increment smallNumber by 1.

Â Â Â Â Â Â Â Â halfPrime = currentPrime/2;

Â Â Â Â Â Â Â Â if (smallNumber >= halfPrime)//If smallNumber == currentPrime then it is a prime number.

Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â printf("%li\n", currentPrime);//Print prime number.

Â Â Â Â Â Â Â Â Â Â currentPrime++;//Increment currentPrime by 1.

Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â }//End If

Â Â Â Â Â Â Â Â }

Â Â Â Â Â else//Do following statement becuase currentPrime is not Prime.

Â Â Â Â Â Â Â Â { Â Â

Â Â Â Â Â Â Â Â currentPrime++;//Increment currentPrime by 1.

Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â }//End If

Â Â Â } Â Â

return (0);// Return a zero value from the function to indicate successful processing

} // End main

--I also translated it into the TI-83 calculator language.--

ClrHome

1.e99->M

2->C

1->S

While C is less than or equal to M

If fPart(C/S)>0

Then

S+1->S

If (S=C)

Then

Disp C

C+1->C

2->S

End

Else

C+1->C

2->S

End

End

Stop

// **************************

// File Name:prime_1.c Â

// Author: Samuel Dockery

// UserID: dockesa

// Date Created: September 28, 2002 Â

// Purpose: Compute Prime Numbers and Display them to stdout. Â

//

// Notes: Compiled using Microsoft Visual C++ 6.0

// **************************

// Include Files

#include <stdio.h> Â

// Global Constants

// None

// Global Types

// None

// Global Variables

// None

// Function Prototypes

// None

//************************

// Function name:main Â

// Author: Samuel Dockery

// Purpose: Compute Prime Numbers and Display them to stdout.

// Â

// Input parameters: none Â

// Output parameters: none Â

// Returns:0 Â

// Error handling: none Â

// Record of changes: none

//************************

int main (void)

{

long halfPrime = 1;//Used to speed up process

unsigned long currentPrime = 2;//Current Prime number computer is calculating.

unsigned long smallNumber = 2;//Number used in checking to see if currentPrime is a prime number.

for (;;)//Loop infinitly.

Â Â Â {

Â Â Â Â Â if ( (currentPrime % smallNumber) != 0)//Checks to see if currentPrime is not evenly divisible.

Â Â Â Â Â Â Â Â { Â Â Â Â Â

Â Â Â Â Â Â Â Â smallNumber++;//Increment smallNumber by 1.

Â Â Â Â Â Â Â Â halfPrime = currentPrime/2;

Â Â Â Â Â Â Â Â if (smallNumber >= halfPrime)//If smallNumber == currentPrime then it is a prime number.

Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â printf("%li\n", currentPrime);//Print prime number.

Â Â Â Â Â Â Â Â Â Â currentPrime++;//Increment

Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â }//End If

Â Â Â Â Â Â Â Â }

Â Â Â Â Â else//Do following statement becuase currentPrime is not Prime.

Â Â Â Â Â Â Â Â { Â Â

Â Â Â Â Â Â Â Â currentPrime++;//Increment

Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â }//End If

Â Â Â } Â Â

return (0);// Return a zero value from the function to indicate successful processing

} // End main

--I also translated it into the TI-83 calculator language.--

ClrHome

1.e99->M

2->C

1->S

While C is less than or equal to M

If fPart(C/S)>0

Then

S+1->S

If (S=C)

Then

Disp C

C+1->C

2->S

End

Else

C+1->C

2->S

End

End

Stop

Math / Science

I just realised that I am an idiot!!!! current will NEVER loop back to 0, it will loop back to 1 and factor should start at 2...so....

main()

{

Â unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

Â unsigned long factor, square;

Â char prime;

Â while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

Â {

Â Â square = sqrt( current ) + 1; Â // function to get the square root (add one to allow for rounding)

Â Â factor = 2;

Â Â prime = 1;

Â Â while( (factor <Â square) &&Â prime ) Â // at this point prime means "possibly prime" :)

Â Â {

Â Â Â if( factor % current )

Â Â Â Â factor++;

Â Â Â else

Â Â Â Â prime = 0; Â Â Â Â Â Â Â Â Â Â Â // Evenly divisible = not prime. D'OH!

Â Â }

Â Â if( prime )

Â Â Â printf("Prime number found ! %d\n", current);

Â Â current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

Â }

}

main()

{

Â unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

Â unsigned long factor, square;

Â char prime;

Â while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

Â {

Â Â square = sqrt( current ) + 1; Â // function to get the square root (add one to allow for rounding)

Â Â factor = 2;

Â Â prime = 1;

Â Â while( (factor <Â square) &&Â prime ) Â // at this point prime means "possibly prime" :)

Â Â {

Â Â Â if( factor % current )

Â Â Â Â factor++;

Â Â Â else

Â Â Â Â prime = 0; Â Â Â Â Â Â Â Â Â Â Â // Evenly divisible = not prime. D'OH!

Â Â }

Â Â if( prime )

Â Â Â printf("Prime number found ! %d\n", current);

Â Â current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

Â }

}

Computational complexity (and hence speed) for finding primes changes with the size. Â For example, one algorithm will be faster for the first 100,000 primes and then a different algorithm will be faster for the next. Â

(To TheBeaver) By the way, even numbers can be prime. Â 2 is the prime even number. Â All other primes are odd.

(To TheBeaver) By the way, even numbers can be prime. Â 2 is the prime even number. Â All other primes are odd.

Unreal7000,

Â There are lots of very very slick ways of finding primes for instance the Miller-Rabin pseudo prime will eliminate most non-primes.

Â There are lots of very very slick ways of finding primes for instance the Miller-Rabin pseudo prime will eliminate most non-primes.

main(){

Â Â int p,d,s;

Â Â printf("%d\n",2);

Â Â for( p=3,s=1;;p+=2 ){

Â Â Â Â if( s*s <Â p ){ s++; }

Â Â Â Â for( d=3;d<=s;d+=2 ){

Â Â Â Â Â Â if( !(p%d) ){ d=s+3; }

Â Â Â Â }

Â Â Â Â if( d <Â s+3 ){

Â Â Â Â Â Â printf("%d\n",p);

Â Â Â Â }

Â Â }

}

Â Â int p,d,s;

Â Â printf("%d\n",2);

Â Â for( p=3,s=1;;p+=2 ){

Â Â Â Â if( s*s <Â p ){ s++; }

Â Â Â Â for( d=3;d<=s;d+=2 ){

Â Â Â Â Â Â if( !(p%d) ){ d=s+3; }

Â Â Â Â }

Â Â Â Â if( d <Â s+3 ){

Â Â Â Â Â Â printf("%d\n",p);

Â Â Â Â }

Â Â }

}

Heres some new code, and its over 37 times faster than the first one I posted. Read the History at the top of code for all new things done.

// *************************************************

// File Name:prime_2.c Â

// Author: Samuel Dockery

// UserID: dockesa

// Date Created: September 28, 2002

// Date Revised: February 23, 2003 Â

// Purpose: Compute Prime Numbers and Display them to stdout. Â

//

// History:

// -Revision 2- Â February 23, 2003

// Â Â Â Â Â Â Â Â -Took halfPrime to the next step.

// Â Â Â Â Â Â Â Â -Took halfPrime out, replaced by currentPrimeValue.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -halfPrime is half, so now currentPrimeVale is third, fourth, fifth, etc...

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Added primeNumberCounter but I have it commented out

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â because it slows the whole thing down significantly.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Maybe I will work this in in Revision 3.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Benchmark to a value of 500,000 on AMD Athlon 1400Mhz:

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 1 Results: 3min. 6 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 2 Results: 0min. 5 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â WOW now Revision 2 is over 37 times faster than Revision 1.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â But as numbers get higher I imagine this value increases (Maybe).

// -Revision 1.1- September 28, 2003

// Â Â Â Â Â Â Â Â -Added halfPrime to speed up processing.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Ex. 101/51, 101/52, is pointless so throw the rest out (53,54,55,etc...).

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -It increases speed but its less than 2 times faster than initial release.

// -Revision 1- Â September 28, 2003

// Â Â Â Â Â Â Â Â -Initial Release.

// *************************************************

// Include Files

#include <stdio.h> Â

// Global Constants

// None

// Global Types

// None

// Global Variables

// None

// Function Prototypes

// None

//****************************************************

// Function name:main Â

// Author: Samuel Dockery

// Purpose: Compute Prime Numbers and Display them to stdout.

// Â

// Input parameters: none Â

// Output parameters: none Â

// Returns:0 Â

// Error handling: none Â

// Record of changes: Revision 2

//****************************************************

Â

int main (void)

{

unsigned long currentPrimeDivision = 0;//Used to speed up processing.

unsigned long currentPrimeCounter = 1;//Tell me what the current prime number count is.

unsigned long currentPrime = 2;//Current Prime number computer is calculating.

unsigned long smallNumber = 2;//Number used in checking to see if currentPrime is a prime number.

printf("2\n");//smallest prime number.

//printf(" Â Â Â Â Â Â Â Â Â Â Â Â 1\n");//Tells that this is first prime number.

for (;;)//Loop infinitly.

Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â if ( (currentPrime % smallNumber) != 0)//Checks to see if currentPrime is not evenly divisible.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = currentPrime / smallNumber;//Assign currentPrimeDivision a value.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â if (currentPrimeDivision <Â smallNumber)//Verify that it is Prime.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â printf("%li\n", currentPrime);//Print prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //currentPrimeCounter++;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //printf(" Â Â Â Â Â Â Â Â Â Â Â Â %li\n",currentPrimeCounter);//Displays what prime number this is.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;//Speeds up process.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime++;//Increment currentPrime by 1.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â else

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber++;//Keep going it could be a possible prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â else//Do following statement becuase currentPrime is not Prime.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime++;//Increment currentPrime by 1.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }//End If

Â Â Â Â Â Â } Â Â Â Â Â

return (0);// Return 0 for successful processing. You probably won't get this far.

} // End main

// **************************

// File Name:prime_2.c Â

// Author: Samuel Dockery

// UserID: dockesa

// Date Created: September 28, 2002

// Date Revised: February 23, 2003 Â

// Purpose: Compute Prime Numbers and Display them to stdout. Â

//

// History:

// -Revision 2- Â February 23, 2003

// Â Â Â Â Â Â Â Â -Took halfPrime to the next step.

// Â Â Â Â Â Â Â Â -Took halfPrime out, replaced by currentPrimeValue.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -halfPrime is half, so now currentPrimeVale is third, fourth, fifth, etc...

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Added primeNumberCounter but I have it commented out

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â because it slows the whole thing down significantly.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Maybe I will work this in in Revision 3.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Benchmark to a value of 500,000 on AMD Athlon 1400Mhz:

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 1 Results: 3min. 6 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 2 Results: 0min. 5 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â WOW now Revision 2 is over 37 times faster than Revision 1.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â But as numbers get higher I imagine this value increases (Maybe).

// -Revision 1.1- September 28, 2003

// Â Â Â Â Â Â Â Â -Added halfPrime to speed up processing.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Ex. 101/51, 101/52, is pointless so throw the rest out (53,54,55,etc...).

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -It increases speed but its less than 2 times faster than initial release.

// -Revision 1- Â September 28, 2003

// Â Â Â Â Â Â Â Â -Initial Release.

// **************************

// Include Files

#include <stdio.h> Â

// Global Constants

// None

// Global Types

// None

// Global Variables

// None

// Function Prototypes

// None

//************************

// Function name:main Â

// Author: Samuel Dockery

// Purpose: Compute Prime Numbers and Display them to stdout.

// Â

// Input parameters: none Â

// Output parameters: none Â

// Returns:0 Â

// Error handling: none Â

// Record of changes: Revision 2

//************************

Â

int main (void)

{

unsigned long currentPrimeDivision = 0;//Used to speed up processing.

unsigned long currentPrimeCounter = 1;//Tell me what the current prime number count is.

unsigned long currentPrime = 2;//Current Prime number computer is calculating.

unsigned long smallNumber = 2;//Number used in checking to see if currentPrime is a prime number.

printf("2\n");//smallest prime number.

//printf(" Â Â Â Â Â Â Â Â Â Â Â Â 1\n");//Tells that this is first prime number.

for (;;)//Loop infinitly.

Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â if ( (currentPrime % smallNumber) != 0)//Checks to see if currentPrime is not evenly divisible.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = currentPrime / smallNumber;//Assign currentPrimeDivision a value.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â if (currentPrimeDivision <Â smallNumber)//Verify that it is Prime.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â printf("%li\n", currentPrime);//Print prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //currentPrimeCounter++;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //printf(" Â Â Â Â Â Â Â Â Â Â Â Â %li\n",currentPrimeCounter

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;//Speeds up process.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime++;//Increment

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â else

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber++;//Keep going it could be a possible prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â else//Do following statement becuase currentPrime is not Prime.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime++;//Increment

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }//End If

Â Â Â Â Â Â } Â Â Â Â Â

return (0);// Return 0 for successful processing. You probably won't get this far.

} // End main

main()

{

unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

unsigned long factor, square;

char prime;

while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

{

Â square = sqrt( current ) + 1; Â // function to get the square root (add one to allow for rounding)

Â factor = 3;

Â prime = 1;

Â while( (factor <Â square) &&Â prime ) Â // at this point prime means "possibly prime" :)

Â {

Â Â if( factor % current )

Â Â Â factor+=2;

Â Â else

Â Â Â prime = 0; Â Â Â Â Â Â Â Â Â Â Â // Evenly divisible = not prime. D'OH!

Â }

Â if( prime )

Â Â printf("Prime number found ! %d\n", current);

Â current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

}

}

Made one change...factors now increment by 2 as I won't ever have even factors since I am not testing even numbers.

{

unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

unsigned long factor, square;

char prime;

while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

{

Â square = sqrt( current ) + 1; Â // function to get the square root (add one to allow for rounding)

Â factor = 3;

Â prime = 1;

Â while( (factor <Â square) &&Â prime ) Â // at this point prime means "possibly prime" :)

Â {

Â Â if( factor % current )

Â Â Â factor+=2;

Â Â else

Â Â Â prime = 0; Â Â Â Â Â Â Â Â Â Â Â // Evenly divisible = not prime. D'OH!

Â }

Â if( prime )

Â Â printf("Prime number found ! %d\n", current);

Â current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

}

}

Made one change...factors now increment by 2 as I won't ever have even factors since I am not testing even numbers.

I wrote a nice program for this a long time ago. It was in VB, but the logic is perfectly portable to any language. I'll post it as soon as I get home.

There an extremely fast probabilistic primality test routine due to Miller and Rabin. Scaled to ``long'' arithmetics such an algorithm is used for real-life cryptographic applications. The version below works with 32-bit integers, and simply prints a sequence of primes to stdout.

/*

Â * Miller-Rabin probabilistic primality test

Â * Probability of false positives is at most 4^(-ROUNDS)

Â */

#define ROUNDS 20

#define TRUE 1

#define FALSE 0

#include <stdlib.h>

#include <time.h>

#include <stdio.h>

typedef unsigned int uint;

uint mulmod(uint x, uint y, uint m) {

Â Â return (x*y)%m;

}

uint powmod(uint x, uint a, uint m) {

Â Â uint r = 1;

Â Â for(;a>0;a>>=1) {

Â Â Â Â if(a&1) r = mulmod(r,x,m);

Â Â Â Â x = mulmod(x,x,m);

Â Â }

Â return r;

}

int isprime(uint p) {

Â Â uint q,i,a;

Â Â if(p<=1) return FALSE;

Â Â if((p&1)==0) return (p==2); Â Â Â /* even p */

Â Â for(q=p-1;(q&1)==0;q>>=1);

Â Â for(i=0;i<ROUNDS;i++) {

Â Â Â Â uint a = 2 + (((rand()<<15)|rand())%(p-2));

Â Â Â Â if(powmod(a,p-1,p)!=1) return FALSE;

Â Â Â Â a = powmod(a,q,p);

Â Â Â Â if(a!=1) {

Â Â Â Â Â Â while(a!=1 &&Â a!=p-1) a=mulmod(a,a,p);

Â Â Â Â Â Â if(a==1) return FALSE;

Â Â Â Â }

Â Â }

Â Â return TRUE;

}

void main() {

Â Â uint t;

Â Â srand(time(0));

Â Â for(t=1;;t++) if(isprime(t)) printf("%d ",t);

}

/*

Â * Miller-Rabin probabilistic primality test

Â * Probability of false positives is at most 4^(-ROUNDS)

Â */

#define ROUNDS 20

#define TRUE 1

#define FALSE 0

#include <stdlib.h>

#include <time.h>

#include <stdio.h>

typedef unsigned int uint;

uint mulmod(uint x, uint y, uint m) {

Â Â return (x*y)%m;

}

uint powmod(uint x, uint a, uint m) {

Â Â uint r = 1;

Â Â for(;a>0;a>>=1) {

Â Â Â Â if(a&1) r = mulmod(r,x,m);

Â Â Â Â x = mulmod(x,x,m);

Â Â }

Â return r;

}

int isprime(uint p) {

Â Â uint q,i,a;

Â Â if(p<=1) return FALSE;

Â Â if((p&1)==0) return (p==2); Â Â Â /* even p */

Â Â for(q=p-1;(q&1)==0;q>>=1);

Â Â for(i=0;i<ROUNDS;i++) {

Â Â Â Â uint a = 2 + (((rand()<<15)|rand())%(p-

Â Â Â Â if(powmod(a,p-1,p)!=1) return FALSE;

Â Â Â Â a = powmod(a,q,p);

Â Â Â Â if(a!=1) {

Â Â Â Â Â Â while(a!=1 &&Â a!=p-1) a=mulmod(a,a,p);

Â Â Â Â Â Â if(a==1) return FALSE;

Â Â Â Â }

Â Â }

Â Â return TRUE;

}

void main() {

Â Â uint t;

Â Â srand(time(0));

Â Â for(t=1;;t++) if(isprime(t)) printf("%d ",t);

}

Some C compilers cannot produce correct result for ``mulmod'' function since an intermediate value requires more than 32 bits. This is a fixed function:

uint mulmod(uint x, uint y, uint m) {

Â Â return ((unsigned long long)x*(unsigned long long)y)%m;

}

uint mulmod(uint x, uint y, uint m) {

Â Â return ((unsigned long long)x*(unsigned long long)y)%m;

}

I'm currently testing the following numbers for primality :-

2^15515089-1 15.72% done

2^16752529-1 29.85% done

2^8051543-1 Â 13.56% done

see www.mersenne.orgÂ and join up!

2^15515089-1 15.72% done

2^16752529-1 29.85% done

2^8051543-1 Â 13.56% done

see www.mersenne.orgÂ and join up!

relf my code and yours seem to run at the same speed.

What is the smallest prime? I've heard 2 was but others say it is not.

Unreal7000, as soon as numbers get larger, your method will take much more time.

The smallest prime number is 2. See

http://mathworld.wolfram.com/PrimeNumber.html

The smallest prime number is 2. See

http://mathworld.wolfram.com/PrimeNumber.html

if you want to display prime numbers, i humbly suggest that you use

if integer%all integers before it != 0, then it must be prime. (%=modulus operator, calculates the remainder of a number.)

make a loop that divides a number Â by each of all the lesser integers.

e.g. integer = 17,

loop:

17/16

17/15

17/14

.

.

.

17/2

find the results. if at any point in this loop the remainder becomes zero, then the number is composite and is thereby eliminated.

but if it can finish the loop up to the divisor=2, then it must be prime.display the prime number!

increment the counter variable

decriment the divisor.

Sorry if i can't make the algorithm right now. but i think you'd understand what i wrote just fine.

Regards,

Centeris

if integer%all integers before it != 0, then it must be prime. (%=modulus operator, calculates the remainder of a number.)

make a loop that divides a number Â by each of all the lesser integers.

e.g. integer = 17,

loop:

17/16

17/15

17/14

.

.

.

17/2

find the results. if at any point in this loop the remainder becomes zero, then the number is composite and is thereby eliminated.

but if it can finish the loop up to the divisor=2, then it must be prime.display the prime number!

increment the counter variable

decriment the divisor.

Sorry if i can't make the algorithm right now. but i think you'd understand what i wrote just fine.

Regards,

Centeris

I had a quick way where you use the square root of the number. I'll try and find it and post.

I think beaver already did that, so I won't do it.

Sorry about the psuedocode, I haven't actually coded in about a year.

CAP=1000;

c;

m;

array primes[CAP+1];

primes[1]=0;

for(c=2,c<=CAP,c++)

{

Â Â primes[c]=1;

Â Â for(m=2, m*c<=CAP,m++)

Â Â {

Â Â Â primes[m*c]=0;

Â Â }

}

for(c=1,c<=CAP,c++)

{

Â Â if(primes[c]=1)

Â Â {

Â Â Â print c;

Â Â }

}

Should beat order N^2. Â CAP is the largest number you want to calculate the primes to.

Let me know if this is any faster!

CAP=1000;

c;

m;

array primes[CAP+1];

primes[1]=0;

for(c=2,c<=CAP,c++)

{

Â Â primes[c]=1;

Â Â for(m=2, m*c<=CAP,m++)

Â Â {

Â Â Â primes[m*c]=0;

Â Â }

}

for(c=1,c<=CAP,c++)

{

Â Â if(primes[c]=1)

Â Â {

Â Â Â print c;

Â Â }

}

Should beat order N^2. Â CAP is the largest number you want to calculate the primes to.

Let me know if this is any faster!

spikemike essentially presented Eratosthenes Prime Number Sieve but not in the optimal form. Optimized version follows.

The major lack of Eratosthenes Prime Number Sieve is memory requirement proportional to maximal prime that we want to reach.

#include <iostream>

#define MAXNUM 1000000

bool prime[MAXNUM];

int main() {

Â // initialize with ``true'' values

Â for(int i=2;i<MAXNUM;i++) prime[i] = true;

Â

Â // sieving

Â for(int i=2;i*i<MAXNUM;i++) if(prime[i]) {

Â Â for(int j=i*i;j<MAXNUM;j+=i) prime[j]=false;

Â }

Â // output primes

Â for(int i=2;i<MAXNUM;i++) if(prime[i])

Â Â std::cout <<Â " "Â << i;

Â

Â Â return 0;

}

The major lack of Eratosthenes Prime Number Sieve is memory requirement proportional to maximal prime that we want to reach.

#include <iostream>

#define MAXNUM 1000000

bool prime[MAXNUM];

int main() {

Â // initialize with ``true'' values

Â for(int i=2;i<MAXNUM;i++) prime[i] = true;

Â

Â // sieving

Â for(int i=2;i*i<MAXNUM;i++) if(prime[i]) {

Â Â for(int j=i*i;j<MAXNUM;j+=i) prime[j]=false;

Â }

Â // output primes

Â for(int i=2;i<MAXNUM;i++) if(prime[i])

Â Â std::cout <<Â " "Â << i;

Â

Â Â return 0;

}

Centeris

But there is no point in dividing for example 100/99, 100/98, 100/97. Numbers higher than half the prime number can never be a possibility in looking for a prime.

The way mine works is lets say you want to find out if 97 is a prime number. It can do it in 8 steps.

97/2 = 48.5 ok prime so far.

97/3 = 32.3 ok prime so far.

97/4 = 24.24 ok prime so far.

97/5 = 19.4 ok prime so far.

97/6 = 16.16 ok prime so far.

97/7 = 13.85 ok prime so far.

97/8 = 12.125 ok prime so far.

97/9 = 10.77 OK Its Prime!

Now I know there is not point in dividing by 4,6,8 because there even, but for some reason it would run slower if I took those out of the code. Go figure.

If the number to the right of the equals sign ever has a decimal its a possiblity the number we are looking for is prime. And we know to stop dividing when the right hand number is bigger than left hand denominator number. IE 10.77 is bigger than 9.

But there is no point in dividing for example 100/99, 100/98, 100/97. Numbers higher than half the prime number can never be a possibility in looking for a prime.

The way mine works is lets say you want to find out if 97 is a prime number. It can do it in 8 steps.

97/2 = 48.5 ok prime so far.

97/3 = 32.3 ok prime so far.

97/4 = 24.24 ok prime so far.

97/5 = 19.4 ok prime so far.

97/6 = 16.16 ok prime so far.

97/7 = 13.85 ok prime so far.

97/8 = 12.125 ok prime so far.

97/9 = 10.77 OK Its Prime!

Now I know there is not point in dividing by 4,6,8 because there even, but for some reason it would run slower if I took those out of the code. Go figure.

If the number to the right of the equals sign ever has a decimal its a possiblity the number we are looking for is prime. And we know to stop dividing when the right hand number is bigger than left hand denominator number. IE 10.77 is bigger than 9.

Relf

Are you saying my code is slower at finding higher prime numbers than other methods?

Are you saying my code is slower at finding higher prime numbers than other methods?

Unreal7000, it is slower.

First off, it is much better to use square root instead of ``halfPrime'' as was suggested by TheBeaver.

But even with such improvement it's still just a ``trial division'' that takes about sqrt(p) division operations (which are expensive) to test primality of p.

In opposite, Eratosthenes Prime Number Sieve doesn't divide at all, and it efficiently produces a whole bunch of primes in the specified interval as soon as there's enough memory to keep a sieve.

Miller-Rabin is a different story. Required number of operations is proportional to the _length_ of p, i.e. log(p).

You cannot feel the difference while p is small. But try to reach number like 10^9 and you'll see.

First off, it is much better to use square root instead of ``halfPrime'' as was suggested by TheBeaver.

But even with such improvement it's still just a ``trial division'' that takes about sqrt(p) division operations (which are expensive) to test primality of p.

In opposite, Eratosthenes Prime Number Sieve doesn't divide at all, and it efficiently produces a whole bunch of primes in the specified interval as soon as there's enough memory to keep a sieve.

Miller-Rabin is a different story. Required number of operations is proportional to the _length_ of p, i.e. log(p).

You cannot feel the difference while p is small. But try to reach number like 10^9 and you'll see.

Relf

Are you saying my code is slower at finding higher prime numbers than other methods?

Are you saying my code is slower at finding higher prime numbers than other methods?

relf:

crap; someone else thought of that first?

That's exactly the way my prime program worked. I guess it's not that bad considering the algorithm was good enough for Eratosthenes to be remembered by a few thousand years after he died. And i swear I'd never heard of Eratosthenes before.

crap; someone else thought of that first?

That's exactly the way my prime program worked. I guess it's not that bad considering the algorithm was good enough for Eratosthenes to be remembered by a few thousand years after he died. And i swear I'd never heard of Eratosthenes before.

Unreal7000,

Â The "best code" depends on the range you wish to find the numbers in. Also the "best code" code for testing if a single number is prime different for fing all primes in a range.

Â Here is a suggestion. Define a range eg say 0 to say 10^6 and specify if the code should find all primes in the range or test if a given single number in the range is prime (for the latter make the range larger eg 10^6 to 10^30 of even higher, there are lists of really long primes that will allow the code to be tested).

Now up the points to say 50 or a 100 and lets see who can write the best code and you can be the judge.

GwynforWeb

Â The "best code" depends on the range you wish to find the numbers in. Also the "best code" code for testing if a single number is prime different for fing all primes in a range.

Â Here is a suggestion. Define a range eg say 0 to say 10^6 and specify if the code should find all primes in the range or test if a given single number in the range is prime (for the latter make the range larger eg 10^6 to 10^30 of even higher, there are lists of really long primes that will allow the code to be tested).

Now up the points to say 50 or a 100 and lets see who can write the best code and you can be the judge.

GwynforWeb

For the kind comment of Unreal7000,

When a number is prime, you are considering several factors (literally)

a. it must be a natural number

b. it can be only be decomposed to exactly two factors: the number 1, and the number itself.

the heart of the algorithm i'm suggesting is to find a third factor of a number (since the number itself and 1 already count as the first two factors) if it exists. if it does, then it is not prime. you will know if the divisor is a factor if the quotient of a number and its divisor has no remainder.

supposing we take the number 12, which is composite.

12/11= 1 remainder 1, therefore 12%11=1 != 0, therefore 11 cannot be a factor.

12/10= 1 rem 2, 12%10=2 !=0, therefore not a factor

now, supposing we take the number 97.

12/9= 1 rem 3, 12%9=3 != 0, not a factor

.

.

.

but 12/6= 2 rem. 0, 12%6=0, therefore the number 12 is evenly divisible by 6. that makes 6 a factor of 12, making composite. now we know that it's definitely not prime.

97/96= 1 remainder 1. therefore, 97%96=1, != 0, therefore, 96 cannot be a factor of 97.

if up to

97/2= x remainder y, y!=0 for all permissible integral values of the divisor, then 97 does not have any other factors which can decompose it, and we therefore conclude that it is prime.

Regards,

Centeris

When a number is prime, you are considering several factors (literally)

a. it must be a natural number

b. it can be only be decomposed to exactly two factors: the number 1, and the number itself.

the heart of the algorithm i'm suggesting is to find a third factor of a number (since the number itself and 1 already count as the first two factors) if it exists. if it does, then it is not prime. you will know if the divisor is a factor if the quotient of a number and its divisor has no remainder.

supposing we take the number 12, which is composite.

12/11= 1 remainder 1, therefore 12%11=1 != 0, therefore 11 cannot be a factor.

12/10= 1 rem 2, 12%10=2 !=0, therefore not a factor

now, supposing we take the number 97.

12/9= 1 rem 3, 12%9=3 != 0, not a factor

.

.

.

but 12/6= 2 rem. 0, 12%6=0, therefore the number 12 is evenly divisible by 6. that makes 6 a factor of 12, making composite. now we know that it's definitely not prime.

97/96= 1 remainder 1. therefore, 97%96=1, != 0, therefore, 96 cannot be a factor of 97.

if up to

97/2= x remainder y, y!=0 for all permissible integral values of the divisor, then 97 does not have any other factors which can decompose it, and we therefore conclude that it is prime.

Regards,

Centeris

and 2 is the smallest prime number.

the smallest counting number is 1. 1 is neither composite nor prime: it only has one factor: which is 1. so the next smaller number is 2: which can only be evenly divisible by itself, and 1. that's why it is called a prime number.

and before starting the program,

initiate a preliminary check such that if

x%2=0 then it is even and it doesn't have to go through the "factor search" algorithm.

:)

Regards,

Centeris

the smallest counting number is 1. 1 is neither composite nor prime: it only has one factor: which is 1. so the next smaller number is 2: which can only be evenly divisible by itself, and 1. that's why it is called a prime number.

and before starting the program,

initiate a preliminary check such that if

x%2=0 then it is even and it doesn't have to go through the "factor search" algorithm.

:)

Regards,

Centeris

Centeris

But your algorithm goes from 97/96 all the way to 97/2.

Wouldn't it be faster to do it the way I have above.

It's only 8 steps.

But your algorithm goes from 97/96 all the way to 97/2.

Wouldn't it be faster to do it the way I have above.

It's only 8 steps.

Relf: "First off, it is much better to use square root instead of ``halfPrime'' as was suggested by TheBeaver..."

May I just ask you where I mentioned halfPrime????

I WAS THE ONE TALKING ABOUT SQUARE ROOTS FIRST IN THE FIRST RESPONSE TO THIS QUESTION!!!

ASS!!!

Anyways, for a 20 point question, the code is in and you guys are just more or less repeating/regurgitating what I submitted for an answer a week ago.

May I just ask you where I mentioned halfPrime????

I WAS THE ONE TALKING ABOUT SQUARE ROOTS FIRST IN THE FIRST RESPONSE TO THIS QUESTION!!!

ASS!!!

Anyways, for a 20 point question, the code is in and you guys are just more or less repeating/regurgitating what I submitted for an answer a week ago.

TheBeaver, calm down!

You didn't explicitly mention ``halfPrime'', but if you compare your program with original Unreal7000's one, they are doing essentially the same thing: trial division. The difference is in the range of divisors. While Unreal7000's program tries all divisors up to half of the number being tested (where the name ``halfPrime'' come from), your program tries all divisors up to square root.

"...you guys are just more or less repeating/regurgitating what I submitted..."

That's not true. Except trial division, there were at least two other methods of discovering primes: Eratosthenes Prime Number Sieve, and Miller-Rabin probabilistic primality test.

You didn't explicitly mention ``halfPrime'', but if you compare your program with original Unreal7000's one, they are doing essentially the same thing: trial division. The difference is in the range of divisors. While Unreal7000's program tries all divisors up to half of the number being tested (where the name ``halfPrime'' come from), your program tries all divisors up to square root.

"...you guys are just more or less repeating/regurgitating what I submitted..."

That's not true. Except trial division, there were at least two other methods of discovering primes: Eratosthenes Prime Number Sieve, and Miller-Rabin probabilistic primality test.

Relf: "First off, it is much better to use square root instead of ``halfPrime'' as was suggested by TheBeaver..."

Under review, I realise that this statement is one of those things that can be taken 2 different ways...

1) TheBeaver suggested using halfPrime, but its much better to use square root

2) As TheBeaver already stated, its better to use square root rather than halfPrime

...so "ASS" is retracted :)

...now back to the question of prime numbers...

I have given this a bit of thought and something just acurred to me. We find primes using only primes. My logic follows...

(1) We can eliminate all even numbers because they have 2 as a factor. This leaves us with odds.

(2) We can eliminate all numbers divisible by 3 (6,9,12,etc)

(3) Numbers divisible by 4 have been elliminated by (1)

(4) We can eliminate all numbers divisible by 5 (10, 15, 20,etc)

(5) Numbers divisible by 6 have been elliminated by (3)

...and so on. As you can see only the prime numbers need be tested as they elliminate all non-primes currently found. So the code for this would be (adapted from my earlier code)...

main()

{

unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

unsigned long count, primes_found=0, primes[10000]; // 10000 primes. Not enough? Then make it bigger. I'm explaining logic here

char prime;

while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

{

Â count = primes_found + 1; Â Â Â Â Â // Could count 0 to primes_found but this is a faster test

Â prime = 1;

Â while( count &&Â prime ) Â // at this point prime means "possibly prime" :)

Â {

Â Â count--;

Â Â prime = primes[count] % current; // Evenly divisible = not prime. D'OH!

Â }

Â if( prime )

Â {

Â Â printf("Prime number found ! %d\n", current);

Â Â primes[primes_found] = current;

Â Â primes_found++;

Â }

Â current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

}

}

This eliminates the need for squareroots, with the bonus of having the minimum of tests for each prime. ie The 16th prime is identified after only 15 factor tests!

Under review, I realise that this statement is one of those things that can be taken 2 different ways...

1) TheBeaver suggested using halfPrime, but its much better to use square root

2) As TheBeaver already stated, its better to use square root rather than halfPrime

...so "ASS" is retracted :)

...now back to the question of prime numbers...

I have given this a bit of thought and something just acurred to me. We find primes using only primes. My logic follows...

(1) We can eliminate all even numbers because they have 2 as a factor. This leaves us with odds.

(2) We can eliminate all numbers divisible by 3 (6,9,12,etc)

(3) Numbers divisible by 4 have been elliminated by (1)

(4) We can eliminate all numbers divisible by 5 (10, 15, 20,etc)

(5) Numbers divisible by 6 have been elliminated by (3)

...and so on. As you can see only the prime numbers need be tested as they elliminate all non-primes currently found. So the code for this would be (adapted from my earlier code)...

main()

{

unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

unsigned long count, primes_found=0, primes[10000]; // 10000 primes. Not enough? Then make it bigger. I'm explaining logic here

char prime;

while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

{

Â count = primes_found + 1; Â Â Â Â Â // Could count 0 to primes_found but this is a faster test

Â prime = 1;

Â while( count &&Â prime ) Â // at this point prime means "possibly prime" :)

Â {

Â Â count--;

Â Â prime = primes[count] % current; // Evenly divisible = not prime. D'OH!

Â }

Â if( prime )

Â {

Â Â printf("Prime number found ! %d\n", current);

Â Â primes[primes_found] = current;

Â Â primes_found++;

Â }

Â current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

}

}

This eliminates the need for squareroots, with the bonus of having the minimum of tests for each prime. ie The 16th prime is identified after only 15 factor tests!

OK. One change is that prime should be a long as well...

...

unsigned long prime;

...

...

unsigned long prime;

...

main()

{

unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

unsigned long prime, count, primes_found=0;

unsigned long primes[10000]; Â Â Â Â // 10000 primes. Not enough? Then make it bigger. I'm explaining logic here

while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

{

count = 0; Â Â Â Â Â

prime = 1;

while( (count <Â primes_found) &&Â prime ) Â // at this point prime means "possibly prime" :)

{

Â count++;

Â prime = primes[count] % current; // Evenly divisible = not prime. D'OH!

}

if( prime )

{

Â printf("Prime number found ! %d\n", current);

Â primes[primes_found] = current;

Â primes_found++;

}

current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

}

}

I've adjusted the code to test primes from the lowest to the highest due to the fact that the lower the potential factor, the higher the possiblity of a match. ie...

2 is a factor of 50% of numbers

3 is a factor of 33% of numbers

5 is a factor of 20% of numbers

and so on.

{

unsigned long current=3; Â Â Â Â Â Â Â // 1 and 2 cant be prime

unsigned long prime, count, primes_found=0;

unsigned long primes[10000]; Â Â Â Â // 10000 primes. Not enough? Then make it bigger. I'm explaining logic here

while( current - 1 ) Â Â Â Â Â Â Â Â Â // Once we loop back to 1, end

{

count = 0; Â Â Â Â Â

prime = 1;

while( (count <Â primes_found) &&Â prime ) Â // at this point prime means "possibly prime" :)

{

Â count++;

Â prime = primes[count] % current; // Evenly divisible = not prime. D'OH!

}

if( prime )

{

Â printf("Prime number found ! %d\n", current);

Â primes[primes_found] = current;

Â primes_found++;

}

current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

}

}

I've adjusted the code to test primes from the lowest to the highest due to the fact that the lower the potential factor, the higher the possiblity of a match. ie...

2 is a factor of 50% of numbers

3 is a factor of 33% of numbers

5 is a factor of 20% of numbers

and so on.

TheBeaver, proceeding in this direction, you'll come to Eratosthenes Prime Number Sieve. First, it sieves out all number divided by 2. Then all divided by 3, and so on.

I DONT USE HALFPRIME ANYMORE. I took it out and now it only goes to the square root of the possible prime number. This is what I have been trying to say. And it is a lot faster than my very first posting. Up to 37 times faster. Just read the revision history below.

// *************************************************

// File Name:prime_2.c Â

// Author: Samuel Dockery

// UserID: dockesa

// Date Created: September 28, 2002

// Date Revised: February 23, 2003 Â

// Purpose: Compute Prime Numbers and Display them to stdout. Â

//

// History:

// -Revision 2.01- February 24, 2003

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -currentPrime now increments by 2, since primes can't be divisible by 2.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Results in no noticable difference in computation.

// -Revision 2.00- February 23, 2003

// Â Â Â Â Â Â Â Â -Took halfPrime to the next step.

// Â Â Â Â Â Â Â Â Â Took halfPrime out, replaced by currentPrimeValue.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â halfPrime is half, so now currentPrimeVale is third, fourth, fifth, etc...

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Added primeNumberCounter but I have it commented out

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â because it slows the whole process down significantly.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Maybe I will work this in in Revision 3.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Benchmark to a value of 500,000 on AMD Athlon 1400Mhz:

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 1.10 Results: 3min. 6 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 2.00 Results: 0min. 5 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â WOW now Revision 2 is over 37 times faster than Revision 1.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â As numbers get higher I imagine this value increases (Maybe).

// -Revision 1.10- September 28, 2003

// Â Â Â Â Â Â Â Â -Added halfPrime to speed up processing.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Ex. 101/51, 101/52, is pointless so throw the rest out (53,54,55,etc...).

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -It increases speed but its less than 2 times faster than initial release.

// -Revision 1.00- September 28, 2003

// Â Â Â Â Â Â Â Â -Initial Release.

// *************************************************

// Include Files

#include <stdio.h> Â

// Global Constants

// None

// Global Types

// None

// Global Variables

// None

// Function Prototypes

// None

//****************************************************

// Function name:main Â

// Author: Samuel Dockery

// Purpose: Compute Prime Numbers and Display them to stdout.

// Â

// Input parameters: none Â

// Output parameters: none Â

// Returns:0 Â

// Error handling: none Â

// Record of changes: Revision 2

//****************************************************

Â

int main (void)

{

unsigned long currentPrimeDivision = 0;//Used to speed up processing.

unsigned long currentPrimeCounter = 1;//Tell what the current prime number count is.

unsigned long currentPrime = 3;//Current Prime number that is being calculating.

unsigned long smallNumber = 2;//Number used in checking to see if currentPrime is a prime number.

printf("2\n");//smallest prime number.

//printf(" Â Â Â Â Â Â Â Â Â Â Â Â 1\n");//Tells that this is first prime number.

for (;;)//Loop infinitly.

Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â if ( (currentPrime % smallNumber) != 0)//Checks to see if currentPrime is not evenly divisible.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = currentPrime / smallNumber;//Assign currentPrimeDivision a value.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â if (currentPrimeDivision <Â smallNumber)//If this is true it is a prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â printf("%li\n", currentPrime);//Print prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //currentPrimeCounter++;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //printf(" Â Â Â Â Â Â Â Â Â Â Â Â %li\n",currentPrimeCounter);//Displays what prime number this is.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime+=2;//Increment currentPrime by 2.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;//Reset currentPrimeDivision to 0 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â else

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber++;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â else//Do following statement becuase currentPrime is not Prime.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime+=2;//Increment currentPrime by 2.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }//End If

Â Â Â Â Â Â } Â Â Â Â Â

return (0);// Return 0 for successful processing. You probably won't get this far.

} // End main

// **************************

// File Name:prime_2.c Â

// Author: Samuel Dockery

// UserID: dockesa

// Date Created: September 28, 2002

// Date Revised: February 23, 2003 Â

// Purpose: Compute Prime Numbers and Display them to stdout. Â

//

// History:

// -Revision 2.01- February 24, 2003

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -currentPrime now increments by 2, since primes can't be divisible by 2.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Results in no noticable difference in computation.

// -Revision 2.00- February 23, 2003

// Â Â Â Â Â Â Â Â -Took halfPrime to the next step.

// Â Â Â Â Â Â Â Â Â Took halfPrime out, replaced by currentPrimeValue.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â halfPrime is half, so now currentPrimeVale is third, fourth, fifth, etc...

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Added primeNumberCounter but I have it commented out

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â because it slows the whole process down significantly.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Maybe I will work this in in Revision 3.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Benchmark to a value of 500,000 on AMD Athlon 1400Mhz:

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 1.10 Results: 3min. 6 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -Revision 2.00 Results: 0min. 5 sec.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â WOW now Revision 2 is over 37 times faster than Revision 1.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â As numbers get higher I imagine this value increases (Maybe).

// -Revision 1.10- September 28, 2003

// Â Â Â Â Â Â Â Â -Added halfPrime to speed up processing.

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Ex. 101/51, 101/52, is pointless so throw the rest out (53,54,55,etc...).

// Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â -It increases speed but its less than 2 times faster than initial release.

// -Revision 1.00- September 28, 2003

// Â Â Â Â Â Â Â Â -Initial Release.

// **************************

// Include Files

#include <stdio.h> Â

// Global Constants

// None

// Global Types

// None

// Global Variables

// None

// Function Prototypes

// None

//************************

// Function name:main Â

// Author: Samuel Dockery

// Purpose: Compute Prime Numbers and Display them to stdout.

// Â

// Input parameters: none Â

// Output parameters: none Â

// Returns:0 Â

// Error handling: none Â

// Record of changes: Revision 2

//************************

Â

int main (void)

{

unsigned long currentPrimeDivision = 0;//Used to speed up processing.

unsigned long currentPrimeCounter = 1;//Tell what the current prime number count is.

unsigned long currentPrime = 3;//Current Prime number that is being calculating.

unsigned long smallNumber = 2;//Number used in checking to see if currentPrime is a prime number.

printf("2\n");//smallest prime number.

//printf(" Â Â Â Â Â Â Â Â Â Â Â Â 1\n");//Tells that this is first prime number.

for (;;)//Loop infinitly.

Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â if ( (currentPrime % smallNumber) != 0)//Checks to see if currentPrime is not evenly divisible.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = currentPrime / smallNumber;//Assign currentPrimeDivision a value.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â if (currentPrimeDivision <Â smallNumber)//If this is true it is a prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â printf("%li\n", currentPrime);//Print prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //currentPrimeCounter++;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â //printf(" Â Â Â Â Â Â Â Â Â Â Â Â %li\n",currentPrimeCounter

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime+=2;//Incremen

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;//Reset currentPrimeDivision to 0 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â else

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber++;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }

Â Â Â Â Â Â Â Â Â Â Â else//Do following statement becuase currentPrime is not Prime.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â { Â Â Â Â Â

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrime+=2;//Incremen

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â smallNumber = 2;//Reset smallNumber to 2 to begin computing next prime number.

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â currentPrimeDivision = 0;

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â }//End If

Â Â Â Â Â Â } Â Â Â Â Â

return (0);// Return 0 for successful processing. You probably won't get this far.

} // End main

For the kind comment of Unreal,

I'm sorry, but no can do man. if you are to get the prime numbers you want, you MUST exhaust all other possible factors. please note that you only need ONE more distinct factor in order to disqualify a number from being prime. since a prime number so far cannot be defined by an expression, i think that's the best shot.

Regards,

Centeris

I'm sorry, but no can do man. if you are to get the prime numbers you want, you MUST exhaust all other possible factors. please note that you only need ONE more distinct factor in order to disqualify a number from being prime. since a prime number so far cannot be defined by an expression, i think that's the best shot.

Regards,

Centeris

you can use the sieve or arastothanes, which says that any number must have a prime divisor less than its square root. Hence, for any number, take its square root, then check all PRIMES less than that to see if they divide into the number. That should make your program a lot faster because there will be a lot less cases to check. This sieve is proven using the fundamental theorem of arithmetic, and I believe it should help.

for the kind comment of pedron:

i agree. but let's consider the number 9.

factors of 9= 1, 3, 9

-->1 is neither prime nor composite.

-->9 is composite

-->therefore, the only prime factor is 3.

->sqrt of 9 = 3

in this case, the prime factor is equal to the number's square root. so we may want to refine the sieve as

"any number must have a prime divisor less than OR EQUAL to its square root"

or unless, anybody else can find a counter example, then we'll refine it some more! :)

Regards,

Centeris

i agree. but let's consider the number 9.

factors of 9= 1, 3, 9

-->1 is neither prime nor composite.

-->9 is composite

-->therefore, the only prime factor is 3.

->sqrt of 9 = 3

in this case, the prime factor is equal to the number's square root. so we may want to refine the sieve as

"any number must have a prime divisor less than OR EQUAL to its square root"

or unless, anybody else can find a counter example, then we'll refine it some more! :)

Regards,

Centeris

Does anyone even know C programming becuase all this stuff ya'll talk about is in my code.

the code isn't commented well but i dont think it matters, i rote the code a while ago and it prompts the user and checks it, and just to let you know, any whole numbers less then the sqrt of a number that go into it make the number not prime

for example:

number being checked: 9

sqrt(9) = 3 u only have to check 3,2, and 1 if none of them go in, then its prime.

also, if u want, use fortran, a lot easier for things like this, also, try assembler, and you'll be happy with the speed

#include <stdio.h>

#define ONE_NUMBER Â 1

#define NUMBERS_ONLY 2

#define A_NUMBER Â Â 3

#define PRIME Â Â Â Â 4

#define NOT_PRIME Â Â 5

int isprime(int TheNumber) {

Â Â long double number;

Â Â int IsOdd;

Â Â int IsPrime;

Â Â long i;

Â Â long intnumber;

Â Â IsOdd = 0;

Â Â if ( TheNumber % 2 != 0 ) {

Â Â Â IsOdd = 1;

Â Â }

Â Â if ( IsOdd ) {

Â Â Â number = TheNumber / 2;

Â Â Â number = number - .5;

Â Â Â intnumber = number;

Â Â Â for ( i = 2; i <Â intnumber; i += 1 ) {

Â Â Â Â Â if ( TheNumber % i == 0 ) {

Â Â Â Â Â Â IsPrime = 0;

Â Â Â Â Â Â return(IsPrime);

Â Â Â Â Â }

Â Â Â }

Â Â } else {

Â Â Â IsPrime = 0;

Â Â }

Â Â return(IsPrime);

}

void HelpAndExit (int HelpKey, int rc, long TheNumber) {

Â Â switch(HelpKey) {

Â Â Â case 1: printf("You can only enter one number at a time.");

Â Â Â Â Â break;

Â Â Â case 2: printf("You must enter a number greater than or equal to 3.");

Â Â Â Â Â break;

Â Â Â case 3: printf("You must enter a number.");

Â Â Â Â Â break;

Â Â Â case 4: printf("%d, is a prime number.", TheNumber);

Â Â Â Â Â break;

Â Â Â case 5: printf("%d, is not a prime number.", TheNumber);

Â Â Â Â Â break;

Â Â }

Â Â printf("\n");

Â Â exit(rc);

}

int main(int argc, char *argv[]) {

Â Â int Â rc;

Â Â long Â TheNumber;

Â Â printf("Welcome to Prime Tester!\n");

Â Â printf("Processing input....\n");

/* printf("argc = %d\n", argc); */

Â Â if (argc == 2) {

Â Â Â TheNumber = atoi(argv[1]);

Â Â Â if (TheNumber <Â 3) {

Â Â Â Â Â HelpAndExit(NUMBERS_ONLY, 0, 0);

Â Â Â } else {

Â Â Â /* Given TheNumber to be checked for prime */

Â Â Â Â if (isprime(TheNumber)) {

Â Â Â Â Â Â HelpAndExit(PRIME, 1, TheNumber);

Â Â Â Â } else {

Â Â Â Â Â HelpAndExit(NOT_PRIME, 1, TheNumber);

Â Â Â Â }

Â Â Â }

Â Â } Â else {

Â Â Â if ( argc >Â 2 ) {

Â Â Â Â Â HelpAndExit(ONE_NUMBER, 0, 0);

Â Â Â } else {

Â Â Â Â HelpAndExit(A_NUMBER, 0, 0);

Â Â Â }

Â Â }

}

Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â

for the kind comment of nasch,

Â Â try removing the factor 1 when dividing in your pseudocode, because a number will always be divisible by 1. because if you don't, it will appear as if every number is composite. the smallest divisor must be 2.

Â Â try removing the factor 1 when dividing in your pseudocode, because a number will always be divisible by 1. because if you don't, it will appear as if every number is composite. the smallest divisor must be 2.

Â Alright, here we go...

Â Â There are three methods of generating primes that I know of: by division, by addition, and by multiplication.

Â Â 1) By division, the straightforward method.

Â Â You need an auxiliary supply of small prime numbers, up to the square root of the number, and you test N for divisibility by 2,3,5,7,11, etc. until P(i)**2 >= N. The small array need not be so large. There are close to n/Ln(n) primes less than n [that's n divided by the natural log (base e) of n], and n = Sqrt(N) but for big N you're in trouble. "Probabilistic" methods make me wince. There are also complex techniques and special case techniques as in Mersenne Primes.

Â Â For the generation of a collection of primes by division, obviously, candidate values for N should be odd, so that N:=N + 2; is the step. However, why not skip multiples of three also? You need steps of 2,4,2,4,etc, so inc:=6 - inc; N:=N + inc; Going beyond this is best done with a table of increments.

Â Â For testing the factors of N, you can stop one prime short of most people's choice of stopper: if f is the highest prime factore so far needed, and f = Prime(i), it will suffice for testing N not just until N >= f**2, but until N = Prime(i + 1)**2. For example, suppose f were 5. It would suffice until N = 49 = 7*7 because there would be no point in dividing by 7 earlier, because the quotient of that division would be less than 7, and that quotient is already being tested for. Thus, although 35 = 7*5, 35 need not be divided by 7 because it will have been tested by 5 already.

Â Â A programme implementing this is a little tricky because you have to get all the tricks to work "in phase" at the startup without the price of slowing the rest of the crunch.

Â Â The effort of this method is about N*Sqrt(N) since N numbers are to be tested, and each test involves effort proportional to Sqrt(N)

Â Â Herewith a Pascal prog...

Program PrimesByDivision; Uses DOS, crt;

{Concocted by R.N.McLean (whom God preserve), Victoria university, NZ.}

Â var N,P2: longint;

Â var hic,i,lpf,np: integer;

Â Const enuff = 8888; {3512 up to 32749, 6542 up to 65521.}

Â var Prime: array[0..enuff] of longint;

Â label nn;

Â BEGIN

Â Write('Primes: 2, 3, 5, ');

Â Prime[0]:=1; Prime[1]:=2; Â Â Â Â Â {These are inconvenient.}

Â Prime[2]:=3; Prime[3]:=5; Â Â Â Â Â {Ease starting arrangements.}

Â np:=3; Â Â Â Â Â Â Â Â Â Â Â {The number of primes so far.}

Â lpf:=2; Â Â Â Â Â Â Â Â Â Â Â {Last prime factor needed so far.}

Â p2:=Prime[lpf + 1]; Â Â Â Â Â {**2, thanks, but no.}

Â p2:=p2*p2; Â Â Â Â Â Â Â Â Â Â Â {No **2, sigh.}

Â n:=5; Â hic:=4; Â Â Â Â Â {Syncopation.}

nn:hic:=6 - hic; Â Â Â Â Â {Thus evade multiples of two and three.}

Â N:=N + hic; Â Â Â Â Â Â Â Â Â Â Â {The candidate.}

Â for i:=3 to lpf do if N mod Prime[i] = 0 then goto nn; Â Â Â Â Â {Next number please.}

Â if N >= P2 then Â Â Â Â Â {Is a further factor needed?}

Â Â begin Â Â Â Â Â Â Â Â Â Â Â {Yes. Actually, N = P2.}

Â Â lpf:=lpf + 1; Â Â Â Â Â {So we need a further prime factor.}

Â Â p2:=Prime[lpf + 1]; Â Â Â Â Â {**2, thanks. But no...}

Â Â p2:=p2*p2; Â Â Â Â Â Â Â Â Â Â Â {More timewasting.}

Â Â goto nn; Â Â Â Â Â Â Â Â Â Â Â {And look further.}

Â Â end; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {So much for finding further work.}

Â np:=np + 1; Â Â Â Â Â Â Â Â Â Â Â {We have another prime.}

Â write(N,', ');

Â if np <= enuff then Prime[np]:=N;

Â if n <Â 66666 then goto nn;

Â END.

Â Â 2) By addition: This is the sieve of Eratoshenes, and the addition aspect is that one has a monster bit array, sets all bits to one, then turn off every second bit beyond 2 (being a number divisible by two), then every third bit beyond 3, every fifth beyond 5, and so on.

Â Â Various tricks are possible, such as only representing the odd numbers for the bit array, and behaving differently as to whether the current step fits within a word's worth of bits (so more than one zap to that word from the array) or is larger.

Â Â Again a supply of small primes is required, again up to Sqrt(N) where N is the highest bit number, however these can be found by scanning the bit array as they are needed. With a bit of organisation, having sieved from 0:N, you can then re-use the bit array to sieve from N + 1 to 2N, and so on, but you will now need an auxiliary supply of primes, such as the disc file into which your programme is pouring the millions of primes.

Â Â This method is best for the mass production of a horde of primes that will be saved and referenced from then on. Thanks to going out for supper and the university building being locked on my return, I once printed all the primes up to 12,000,000 before the printer ran out of paper...

Â Â The effort of this scheme is proportional to the summation of N/Prime[i] for i = 1,2,... n

Â Â where n is such that Prime[n]**2 >Â N

Â Â since n passes need be made, each involving N/Prime[i] deeds.

Â Â The summation 1/i for i = 1...n differs from Ln(n) by Euler's constant, .577623 (and I don't recall which way) so the effort of this process, given that n is usually well above a hundred, is proportional to N/Ln(n)

Â Â Herewith another Pascal prog.

Program PrimesByAddition; Uses DOS, crt;

{Concocted by R.N.McLean (whom God preserve), Victoria university, NZ.}

Â var NBase,NLast,P: longint;

Â var i,l,np: integer;

Â Const enuff = 8888; lots = 2500;

Â var Prime: array[0..enuff] of longint;

Â var Sieve: array[0..lots] of boolean; {Eratoshenes lives on.}

Â label Next;

Â BEGIN

Â ClrScr;

Â Prime[0]:=1; Prime[1]:=2; Prime[2]:=3; np:=2;

Â for i:=0 to lots do Sieve[i]:=true; {Such blythe hope will soon be punctured.}

Â NBase:=3; Â Â Â Â Â Â Â Â Â Â {Corresponds to Sieve[0] and *must* be odd.}

Â NLast:=NBase + lots*2; Â Â Â {Only odd numbers are aligned with the sieve.}

Â l:=2; Â Â Â Â Â Â Â Â Â Â Â Â {So multiples of two are not reprresented.}

Â Repeat Â Â Â Â Â Â Â Â Â Â Â {Work my way through the primes.}

Â Â if l <= np then P:=Prime[l] {The next stepper.}

Â Â else Â Â Â Â Â Â Â Â Â Â Â {But my list may be inadequate.}

Â Â Â begin Â Â Â Â Â Â Â Â Â Â {So search for a further prime.}

Â Â Â i:=(Prime[np] - NBase) div 2; Â {Locate the last prime in my sieve.}

Â Â Â repeat i:=i + 1 until Sieve[i]; {Find the next survivor.}

Â Â Â P:=NBase + i*2; Â Â Â Â Â {Its corresponding number must be prime.}

Â Â Â np:=np + 1; Â Â Â Â Â Â Â {So I can augment my list.}

Â Â Â Prime[np]:=P; Â Â Â Â Â Â {With this.}

Â Â Â end; Â Â Â Â Â Â Â Â Â Â Â {And carry on.}

Â Â i:=(P - NBase) div 2 + P; Â {P is itself prime, while 2P is the first multiple not...}

Â {for i:=i to lots step p do Sieve[i]:=false; not available in pascal, sigh.}

Â Â while i <= lots do Â Â Â Â Â {And assuming i >= 0, as it is with NBase = 3.}

Â Â begin

Â Â Â Sieve[i]:=false; Â Â Â Â Â {Deny a non-prime.}

Â Â Â i:=i + p; Â Â Â Â Â Â Â Â {This is the *addition* of the method.}

Â Â end;

Â Â l:=l + 1; Â Â Â Â Â Â Â Â Â {Ready to use the next prime.}

Â Until P*P >Â NLast; Â Â Â Â Â {Remember that lesser primes *have* been used.}

Â Write ('Primes: 2, '); Â Â Â {Not seen by the sieve.}

Â for i:=0 to lots do if Sieve[i] then write(NBase + i*2,', ');

Â END.

Â Â 3) By multiplication... Again a bit array is set up, but notice that when the 5-step is in progress, it will zap bits such as 15 that have already been zapped by the 3-step. Indeed, every composite number will have its bit zapped by each of its prime factors. Aha! Why not contrive a scheme that generates all composite numbers? (since we can't find a scheme that itself generates all prime numbers) Once only. By the fundamental theorem of arithmetic, all numbers have onne prime factorisation.

Â Â As before, a small table of primes is needed for convenience. Clearly, given 2, all further powers of 2 need to be generated, and also, all products of 2 with 3, and 5, and so on; and for each power of two all products of that power with higher primes such as 3, 5, etc...

Â Â Then, given 3, etc.

Â Â I have not been able to concoct a satisfactory scheme, and have been pursuing other rabbits, but here a version of the key procedure:

Â Procedure Breed(n: longint; i: integer);

Â var m,p: longint;

Â Begin

Â Â if (n >Â lots) or (i >Â enuff) then exit;

Â Â p:=Prime(i);

Â Â m:=n*p; Â Â {n is the "base" number.}

Â Â while m <= lots do

Â Â begin

Â Â Â bit[m]:=true;

Â Â Â m:=m*p;

Â Â end;

Â Â m:=n;

Â Â while m <Â lots do

Â Â begin

Â Â Â Breed(m,i + 1); {Thus zap all products derived from my composites.}

Â Â Â m:=m*p;

Â Â end;

Â End;

Â Â Which is initiated by

Â Â for i:=1 to enuff do Breed(i,Prime[i]);

Â Â I can't decide whether the n of Breed ought be an already-zapped bit (from the previous level), or a possibly non-composite number (ie, a prime starter), or other indecisions.

Â Â It is a bit tricky. If you think to avoid generating the sequence m:=m*p; twice, you have to remember that the Breeds that are invoked will need higher primes, and if you haven't zapped all the powers of 2 first, the Get-the-next-prime from the bit array will find an as-yet unzapped composite number as its next prime number.

Â Â Conceptually simple, if there is already a supply of primes available whenever needed. But when the supply is derived from the same bit array as is in preparation, confusions arise. If an auxiliary array of primes is supplied, it must go up to N/2, not Sqrt(N)...

Â Â So, I have to stop, and think carefully...

Â Â Cheers,

Â Â Nicky.

Â Â There are three methods of generating primes that I know of: by division, by addition, and by multiplication.

Â Â 1) By division, the straightforward method.

Â Â You need an auxiliary supply of small prime numbers, up to the square root of the number, and you test N for divisibility by 2,3,5,7,11, etc. until P(i)**2 >= N. The small array need not be so large. There are close to n/Ln(n) primes less than n [that's n divided by the natural log (base e) of n], and n = Sqrt(N) but for big N you're in trouble. "Probabilistic" methods make me wince. There are also complex techniques and special case techniques as in Mersenne Primes.

Â Â For the generation of a collection of primes by division, obviously, candidate values for N should be odd, so that N:=N + 2; is the step. However, why not skip multiples of three also? You need steps of 2,4,2,4,etc, so inc:=6 - inc; N:=N + inc; Going beyond this is best done with a table of increments.

Â Â For testing the factors of N, you can stop one prime short of most people's choice of stopper: if f is the highest prime factore so far needed, and f = Prime(i), it will suffice for testing N not just until N >= f**2, but until N = Prime(i + 1)**2. For example, suppose f were 5. It would suffice until N = 49 = 7*7 because there would be no point in dividing by 7 earlier, because the quotient of that division would be less than 7, and that quotient is already being tested for. Thus, although 35 = 7*5, 35 need not be divided by 7 because it will have been tested by 5 already.

Â Â A programme implementing this is a little tricky because you have to get all the tricks to work "in phase" at the startup without the price of slowing the rest of the crunch.

Â Â The effort of this method is about N*Sqrt(N) since N numbers are to be tested, and each test involves effort proportional to Sqrt(N)

Â Â Herewith a Pascal prog...

Program PrimesByDivision; Uses DOS, crt;

{Concocted by R.N.McLean (whom God preserve), Victoria university, NZ.}

Â var N,P2: longint;

Â var hic,i,lpf,np: integer;

Â Const enuff = 8888; {3512 up to 32749, 6542 up to 65521.}

Â var Prime: array[0..enuff] of longint;

Â label nn;

Â BEGIN

Â Write('Primes: 2, 3, 5, ');

Â Prime[0]:=1; Prime[1]:=2; Â Â Â Â Â {These are inconvenient.}

Â Prime[2]:=3; Prime[3]:=5; Â Â Â Â Â {Ease starting arrangements.}

Â np:=3; Â Â Â Â Â Â Â Â Â Â Â {The number of primes so far.}

Â lpf:=2; Â Â Â Â Â Â Â Â Â Â Â {Last prime factor needed so far.}

Â p2:=Prime[lpf + 1]; Â Â Â Â Â {**2, thanks, but no.}

Â p2:=p2*p2; Â Â Â Â Â Â Â Â Â Â Â {No **2, sigh.}

Â n:=5; Â hic:=4; Â Â Â Â Â {Syncopation.}

nn:hic:=6 - hic; Â Â Â Â Â {Thus evade multiples of two and three.}

Â N:=N + hic; Â Â Â Â Â Â Â Â Â Â Â {The candidate.}

Â for i:=3 to lpf do if N mod Prime[i] = 0 then goto nn; Â Â Â Â Â {Next number please.}

Â if N >= P2 then Â Â Â Â Â {Is a further factor needed?}

Â Â begin Â Â Â Â Â Â Â Â Â Â Â {Yes. Actually, N = P2.}

Â Â lpf:=lpf + 1; Â Â Â Â Â {So we need a further prime factor.}

Â Â p2:=Prime[lpf + 1]; Â Â Â Â Â {**2, thanks. But no...}

Â Â p2:=p2*p2; Â Â Â Â Â Â Â Â Â Â Â {More timewasting.}

Â Â goto nn; Â Â Â Â Â Â Â Â Â Â Â {And look further.}

Â Â end; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {So much for finding further work.}

Â np:=np + 1; Â Â Â Â Â Â Â Â Â Â Â {We have another prime.}

Â write(N,', ');

Â if np <= enuff then Prime[np]:=N;

Â if n <Â 66666 then goto nn;

Â END.

Â Â 2) By addition: This is the sieve of Eratoshenes, and the addition aspect is that one has a monster bit array, sets all bits to one, then turn off every second bit beyond 2 (being a number divisible by two), then every third bit beyond 3, every fifth beyond 5, and so on.

Â Â Various tricks are possible, such as only representing the odd numbers for the bit array, and behaving differently as to whether the current step fits within a word's worth of bits (so more than one zap to that word from the array) or is larger.

Â Â Again a supply of small primes is required, again up to Sqrt(N) where N is the highest bit number, however these can be found by scanning the bit array as they are needed. With a bit of organisation, having sieved from 0:N, you can then re-use the bit array to sieve from N + 1 to 2N, and so on, but you will now need an auxiliary supply of primes, such as the disc file into which your programme is pouring the millions of primes.

Â Â This method is best for the mass production of a horde of primes that will be saved and referenced from then on. Thanks to going out for supper and the university building being locked on my return, I once printed all the primes up to 12,000,000 before the printer ran out of paper...

Â Â The effort of this scheme is proportional to the summation of N/Prime[i] for i = 1,2,... n

Â Â where n is such that Prime[n]**2 >Â N

Â Â since n passes need be made, each involving N/Prime[i] deeds.

Â Â The summation 1/i for i = 1...n differs from Ln(n) by Euler's constant, .577623 (and I don't recall which way) so the effort of this process, given that n is usually well above a hundred, is proportional to N/Ln(n)

Â Â Herewith another Pascal prog.

Program PrimesByAddition; Uses DOS, crt;

{Concocted by R.N.McLean (whom God preserve), Victoria university, NZ.}

Â var NBase,NLast,P: longint;

Â var i,l,np: integer;

Â Const enuff = 8888; lots = 2500;

Â var Prime: array[0..enuff] of longint;

Â var Sieve: array[0..lots] of boolean; {Eratoshenes lives on.}

Â label Next;

Â BEGIN

Â ClrScr;

Â Prime[0]:=1; Prime[1]:=2; Prime[2]:=3; np:=2;

Â for i:=0 to lots do Sieve[i]:=true; {Such blythe hope will soon be punctured.}

Â NBase:=3; Â Â Â Â Â Â Â Â Â Â {Corresponds to Sieve[0] and *must* be odd.}

Â NLast:=NBase + lots*2; Â Â Â {Only odd numbers are aligned with the sieve.}

Â l:=2; Â Â Â Â Â Â Â Â Â Â Â Â {So multiples of two are not reprresented.}

Â Repeat Â Â Â Â Â Â Â Â Â Â Â {Work my way through the primes.}

Â Â if l <= np then P:=Prime[l] {The next stepper.}

Â Â else Â Â Â Â Â Â Â Â Â Â Â {But my list may be inadequate.}

Â Â Â begin Â Â Â Â Â Â Â Â Â Â {So search for a further prime.}

Â Â Â i:=(Prime[np] - NBase) div 2; Â {Locate the last prime in my sieve.}

Â Â Â repeat i:=i + 1 until Sieve[i]; {Find the next survivor.}

Â Â Â P:=NBase + i*2; Â Â Â Â Â {Its corresponding number must be prime.}

Â Â Â np:=np + 1; Â Â Â Â Â Â Â {So I can augment my list.}

Â Â Â Prime[np]:=P; Â Â Â Â Â Â {With this.}

Â Â Â end; Â Â Â Â Â Â Â Â Â Â Â {And carry on.}

Â Â i:=(P - NBase) div 2 + P; Â {P is itself prime, while 2P is the first multiple not...}

Â {for i:=i to lots step p do Sieve[i]:=false; not available in pascal, sigh.}

Â Â while i <= lots do Â Â Â Â Â {And assuming i >= 0, as it is with NBase = 3.}

Â Â begin

Â Â Â Sieve[i]:=false; Â Â Â Â Â {Deny a non-prime.}

Â Â Â i:=i + p; Â Â Â Â Â Â Â Â {This is the *addition* of the method.}

Â Â end;

Â Â l:=l + 1; Â Â Â Â Â Â Â Â Â {Ready to use the next prime.}

Â Until P*P >Â NLast; Â Â Â Â Â {Remember that lesser primes *have* been used.}

Â Write ('Primes: 2, '); Â Â Â {Not seen by the sieve.}

Â for i:=0 to lots do if Sieve[i] then write(NBase + i*2,', ');

Â END.

Â Â 3) By multiplication... Again a bit array is set up, but notice that when the 5-step is in progress, it will zap bits such as 15 that have already been zapped by the 3-step. Indeed, every composite number will have its bit zapped by each of its prime factors. Aha! Why not contrive a scheme that generates all composite numbers? (since we can't find a scheme that itself generates all prime numbers) Once only. By the fundamental theorem of arithmetic, all numbers have onne prime factorisation.

Â Â As before, a small table of primes is needed for convenience. Clearly, given 2, all further powers of 2 need to be generated, and also, all products of 2 with 3, and 5, and so on; and for each power of two all products of that power with higher primes such as 3, 5, etc...

Â Â Then, given 3, etc.

Â Â I have not been able to concoct a satisfactory scheme, and have been pursuing other rabbits, but here a version of the key procedure:

Â Procedure Breed(n: longint; i: integer);

Â var m,p: longint;

Â Begin

Â Â if (n >Â lots) or (i >Â enuff) then exit;

Â Â p:=Prime(i);

Â Â m:=n*p; Â Â {n is the "base" number.}

Â Â while m <= lots do

Â Â begin

Â Â Â bit[m]:=true;

Â Â Â m:=m*p;

Â Â end;

Â Â m:=n;

Â Â while m <Â lots do

Â Â begin

Â Â Â Breed(m,i + 1); {Thus zap all products derived from my composites.}

Â Â Â m:=m*p;

Â Â end;

Â End;

Â Â Which is initiated by

Â Â for i:=1 to enuff do Breed(i,Prime[i]);

Â Â I can't decide whether the n of Breed ought be an already-zapped bit (from the previous level), or a possibly non-composite number (ie, a prime starter), or other indecisions.

Â Â It is a bit tricky. If you think to avoid generating the sequence m:=m*p; twice, you have to remember that the Breeds that are invoked will need higher primes, and if you haven't zapped all the powers of 2 first, the Get-the-next-prime from the bit array will find an as-yet unzapped composite number as its next prime number.

Â Â Conceptually simple, if there is already a supply of primes available whenever needed. But when the supply is derived from the same bit array as is in preparation, confusions arise. If an auxiliary array of primes is supplied, it must go up to N/2, not Sqrt(N)...

Â Â So, I have to stop, and think carefully...

Â Â Cheers,

Â Â Nicky.

This is a different approach but seems very fast, you can increase the array size to whatever you want or the languge will allow. Â Its Visual Basic but will transcribe easily. Â Note I am of the old school that believes 1 is a prime number whereas in schools in the UK they now teach that primes begin with 2.

Private Sub Form_Load()

Dim lngP(2000) As Long

Dim int1 As Integer, int2 As Integer

For int1 = 1 To 2000

Â Â lngP(int1) = int1

Next int1

For int1 = 2 To 1999

Â Â For int2 = int1 + 1 To 2000

Â Â Â Â If lngP(int2) <>Â 0 Then

Â Â Â Â Â Â If lngP(int2) Mod int1 = 0 Then lngP(int2) = 0

Â Â Â Â End If

Â Â Next int2

Next int1

For int1 = 1 To 2000

Â Â If lngP(int1) <>Â 0 Then lstPrimes.AddItem Str(lngP(int1))

Next int1

End Sub

Private Sub Form_Load()

Dim lngP(2000) As Long

Dim int1 As Integer, int2 As Integer

For int1 = 1 To 2000

Â Â lngP(int1) = int1

Next int1

For int1 = 2 To 1999

Â Â For int2 = int1 + 1 To 2000

Â Â Â Â If lngP(int2) <>Â 0 Then

Â Â Â Â Â Â If lngP(int2) Mod int1 = 0 Then lngP(int2) = 0

Â Â Â Â End If

Â Â Next int2

Next int1

For int1 = 1 To 2000

Â Â If lngP(int1) <>Â 0 Then lstPrimes.AddItem Str(lngP(int1))

Next int1

End Sub

Â Well, I too recall arguments with a friend over the admission of one to the set of Primes, and we finally compromised on it being the zero'th prime. There are many theorems about primes, and admitting one makes for messy statements. For instance, in "Every number has an unique prime factorisation", allowing one to be included means that there are an arbitrary number of factorisations, insignificantly differing by only the number of ones admitted.

Â Â I have been meaning to try a plot of N/Ln(N) vs. the number of primes up to N, but haven't got around to it...

Â Â With regard to your programme, it is fast only because the upper bound is just 2,000. It is an order N**2 method, because the outer loop requires N goes, and the inner loop requires on average N/2 goes.

Â Â Why are you using an array of numbers? The array elements contain either zero, or the same value as their index, so to acquire the value of the array element is easy. That is, make the array booleans and if a value is true, use the index's value, while if false, use zero, or similar.

Â Â I may have misinterpreted your prog, so here it is in Pascal, one version using numbers, the second using booleans...

Â const trialspan = 20000;

Â Procedure Trialw(Limit: longint);

Â var lngP: array[1..trialspan] of word;

Â var int1,int2: word;

Â Begin

Â Â if limit >Â trialspan then exit;

Â Â For int1:=1 to Limit do lngP[int1]:=int1;

Â Â For int1:=2 To Limit - 1 do

Â Â For int2:= int1 + 1 To limit do

Â Â Â Â If lngP[int2] <>Â 0 Then

Â Â Â Â Â Â If lngP[int2] Mod int1 = 0 Then lngP[int2]:=0;

Â Â if show then

Â Â For int1:= 1 To Limit do

Â Â Â If lngP[int1] <>Â 0 Then Write(lngP[int1],', ');

Â End;

Â Procedure Trialb(Limit: longint);

Â var lngP: array[1..trialspan] of boolean;

Â var int1,int2: word;

Â Begin

Â Â if limit >Â trialspan then exit;

Â Â For int1:=1 to Limit do lngP[int1]:=true;

Â Â For int1:=2 To Limit - 1 do

Â Â For int2:= int1 + 1 To limit do

Â Â Â Â If lngP[int2] Then

Â Â Â Â Â Â If int2 Mod int1 = 0 Then lngP[int2]:=false;

Â Â if show then

Â Â For int1:= 1 To Limit do

Â Â Â If lngP[int1] Then Write(int1,', ');

Â End;

Â Â And for comparison, I concocted a test prog. that ran and timed these methods, a tedious business. The addition (sieve) and division methods were so fast that they were run a hundred times each, but your method once only so that as the cpu timer available yielded only hundredths of a second, the lst digit is always zero. I didn't have confidence in the Multiplication method as I haven't settled on a good version. Anyway, the results, on a 266Mc Pentium II:

Seconds taken by various methods for finding prime numbers.

Â Up To,Addition,Division, Â TestW, Â TestB.

Â 1000 Â Â 0.001 Â Â 0.001 Â Â 0.220 Â Â 0.160

Â 2000 Â Â 0.001 Â Â 0.003 Â Â 0.820 Â Â 0.720

Â 3000 Â Â 0.002 Â Â 0.005 Â Â 1.810 Â Â 1.590

Â 4000 Â Â 0.002 Â Â 0.007 Â Â 3.130 Â Â 2.850

Â 5000 Â Â 0.003 Â Â 0.009 Â Â 5.000 Â Â 4.450

Â 6000 Â Â 0.003 Â Â 0.012 Â Â 7.140 Â Â 6.420

Â 7000 Â Â 0.004 Â Â 0.014 Â Â 9.780 Â Â 8.730

Â 8000 Â Â 0.005 Â Â 0.018 Â 12.630 Â 11.320

Â 9000 Â Â 0.005 Â Â 0.020 Â 15.930 Â 14.390

Â 10000 Â Â 0.006 Â Â 0.023 Â 19.610 Â 17.570

Â 11000 Â Â 0.007 Â Â 0.026 Â 23.950 Â 21.420

Â 12000 Â Â 0.007 Â Â 0.029 Â 28.450 Â 25.490

Â 13000 Â Â 0.008 Â Â 0.033 Â 33.290 Â 29.760

Â 14000 Â Â 0.009 Â Â 0.036 Â 38.560 Â 34.550

Â 15000 Â Â 0.009 Â Â 0.039 Â 44.100 Â 39.930

Â 16000 Â Â 0.009 Â Â 0.042 Â 50.200 Â 45.090

Â 17000 Â Â 0.011 Â Â 0.046 Â 56.790 Â 51.300

Â 18000 Â Â 0.011 Â Â 0.050 Â 63.660 Â 57.450

Â 19000 Â Â 0.012 Â Â 0.053 Â 70.750 Â 63.870

Â 20000 Â Â 0.012 Â Â 0.057 Â 78.430 Â 70.850

Â Â The division method is less than an order N*Sqrt(N) rate method, whereas the Addition (sieve) method is faster still, as I mentioned earlier. It is clearly order N+ since N items are being considered, but the factor is a bit difficult. The effort is proportional to the number of passes, and the work done by each pass. For a pass with prime p, the effort for that pass is proportional to 1/p since fewer steps are made with Â larger primes.

Â Â The summation of the series 1/p(i) has been considered in a paper by G.W. Hardy as I recall, but I have no access to it just now.

Â Â It is certainly less than the summation of the series 1/i, and (sum 1:n) of 1/i approaches Ln(n)

Â Â Your prog. is certainly short, but it divides too many numbers by too many possible factors...

Â Â I have been meaning to try a plot of N/Ln(N) vs. the number of primes up to N, but haven't got around to it...

Â Â With regard to your programme, it is fast only because the upper bound is just 2,000. It is an order N**2 method, because the outer loop requires N goes, and the inner loop requires on average N/2 goes.

Â Â Why are you using an array of numbers? The array elements contain either zero, or the same value as their index, so to acquire the value of the array element is easy. That is, make the array booleans and if a value is true, use the index's value, while if false, use zero, or similar.

Â Â I may have misinterpreted your prog, so here it is in Pascal, one version using numbers, the second using booleans...

Â const trialspan = 20000;

Â Procedure Trialw(Limit: longint);

Â var lngP: array[1..trialspan] of word;

Â var int1,int2: word;

Â Begin

Â Â if limit >Â trialspan then exit;

Â Â For int1:=1 to Limit do lngP[int1]:=int1;

Â Â For int1:=2 To Limit - 1 do

Â Â For int2:= int1 + 1 To limit do

Â Â Â Â If lngP[int2] <>Â 0 Then

Â Â Â Â Â Â If lngP[int2] Mod int1 = 0 Then lngP[int2]:=0;

Â Â if show then

Â Â For int1:= 1 To Limit do

Â Â Â If lngP[int1] <>Â 0 Then Write(lngP[int1],', ');

Â End;

Â Procedure Trialb(Limit: longint);

Â var lngP: array[1..trialspan] of boolean;

Â var int1,int2: word;

Â Begin

Â Â if limit >Â trialspan then exit;

Â Â For int1:=1 to Limit do lngP[int1]:=true;

Â Â For int1:=2 To Limit - 1 do

Â Â For int2:= int1 + 1 To limit do

Â Â Â Â If lngP[int2] Then

Â Â Â Â Â Â If int2 Mod int1 = 0 Then lngP[int2]:=false;

Â Â if show then

Â Â For int1:= 1 To Limit do

Â Â Â If lngP[int1] Then Write(int1,', ');

Â End;

Â Â And for comparison, I concocted a test prog. that ran and timed these methods, a tedious business. The addition (sieve) and division methods were so fast that they were run a hundred times each, but your method once only so that as the cpu timer available yielded only hundredths of a second, the lst digit is always zero. I didn't have confidence in the Multiplication method as I haven't settled on a good version. Anyway, the results, on a 266Mc Pentium II:

Seconds taken by various methods for finding prime numbers.

Â Up To,Addition,Division, Â TestW, Â TestB.

Â 1000 Â Â 0.001 Â Â 0.001 Â Â 0.220 Â Â 0.160

Â 2000 Â Â 0.001 Â Â 0.003 Â Â 0.820 Â Â 0.720

Â 3000 Â Â 0.002 Â Â 0.005 Â Â 1.810 Â Â 1.590

Â 4000 Â Â 0.002 Â Â 0.007 Â Â 3.130 Â Â 2.850

Â 5000 Â Â 0.003 Â Â 0.009 Â Â 5.000 Â Â 4.450

Â 6000 Â Â 0.003 Â Â 0.012 Â Â 7.140 Â Â 6.420

Â 7000 Â Â 0.004 Â Â 0.014 Â Â 9.780 Â Â 8.730

Â 8000 Â Â 0.005 Â Â 0.018 Â 12.630 Â 11.320

Â 9000 Â Â 0.005 Â Â 0.020 Â 15.930 Â 14.390

Â 10000 Â Â 0.006 Â Â 0.023 Â 19.610 Â 17.570

Â 11000 Â Â 0.007 Â Â 0.026 Â 23.950 Â 21.420

Â 12000 Â Â 0.007 Â Â 0.029 Â 28.450 Â 25.490

Â 13000 Â Â 0.008 Â Â 0.033 Â 33.290 Â 29.760

Â 14000 Â Â 0.009 Â Â 0.036 Â 38.560 Â 34.550

Â 15000 Â Â 0.009 Â Â 0.039 Â 44.100 Â 39.930

Â 16000 Â Â 0.009 Â Â 0.042 Â 50.200 Â 45.090

Â 17000 Â Â 0.011 Â Â 0.046 Â 56.790 Â 51.300

Â 18000 Â Â 0.011 Â Â 0.050 Â 63.660 Â 57.450

Â 19000 Â Â 0.012 Â Â 0.053 Â 70.750 Â 63.870

Â 20000 Â Â 0.012 Â Â 0.057 Â 78.430 Â 70.850

Â Â The division method is less than an order N*Sqrt(N) rate method, whereas the Addition (sieve) method is faster still, as I mentioned earlier. It is clearly order N+ since N items are being considered, but the factor is a bit difficult. The effort is proportional to the number of passes, and the work done by each pass. For a pass with prime p, the effort for that pass is proportional to 1/p since fewer steps are made with Â larger primes.

Â Â The summation of the series 1/p(i) has been considered in a paper by G.W. Hardy as I recall, but I have no access to it just now.

Â Â It is certainly less than the summation of the series 1/i, and (sum 1:n) of 1/i approaches Ln(n)

Â Â Your prog. is certainly short, but it divides too many numbers by too many possible factors...

Visual Basic again. Â This is the fastest I have managed, it takes 5.5 seconds for 100,000 primes running on an Athlon XP 1700+ under Windows 2000 Pro.

Private Sub Form_Load()

Const Max = 100000 ' Number of primes to be found

Dim lngP(Max) As Long ' Array to store primes

Dim lngJ As Long, lngK As Long ' indices

Dim lngN As Long, lngY As Long ' value to be tested for prime

Dim sngT1 As Single, sngT2 As Single, sngT As Single ' Times

Dim bytNt As Byte ' Used like a boolean flag

sngT1 = Timer ' Store start time in seconds from midnight

lngP(1) = 1

lngP(2) = 2

lngP(3) = 3

lngJ = 3

lngN = 5

Do While lngJ <Â Max

Â ' Use square root to limit range of tests

Â ' if primes up to the square root of lngN

Â ' will not divide into lngN then no higher primes will

Â Â lngY = Sqr(lngN) ' Square root of number to be tested

Â Â bytNt = 0 ' Set to 1 if a divide is successful

Â Â For lngK = 3 To lngJ

Â Â Â Â If lngN Mod lngP(lngK) = 0 Then bytNt = 1

Â Â Â Â If (lngP(lngK) >= lngY) Or (bytNt = 1) Then lngK = lngJ

Â Â Next lngK

Â Â If bytNt = 0 Then

Â Â Â lngJ = lngJ + 1

Â Â Â lngP(lngJ) = lngN

Â Â End If

Â Â lngN = lngN + 2 ' Next number to be tested (must be odd)

Loop

sngT2 = Timer ' Store end time in seconds from midnight

sngT = sngT2 - sngT1 ' Time taken in seconds

' Now display the primes in List Box lstPrimes

For lngJ = 1 To Max

Â Â lstPrimes.AddItem Str(lngP(lngJ))

Next lngJ

' Now display the time taken in List Box lstPrimes

Â Â lstPrimes.AddItem ""

Â Â lstPrimes.AddItem "TimeTaken"

Â Â strWork = Str(Format(sngT, "##0.000")) &Â " secs"

Â Â lstPrimes.AddItem strWork

End Sub

Private Sub Form_Load()

Const Max = 100000 ' Number of primes to be found

Dim lngP(Max) As Long ' Array to store primes

Dim lngJ As Long, lngK As Long ' indices

Dim lngN As Long, lngY As Long ' value to be tested for prime

Dim sngT1 As Single, sngT2 As Single, sngT As Single ' Times

Dim bytNt As Byte ' Used like a boolean flag

sngT1 = Timer ' Store start time in seconds from midnight

lngP(1) = 1

lngP(2) = 2

lngP(3) = 3

lngJ = 3

lngN = 5

Do While lngJ <Â Max

Â ' Use square root to limit range of tests

Â ' if primes up to the square root of lngN

Â ' will not divide into lngN then no higher primes will

Â Â lngY = Sqr(lngN) ' Square root of number to be tested

Â Â bytNt = 0 ' Set to 1 if a divide is successful

Â Â For lngK = 3 To lngJ

Â Â Â Â If lngN Mod lngP(lngK) = 0 Then bytNt = 1

Â Â Â Â If (lngP(lngK) >= lngY) Or (bytNt = 1) Then lngK = lngJ

Â Â Next lngK

Â Â If bytNt = 0 Then

Â Â Â lngJ = lngJ + 1

Â Â Â lngP(lngJ) = lngN

Â Â End If

Â Â lngN = lngN + 2 ' Next number to be tested (must be odd)

Loop

sngT2 = Timer ' Store end time in seconds from midnight

sngT = sngT2 - sngT1 ' Time taken in seconds

' Now display the primes in List Box lstPrimes

For lngJ = 1 To Max

Â Â lstPrimes.AddItem Str(lngP(lngJ))

Next lngJ

' Now display the time taken in List Box lstPrimes

Â Â lstPrimes.AddItem ""

Â Â lstPrimes.AddItem "TimeTaken"

Â Â strWork = Str(Format(sngT, "##0.000")) &Â " secs"

Â Â lstPrimes.AddItem strWork

End Sub

Â Well, I transliterated your prog. to Pascal and ran another comparison. Since Turbo Pascal 7 won't allow me an array with 100,000 elements, I changed the prog. to write the primes to a disc file as they were found, and only store primes up to the highest likely to be needed. With 32-bit integers, the highest needed value is a 16-bit integer (be careful about signed or not), and note that Prime[3512]=32749, Prime[3513]=32771, Prime[6542]=65521, Prime[6543]=65537, though this is counting with Prime[0] = 1. I modified your programme to use the global array Prime as do the other routines, rather than have its own.

Â Â The 100,000'th (+-1) prime is 1,299,709 and I further modified your prog. to stop when N exceeds the given limit, not when Number of Primes exceeds the given limit so as to correspond with the scheme used by the other trial routines. The sieve for instance doesn't know the number of primes it will find until after all the work is done.

Â Â The comparison was as follows:

Seconds taken by various methods for finding prime numbers.

Â Â Up To,Addition,Division, Â TestX.

Â 1299709 Â Â 1.870 Â 14.830 Â 22.250

Â Â All three routines produced output files with the same content, except for quibbles over starting with one, or ending with one prime number more or less than another routine.

Â Â The Pascal version I used is as follows:

Â Procedure Trialx(Limit: longint);

Â {var lngP: array[1..100000] of longint; Â Â Â Â Â {Changed to use Prime. Array to store primes}

Â var lngJ,lngK: longint; Â Â Â Â Â {indices}

Â var lngN,lngY: Longint; Â Â Â Â Â {value to be tested for prime}

Â var bytNt: Byte; Â Â Â Â Â Â Â Â Â Â Â {Used like a boolean flag}

Â Begin

Â Â if dump then begin assign(outx,'PrimesX.txt'); rewrite(outx); end;

Â Â Prime[1]:=1;

Â Â Prime[2]:=2;

Â Â Prime[3]:=3;

Â Â if dump then write(outx,'?Primes: 1, 2, 3, ');

Â Â lngJ:=3;

Â Â lngN:=5;

Â Â While lngN <Â Limit do{Was lngJ: Adjusted to control via N, not Number of primes}

Â Â begin Â Â Â Â Â {Use square root to limit range of tests}

Â Â Â Â Â Â Â Â Â Â Â {if primes up to the square root of lngN}

Â Â Â Â Â Â Â Â Â Â Â { will not divide into lngN then no higher primes will}

Â Â Â lngY:=Round(Sqrt(lngN)); {Square root of number to be tested}

Â Â Â bytNt:=0; Â Â Â Â Â {Set to 1 if a divide is successful}

Â Â Â For lngK:=3 To lngJ do

Â Â Â begin

Â Â Â Â If lngN Mod Prime[lngK] = 0 Then bytNt:=1;

Â Â Â Â If (Prime[lngK] >= lngY) Or (bytNt = 1) Then lngK:=lngJ;

Â Â Â end;

Â Â Â If bytNt = 0 Then

Â Â Â begin

Â Â Â Â lngJ:= lngJ + 1;

Â Â Â Â if lngJ <= PrimeLast then Â Prime[lngJ]:=lngN;

Â Â Â Â if dump then write (outx,lngN,', ');

Â Â Â end;

Â Â Â lngN:=lngN + 2; Â Â Â Â Â {Next number to be tested (must be odd)}

Â Â end; Â Â Â Â Â {of the While.}

Â Â if dump then close(outx);

Â Â End;

Â Â Your method is the division method, and although it is dodging the even numbers, it is not dodging multiples of three as well. Its rule for determining the last factor to try involves computing a square root. This is grievous. Rather than "if sqrt(x) <Â y", use "if x <Â y**2". And this limit can be improved upon slightly, by using the next prime along as I have mentioned in previous posts, and is to be seen in the Pascal source for the division method I posted.

Â Â But the addition method (the sieve of Eratoshenes) is way ahead in speed. Even though I am running the test on a Gateway 2000 w98se 266Mc PII system, it ran about twice as fast as your division method on your Athlon 1700Mc system. I have not attempted to run with array bound checking suppressed, etc, or other games.

Â Â It should be clear that the sieve method is the fastest, though it is not so easy to get to work...

Â Â Since again, TP7 does not allow monster arrays, the sieve routine made do with a boolean array of only 0:30000 elements, and worked in batches. Since even numbers are not involved, the first batch was for 3 to 60003, the second for 60005 to 120005, and so on. Setting this up is a bit tricky, especially the requirement to extract a list of primes for use in the sieving process as the sieving process runs, and in test runs, the sieve length may not be so great as to provide all of the possibly needed primes in the first batch.

Â Â The previous source was for one run of the sieve, but the multi-batch extension is easily bungled, so, here is the Pascal source...

Â Const PrimeLast = 6542; {Prime[3512]=32749,Prime[3513]=32771,Prime[6542]=65521,Prime[6543]=65537.}

Â Const SieveLast = 30000;{Thus span 2*(Sievelast + 1)}

Â var Prime: array[0..PrimeLast] of word;

Â var Sieve: array[0..SieveLast] of boolean; {Eratoshenes lives on.}

Â Procedure Addition(Limit: longint); Â Â Â Â Â {The sieve of Eratoshenes.}

Â var NBase,NLast,SLast,P: longint; Â Â Â Â Â {Possibly in many pieces.}

Â var NSavedP: word; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {The number of primes discovered as the need arose.}

Â Procedure ExtractAnother; Â Â Â Â Â {Finds a needed prime in the current sieve.}

Â {Perpetrated by R.N.McLean (whom God preserve), N.Z.}

Â Â var i: longint; Â Â Â Â Â Â Â {Big numbers can be involved.}

Â Â Begin

Â Â i:=Prime[NSavedP] - NBase; Â Â Â Â Â {How stands my last saved prime against the sieve?}

Â Â if i <Â 0 then i:=0 else i:=i div 2 + 1; Â Â Â Â Â {The next candidate in the current sieve.}

Â Â while (i <= SLast) and (not sieve[i]) do inc(i); Â Â Â Â Â {Find the next survivor. Presume there is one!}

Â Â if i >Â SLast then p:=0 Â Â Â Â Â {Cause trouble if nothing was found and a P is needed.}

Â Â Â else Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {But otherwise, we have a survivor.}

Â Â Â begin Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {So add it to my list.}

Â Â Â P:=NBase + i*2; Â Â Â Â Â Â Â Â Â Â Â {Its corresponding number must be prime.}

Â Â Â Inc(NSavedP); Â Â Â Â Â Â Â Â Â Â Â {So I can augment my list.}

Â Â Â Prime[NSavedP]:=P; Â Â Â Â Â {With this.}

Â Â Â end;

Â Â End;

Â var i,l: word;

Â var j: longint;

Â label Next;

Â Begin

Â Â if dump then begin assign(outa,'PrimesA.txt'); rewrite(outa); end;

Â Â if limit mod 2 = 0 then dec(limit); Â Â Â Â Â {Only odd numbers are of interest.}

Â Â Prime[0]:=1; Prime[1]:=2; Prime[2]:=3; NSavedP:=2; Â Â Â Â Â {Simplify.}

Â Â NBase:=3; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Corresponds to Sieve[0] and *must* be odd.}

Â Â Repeat Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Crank up multiple sieves.}

Â Â NLast:=NBase + SieveLast*2; Â Â Â Â Â Â Â Â Â Â Â {Only odd numbers are aligned with the sieve.}

Â Â if NLast >Â Limit then NLast:=Limit; Â Â Â Â Â {Possibly a truncated last span.}

Â Â SLast:=(NLast - NBase) div 2; Â Â Â Â Â {NB! The sieve starts with element zero.}

Â Â for i:=0 to SLast do Sieve[i]:=true;{Such blythe hope will soon be punctured.}

Â Â l:=2; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Start with 3. Multiples of two are not represented.}

Â Â Repeat Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Work my way through the primes.}

Â Â Â if l <= NSavedP then P:=Prime[l]{Locate the next stepper.}

Â Â Â else ExtractAnother; Â Â Â Â Â {But my list may be inadequate.}

Â Â Â j:=((NBase - 1) div P) + 1;{P*j is the first multiple of P within the sieve span.}

Â Â Â if j <= 1 then j:=3; Â Â Â Â Â {Don't zap P itself! eg, if NBase = 3, P = 5.}

Â Â Â if j mod 2 = 0 then inc(j);{Advance the jump point past any even multiple.}

Â Â Â i:=(P*j - NBase) div 2; Â Â Â Â Â {Thus, the first to zap. Overflow if big P and big sieve.}

Â Â {for i:=i to SLast step p do Sieve[i]:=false; not available in pascal, sigh.}

Â Â Â while i <= SLast do Â Â Â Â Â {So I'm stuck with this drivel.}

Â Â Â begin

Â Â Â Â Sieve[i]:=false; Â Â Â Â Â Â Â Â Â Â Â {Deny a non-prime.}

Â Â Â Â i:=i + p; Â Â Â Â Â Â Â Â Â Â Â {This is the *addition* of the method.}

Â Â Â end;

Â Â Â inc(l); Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Ready to use the next prime.}

Â Â Â if l <= NSavedP then P:=Prime[l] else ExtractAnother;

Â Â Until P*P >Â NLast; Â Â Â Â Â {If needed.}

Â Â While (NSavedP <Â PrimeLast) and (P >Â 0) do ExtractAnother; Â Â Â Â Â {The first decent-sized sieve will fill.}

Â Â if dump then

Â Â Â begin

Â Â Â if NBase = 3 then Write(outa,'+Primes: 2, '); Â Â {Not seen by the sieve.}

Â Â Â for i:=0 to Slast do

Â Â Â Â if Sieve[i] then

Â Â Â Â write(outa,NBase + LongInt(i)*2,', ');

Â Â Â end;

Â Â NBase:=NLast + 2; Â Â Â Â Â Â Â Â Â Â Â {The first of the next span. Must be odd!}

Â Â until NBase >Â Limit; Â Â Â Â Â Â Â Â Â Â Â {Slog on until the assigned mark.}

Â Â if dump then close (outa); Â Â Â Â Â {Why must *I* do this?}

Â End;

Â Â Much to my annoyance, in TP7 a bit array is not packed but employs one bit per byte, even if you specify "Packed", presumably because the overhead of access is not worth the space saving for small arrays. This however is a big array, and it may be worth working with a size other than a byte and explicitly coding a bit packing scheme. This would then lead to distinction between stepping with P <Â 2*bitsize and P >Â 2*bitsize and so yet more variants could be tested. So time bleeds away.

Â Â Obviously, this method could be augmented with a count of the primes as they are produced (getting back to the original issue), and why not, of prime number pairs also. Similarly, the primes could be saved in a disc file (in binary), and possibly, a compaction scheme used to save disc space if a large collection is needed. Other programmes could then invoke a function to return the i'th prime without repeating all the work...

Â Â The 100,000'th (+-1) prime is 1,299,709 and I further modified your prog. to stop when N exceeds the given limit, not when Number of Primes exceeds the given limit so as to correspond with the scheme used by the other trial routines. The sieve for instance doesn't know the number of primes it will find until after all the work is done.

Â Â The comparison was as follows:

Seconds taken by various methods for finding prime numbers.

Â Â Up To,Addition,Division, Â TestX.

Â 1299709 Â Â 1.870 Â 14.830 Â 22.250

Â Â All three routines produced output files with the same content, except for quibbles over starting with one, or ending with one prime number more or less than another routine.

Â Â The Pascal version I used is as follows:

Â Procedure Trialx(Limit: longint);

Â {var lngP: array[1..100000] of longint; Â Â Â Â Â {Changed to use Prime. Array to store primes}

Â var lngJ,lngK: longint; Â Â Â Â Â {indices}

Â var lngN,lngY: Longint; Â Â Â Â Â {value to be tested for prime}

Â var bytNt: Byte; Â Â Â Â Â Â Â Â Â Â Â {Used like a boolean flag}

Â Begin

Â Â if dump then begin assign(outx,'PrimesX.txt')

Â Â Prime[1]:=1;

Â Â Prime[2]:=2;

Â Â Prime[3]:=3;

Â Â if dump then write(outx,'?Primes: 1, 2, 3, ');

Â Â lngJ:=3;

Â Â lngN:=5;

Â Â While lngN <Â Limit do{Was lngJ: Adjusted to control via N, not Number of primes}

Â Â begin Â Â Â Â Â {Use square root to limit range of tests}

Â Â Â Â Â Â Â Â Â Â Â {if primes up to the square root of lngN}

Â Â Â Â Â Â Â Â Â Â Â { will not divide into lngN then no higher primes will}

Â Â Â lngY:=Round(Sqrt(lngN)); {Square root of number to be tested}

Â Â Â bytNt:=0; Â Â Â Â Â {Set to 1 if a divide is successful}

Â Â Â For lngK:=3 To lngJ do

Â Â Â begin

Â Â Â Â If lngN Mod Prime[lngK] = 0 Then bytNt:=1;

Â Â Â Â If (Prime[lngK] >= lngY) Or (bytNt = 1) Then lngK:=lngJ;

Â Â Â end;

Â Â Â If bytNt = 0 Then

Â Â Â begin

Â Â Â Â lngJ:= lngJ + 1;

Â Â Â Â if lngJ <= PrimeLast then Â Prime[lngJ]:=lngN;

Â Â Â Â if dump then write (outx,lngN,', ');

Â Â Â end;

Â Â Â lngN:=lngN + 2; Â Â Â Â Â {Next number to be tested (must be odd)}

Â Â end; Â Â Â Â Â {of the While.}

Â Â if dump then close(outx);

Â Â End;

Â Â Your method is the division method, and although it is dodging the even numbers, it is not dodging multiples of three as well. Its rule for determining the last factor to try involves computing a square root. This is grievous. Rather than "if sqrt(x) <Â y", use "if x <Â y**2". And this limit can be improved upon slightly, by using the next prime along as I have mentioned in previous posts, and is to be seen in the Pascal source for the division method I posted.

Â Â But the addition method (the sieve of Eratoshenes) is way ahead in speed. Even though I am running the test on a Gateway 2000 w98se 266Mc PII system, it ran about twice as fast as your division method on your Athlon 1700Mc system. I have not attempted to run with array bound checking suppressed, etc, or other games.

Â Â It should be clear that the sieve method is the fastest, though it is not so easy to get to work...

Â Â Since again, TP7 does not allow monster arrays, the sieve routine made do with a boolean array of only 0:30000 elements, and worked in batches. Since even numbers are not involved, the first batch was for 3 to 60003, the second for 60005 to 120005, and so on. Setting this up is a bit tricky, especially the requirement to extract a list of primes for use in the sieving process as the sieving process runs, and in test runs, the sieve length may not be so great as to provide all of the possibly needed primes in the first batch.

Â Â The previous source was for one run of the sieve, but the multi-batch extension is easily bungled, so, here is the Pascal source...

Â Const PrimeLast = 6542; {Prime[3512]=32749,Prime[3

Â Const SieveLast = 30000;{Thus span 2*(Sievelast + 1)}

Â var Prime: array[0..PrimeLast] of word;

Â var Sieve: array[0..SieveLast] of boolean; {Eratoshenes lives on.}

Â Procedure Addition(Limit: longint); Â Â Â Â Â {The sieve of Eratoshenes.}

Â var NBase,NLast,SLast,P: longint; Â Â Â Â Â {Possibly in many pieces.}

Â var NSavedP: word; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {The number of primes discovered as the need arose.}

Â Procedure ExtractAnother; Â Â Â Â Â {Finds a needed prime in the current sieve.}

Â {Perpetrated by R.N.McLean (whom God preserve), N.Z.}

Â Â var i: longint; Â Â Â Â Â Â Â {Big numbers can be involved.}

Â Â Begin

Â Â i:=Prime[NSavedP] - NBase; Â Â Â Â Â {How stands my last saved prime against the sieve?}

Â Â if i <Â 0 then i:=0 else i:=i div 2 + 1; Â Â Â Â Â {The next candidate in the current sieve.}

Â Â while (i <= SLast) and (not sieve[i]) do inc(i); Â Â Â Â Â {Find the next survivor. Presume there is one!}

Â Â if i >Â SLast then p:=0 Â Â Â Â Â {Cause trouble if nothing was found and a P is needed.}

Â Â Â else Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {But otherwise, we have a survivor.}

Â Â Â begin Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {So add it to my list.}

Â Â Â P:=NBase + i*2; Â Â Â Â Â Â Â Â Â Â Â {Its corresponding number must be prime.}

Â Â Â Inc(NSavedP); Â Â Â Â Â Â Â Â Â Â Â {So I can augment my list.}

Â Â Â Prime[NSavedP]:=P; Â Â Â Â Â {With this.}

Â Â Â end;

Â Â End;

Â var i,l: word;

Â var j: longint;

Â label Next;

Â Begin

Â Â if dump then begin assign(outa,'PrimesA.txt')

Â Â if limit mod 2 = 0 then dec(limit); Â Â Â Â Â {Only odd numbers are of interest.}

Â Â Prime[0]:=1; Prime[1]:=2; Prime[2]:=3; NSavedP:=2; Â Â Â Â Â {Simplify.}

Â Â NBase:=3; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Corresponds to Sieve[0] and *must* be odd.}

Â Â Repeat Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Crank up multiple sieves.}

Â Â NLast:=NBase + SieveLast*2; Â Â Â Â Â Â Â Â Â Â Â {Only odd numbers are aligned with the sieve.}

Â Â if NLast >Â Limit then NLast:=Limit; Â Â Â Â Â {Possibly a truncated last span.}

Â Â SLast:=(NLast - NBase) div 2; Â Â Â Â Â {NB! The sieve starts with element zero.}

Â Â for i:=0 to SLast do Sieve[i]:=true;{Such blythe hope will soon be punctured.}

Â Â l:=2; Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Start with 3. Multiples of two are not represented.}

Â Â Repeat Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Work my way through the primes.}

Â Â Â if l <= NSavedP then P:=Prime[l]{Locate the next stepper.}

Â Â Â else ExtractAnother; Â Â Â Â Â {But my list may be inadequate.}

Â Â Â j:=((NBase - 1) div P) + 1;{P*j is the first multiple of P within the sieve span.}

Â Â Â if j <= 1 then j:=3; Â Â Â Â Â {Don't zap P itself! eg, if NBase = 3, P = 5.}

Â Â Â if j mod 2 = 0 then inc(j);{Advance the jump point past any even multiple.}

Â Â Â i:=(P*j - NBase) div 2; Â Â Â Â Â {Thus, the first to zap. Overflow if big P and big sieve.}

Â Â {for i:=i to SLast step p do Sieve[i]:=false; not available in pascal, sigh.}

Â Â Â while i <= SLast do Â Â Â Â Â {So I'm stuck with this drivel.}

Â Â Â begin

Â Â Â Â Sieve[i]:=false; Â Â Â Â Â Â Â Â Â Â Â {Deny a non-prime.}

Â Â Â Â i:=i + p; Â Â Â Â Â Â Â Â Â Â Â {This is the *addition* of the method.}

Â Â Â end;

Â Â Â inc(l); Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â Â {Ready to use the next prime.}

Â Â Â if l <= NSavedP then P:=Prime[l] else ExtractAnother;

Â Â Until P*P >Â NLast; Â Â Â Â Â {If needed.}

Â Â While (NSavedP <Â PrimeLast) and (P >Â 0) do ExtractAnother; Â Â Â Â Â {The first decent-sized sieve will fill.}

Â Â if dump then

Â Â Â begin

Â Â Â if NBase = 3 then Write(outa,'+Primes: 2, '); Â Â {Not seen by the sieve.}

Â Â Â for i:=0 to Slast do

Â Â Â Â if Sieve[i] then

Â Â Â Â write(outa,NBase + LongInt(i)*2,', ');

Â Â Â end;

Â Â NBase:=NLast + 2; Â Â Â Â Â Â Â Â Â Â Â {The first of the next span. Must be odd!}

Â Â until NBase >Â Limit; Â Â Â Â Â Â Â Â Â Â Â {Slog on until the assigned mark.}

Â Â if dump then close (outa); Â Â Â Â Â {Why must *I* do this?}

Â End;

Â Â Much to my annoyance, in TP7 a bit array is not packed but employs one bit per byte, even if you specify "Packed", presumably because the overhead of access is not worth the space saving for small arrays. This however is a big array, and it may be worth working with a size other than a byte and explicitly coding a bit packing scheme. This would then lead to distinction between stepping with P <Â 2*bitsize and P >Â 2*bitsize and so yet more variants could be tested. So time bleeds away.

Â Â Obviously, this method could be augmented with a count of the primes as they are produced (getting back to the original issue), and why not, of prime number pairs also. Similarly, the primes could be saved in a disc file (in binary), and possibly, a compaction scheme used to save disc space if a large collection is needed. Other programmes could then invoke a function to return the i'th prime without repeating all the work...

RNMcLean,

You've clearly studied primes far more than I. Â I'd never heard of Eratoshenes but, having looked at your code, I remember reading about the method some years ago although I had forgotten it.

I only have Turbo Pascal 4 (from my DOS days). Â I copied your code into a .PAS file and tried it. TP 4 didn't like your unused label Next but having fixed that it ran fine.

I then transcribed it into VB and increased enuff and lots so that it found 100,000 primes. Â I compared the outputs of your program and mine and they matched.

Your "addition" method took 2.4 secs compared with my programs 5.5 seconds. Â Clearly a winner.

Re using a boolean array instead of an integer array. Â I think it is more elegant but has little impact on speed. Â VB like most languages holds booleans as an integer or byte with False being 0 and True being -1 so either way one is setting an integer (or a byte).

Finally, my grievous use of the square root. Â The form

x <Â y**2 is painfully slow compared with x <Â y*y but even using the latter makes the program slower when finding a large number of primes (I tried it). Â The reason is that the sqrt is outside the inner loop whereas x <Â y*y has to be inside the loop. Â (You can't win them all!).

You've clearly studied primes far more than I. Â I'd never heard of Eratoshenes but, having looked at your code, I remember reading about the method some years ago although I had forgotten it.

I only have Turbo Pascal 4 (from my DOS days). Â I copied your code into a .PAS file and tried it. TP 4 didn't like your unused label Next but having fixed that it ran fine.

I then transcribed it into VB and increased enuff and lots so that it found 100,000 primes. Â I compared the outputs of your program and mine and they matched.

Your "addition" method took 2.4 secs compared with my programs 5.5 seconds. Â Clearly a winner.

Re using a boolean array instead of an integer array. Â I think it is more elegant but has little impact on speed. Â VB like most languages holds booleans as an integer or byte with False being 0 and True being -1 so either way one is setting an integer (or a byte).

Finally, my grievous use of the square root. Â The form

x <Â y**2 is painfully slow compared with x <Â y*y but even using the latter makes the program slower when finding a large number of primes (I tried it). Â The reason is that the sqrt is outside the inner loop whereas x <Â y*y has to be inside the loop. Â (You can't win them all!).

Â Apologies about the unused label "Next"; clearly a vestige from an earlier version. At least it is less ugly than the digits-only labels of proper Pascal. I think I used the division method source file as a template, rather than retype the whole source file. TP7 is more lax. I shall fix it.

Â Â Likewise, I misread the square root. I was remembering the Mandelbrot calculation, where Sqrt(x**2 + y**2) >Â 2 is the key, and of course, this can most advantageously be x**2 + y**2 >Â 4 and *is* in the inner loop.

Â Â

Â Â However, the fact that your much faster computer runs a VB implementation at about the speed of this computer running compiled code for the sieve rather suggests that VB runs on your computer as interpreted code. If you wish, I could pass on a copy of TP7 for your messing about. It runs rather better under wunduhs than tp4. It is 2.27MB and thus beyond the Yahoo e-mail limit on attachments, but a scheme can be worked out. (RMcLean@yahoo.com)

Â Â One peculiarity of these computers is that more memory does not necessarily mean more speed. Expanding the sieve size to 50,000 led to a slightly slower run. I suspect that the cpu chip's internal memory sizes are involved in this. Thus, the addition method has been set up to deliver primes in batches (necessary anyway as no matter what limit, there are primes beyond it) and those batches are not as big as they could be even with TP7's limitations. I have not experimented to find an optimum size, as that would presumably be different on a different computer anyway.

Â Â When chasing larger primes, the step size will eventually be so large that only one zap per batch for that large P could be involved, because the prime is larger than the sieve span. There are so many twiddles to consider...

Â Â As for the sizes of Boolean variables, in Fortran, the language specification states that a version equivalent in size to the default integer is to be allowable. Thus, not a byte, but a 32-bit whack is involved. This sort of thing is important when devising complex record arrangements, but it is painful to waste so much. The values of boolean variables in their representation may be also 0 and 1, as well as other odd arrangements. PL/1 on the IBM mainframes used the high order bit of a word, whereas fortan used the low (or was it the other way about?) Attempting to use the value can lead to startlement, as did invoking a fortran routine from PL/1...

Â Â Again, I haven't experimented with the packing and unpacking of bits for the sieve on this machine. t was most worthwhile on a B6700 where the word size was 48 bits and the waste was painful. In another context it made a great improvement (a prog to play "hexapawn" on a NxM size board) because although the pack/unpack took time, the resulting compaction of disc file records enabled the disc buffering to hold many more records at once, and that made a great difference.

Â Â Likewise, I misread the square root. I was remembering the Mandelbrot calculation, where Sqrt(x**2 + y**2) >Â 2 is the key, and of course, this can most advantageously be x**2 + y**2 >Â 4 and *is* in the inner loop.

Â Â

Â Â However, the fact that your much faster computer runs a VB implementation at about the speed of this computer running compiled code for the sieve rather suggests that VB runs on your computer as interpreted code. If you wish, I could pass on a copy of TP7 for your messing about. It runs rather better under wunduhs than tp4. It is 2.27MB and thus beyond the Yahoo e-mail limit on attachments, but a scheme can be worked out. (RMcLean@yahoo.com)

Â Â One peculiarity of these computers is that more memory does not necessarily mean more speed. Expanding the sieve size to 50,000 led to a slightly slower run. I suspect that the cpu chip's internal memory sizes are involved in this. Thus, the addition method has been set up to deliver primes in batches (necessary anyway as no matter what limit, there are primes beyond it) and those batches are not as big as they could be even with TP7's limitations. I have not experimented to find an optimum size, as that would presumably be different on a different computer anyway.

Â Â When chasing larger primes, the step size will eventually be so large that only one zap per batch for that large P could be involved, because the prime is larger than the sieve span. There are so many twiddles to consider...

Â Â As for the sizes of Boolean variables, in Fortran, the language specification states that a version equivalent in size to the default integer is to be allowable. Thus, not a byte, but a 32-bit whack is involved. This sort of thing is important when devising complex record arrangements, but it is painful to waste so much. The values of boolean variables in their representation may be also 0 and 1, as well as other odd arrangements. PL/1 on the IBM mainframes used the high order bit of a word, whereas fortan used the low (or was it the other way about?) Attempting to use the value can lead to startlement, as did invoking a fortran routine from PL/1...

Â Â Again, I haven't experimented with the packing and unpacking of bits for the sieve on this machine. t was most worthwhile on a B6700 where the word size was 48 bits and the waste was painful. In another context it made a great improvement (a prog to play "hexapawn" on a NxM size board) because although the pack/unpack took time, the resulting compaction of disc file records enabled the disc buffering to hold many more records at once, and that made a great difference.

int is_prime(int n)

{

Â Â int r,i;

Â Â if ((n % 2) == 0) return 0;

Â Â r=(int)(sqrt(n));

Â Â for(i=3;i<=r;i+=2) if ((n % i) == 0) return 0;

Â Â return 1;

}

it's small, but i haven't gone through it looking for some way to make it look long and complicated

one of the assembly people here could probably optimize it to Â Â be thousands of clock cycles faster

~nasch

{

Â Â int r,i;

Â Â if ((n % 2) == 0) return 0;

Â Â r=(int)(sqrt(n));

Â Â for(i=3;i<=r;i+=2) if ((n % i) == 0) return 0;

Â Â return 1;

}

it's small, but i haven't gone through it looking for some way to make it look long and complicated

one of the assembly people here could probably optimize it to Â Â be thousands of clock cycles faster

~nasch

Â Not to be long and complicated about it, just because a programme is specified briefly does not mean that it will complete in a brief time. Or would you advocate a bubblesort over a Quicksort median-of-three with insertionsort for partitions less than 9, just because it is less complicated?

Â Â The sieve method can be written briefly too, but that brevity makes its use inconvenient. All sieve implementations will outrun division methods as N is increased, and by ever larger margins too.

Â Â In your implementation of the division method, instead of dividing by all odd numbers, you could skip divisions by multiples of three as well, by generating increments of 2,4,2,4,2,4,... rather than always 2 and this added complexity will enable faster running. This is easily enough achieved by i+=2 being replaced by i+=(inc:=6 - inc), if C handles this sort of trick (I'm more used to Algol), and yet further complexity might involve Â atable of increments so as to dodge multiples of five, and so on, though for lesser gain.

Â Â Going the other way, you have accepted the extra complexity of an explicit division by two test (which might better be an AND with 1) so as to avoid all further divisions by even numbers, presumably because the improvement in speed was worth it.

Â Â Implementing your version of the algorithm in assembler won't make all that much difference either. Â For this task, assembler advantages are to be found in access to both the quotient and the remainder after one division operation, which is routine in machine code but not available in higher-level languages, thus a C programme can't specify it.

Â Â The major improvement in the method is to divide only by prime numbers (from a handy array previously computed), and if you want to be careful, sqrt(n) might be computed on the low side of an integral value, as in sqrt(49) = 6.99999, so adding one, or some other such check, would be a precaution.

Â Â The sieve method can be written briefly too, but that brevity makes its use inconvenient. All sieve implementations will outrun division methods as N is increased, and by ever larger margins too.

Â Â In your implementation of the division method, instead of dividing by all odd numbers, you could skip divisions by multiples of three as well, by generating increments of 2,4,2,4,2,4,... rather than always 2 and this added complexity will enable faster running. This is easily enough achieved by i+=2 being replaced by i+=(inc:=6 - inc), if C handles this sort of trick (I'm more used to Algol), and yet further complexity might involve Â atable of increments so as to dodge multiples of five, and so on, though for lesser gain.

Â Â Going the other way, you have accepted the extra complexity of an explicit division by two test (which might better be an AND with 1) so as to avoid all further divisions by even numbers, presumably because the improvement in speed was worth it.

Â Â Implementing your version of the algorithm in assembler won't make all that much difference either. Â For this task, assembler advantages are to be found in access to both the quotient and the remainder after one division operation, which is routine in machine code but not available in higher-level languages, thus a C programme can't specify it.

Â Â The major improvement in the method is to divide only by prime numbers (from a handy array previously computed), and if you want to be careful, sqrt(n) might be computed on the low side of an integral value, as in sqrt(49) = 6.99999, so adding one, or some other such check, would be a precaution.

View this solution by signing up for a free trial.

Members can start a 7-Day free trial and enjoy unlimited access to the platform.

{

Â unsigned long current=1;

Â unsigned long factor, square;

Â char prime;

Â while( current ) Â Â Â Â Â Â Â Â Â // Once we loop back to 0, end

Â {

Â Â square = sqrt( current ) + 1; Â // function to get the square root (add one to allow for rounding)

Â Â factor = 1;

Â Â prime = 1;

Â Â while( (factor <Â square) &&Â prime ) Â // at this point prime means "possibly prime" :)

Â Â {

Â Â Â if( factor % current )

Â Â Â Â factor++;

Â Â Â else

Â Â Â Â prime = 0; Â Â Â Â Â Â Â Â Â Â Â // Evenly divisible = not prime. D'OH!

Â Â }

Â Â if( prime )

Â Â Â printf("Prime number found ! %d\n", current);

Â Â current += 2; Â Â Â Â Â Â Â Â Â // even numbers cant be prime

Â }

}

... if the sqrt fuction turns out to be too slow, you could get rid of square and instead of testing for (factor <Â square) you could do this ( current >Â (factor * factor)). But that is an extra calculation per factor per number