Avatar of Unreal7000
Unreal7000
 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
Math / Science

Avatar of undefined
Last Comment
modulo

8/22/2022 - Mon
TheBeaver

main()
{
  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
TheBeaver

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
 }
}
integer

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.
Experts Exchange is like having an extremely knowledgeable team sitting and waiting for your call. Couldn't do my job half as well as I do without it!
James Murphy
GwynforWeb

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

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);
        }
    }
}
Unreal7000

ASKER
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
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
TheBeaver

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

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

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);
}
All of life is about relationships, and EE has made a viirtual community a real community. It lifts everyone's boat
William Peck
relf

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;
}
BigRat

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!
Unreal7000

ASKER
relf my code and yours seem to run at the same speed.
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
Unreal7000

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

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
Centeris

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
Experts Exchange has (a) saved my job multiple times, (b) saved me hours, days, and even weeks of work, and often (c) makes me look like a superhero! This place is MAGIC!
Walt Forbes
EVelasco

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

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

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!
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
relf

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;
}
Unreal7000

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

ASKER
Relf
Are you saying my code is slower at finding higher prime numbers than other methods?
Your help has saved me hundreds of hours of internet surfing.
fblack61
relf

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

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

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.
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
GwynforWeb

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
Centeris

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
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
I started with Experts Exchange in 2004 and it's been a mainstay of my professional computing life since. It helped me launch a career as a programmer / Oracle data analyst
William Peck
Unreal7000

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

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

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.
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
TheBeaver

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!
TheBeaver

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

...
unsigned long prime;
...
TheBeaver

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.
This is the best money I have ever spent. I cannot not tell you how many times these folks have saved my bacon. I learn so much from the contributors.
rwheeler23
relf

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

ASKER
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
Centeris

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
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
pedron

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

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
Unreal7000

ASKER
Does anyone even know C programming becuase all this stuff ya'll talk about is in my code.
Experts Exchange is like having an extremely knowledgeable team sitting and waiting for your call. Couldn't do my job half as well as I do without it!
James Murphy
nasch


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

                                                                                                           
Centeris

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

  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.
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
Lengore

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
RNMcLean

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

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
All of life is about relationships, and EE has made a viirtual community a real community. It lifts everyone's boat
William Peck
RNMcLean

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

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

RNMcLean

  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.
âš¡ FREE TRIAL OFFER
Try out a week of full access for free.
Find out why thousands trust the EE community with their toughest problems.
nasch

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
RNMcLean

  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.
ASKER CERTIFIED SOLUTION
modulo

THIS SOLUTION ONLY AVAILABLE TO MEMBERS.
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.
See Pricing Options
Start Free Trial
GET A PERSONALIZED SOLUTION
Ask your own question & get feedback from real experts
Find out why thousands trust the EE community with their toughest problems.