Prime Number Algorithm...

YT2095 · June 7, 2007

I wrote a prime number generator (sequential) about 20 years on a BBC micro.

naturally it`s in Basic, can anyone check it out on a more Modern system for me to make sure it`s a good prog?

I`ve managed to resurrect it from a 5.25 inch floppy and it`s working

but I don`t know if it was ever any good to start with!

10 FOR A%=1 TO 65535 STEP 2

20 FOR B%=2 TO A%/3

30 C = A%/B% : D%=C

40 IF D%=C GOTO 70

50 NEXT

60 PRINT A%;

70 NEXT A%

as you know the % forces Integer.

can this be improved upon in anyway? or made faster?

I chose 65535 as that`s a nice 16 bit number and wasn`t interested in going beyond that in primes, I use this prog and stop watch as a benchmark for CPU speed.

it`s not my Best "coding" as it was 20 years ago that I wrote it, and am Sure it can be improved upon!

any help?

timo · June 7, 2007

- You can increase the range. On modern computers, an integer is either 32 or 64 bit which is quite a lot more than 16 bit range. For even larger numbers, you can either check if the language you use supports types with a longer range, or write structures that do so yourself.

- You only need to check if a number N divides by any of the primes <N. There's no need to check if a number divides by 9 if it didn't divide by 3. You could keep a list of the primes already found and only check them. In other words, you could replace B% being all numbers from 2 to A%/3 to being all the primes you already found.

- You could try finding a function that gives you the current system time before start and after end => saves you using a stop watch.

bascule · June 7, 2007

I wrote one in Erlang last week It's an implementation of the Sieve of Eratosthenes

sieve(N) -> sieve(lists:seq(2, N), []).
sieve([], L) -> lists:reverse(L);
sieve([Prime|T], L) -> sieve([X || X <- T, X rem Prime /= 0], [Prime|L]).

Here's a quick English translation of each line:

sieve(N) -> sieve(lists:seq(2, N), []).

To find all of the primes up to N, start by generating a list of numbers between 2 and N (since 2 is the smallest prime larger than 1). Use this list of numbers to seed the algorithm.

sieve([Prime|T], L) -> sieve([X || X <- T, X rem Prime /= 0], [Prime|L]).

Knowing the number at the start of the list is a prime, go through the list and remove everything which is evenly divisible by the known prime. After that, store the prime in a list of known primes, then recursively call the same function, with the knowledge that the new number at the head of the input list is prime since it was not divisible by any number before it (besides 1).

sieve([], L) -> lists:reverse(L);

When the input list is empty, we've successfully removed all the non-primes and moved all the prime numbers into the primes list.

Dak · June 7, 2007

sieve(N) -> sieve(lists:seq(2, N), []).

To find all of the primes up to N, start by generating a list of numbers between 2 and N (since 2 is the smallest prime larger than 1). Use this list of numbers to seed the algorithm.

for sieve(N) -> sieve(lists:seq(2, N), []), assuming that '[]' is the incriment/step, and [] defaults to 1, wouldn't sieve(N) -> sieve(lists:seq(3, N), [2]) be faster by skipping all the even (and thus non-prime) numbers?

then just whack '2' at the begining of the generated list.

bascule · June 7, 2007

for sieve(N) -> sieve(lists:seq(2, N), []), assuming that '[]' is the incriment/step, and [] defaults to 1, wouldn't sieve(N) -> sieve(lists:seq(3, N), [2]) be faster by skipping all the even (and thus non-prime) numbers?

then just whack '2' at the begining of the generated list.

I could do that, but I'd have to start with a list of only odd numbers, which is a potential optimization.

To do that, the function could be rewritten as:

sieve(N) -> sieve([X * 2 - 1 || X <- lists:seq(2, N div 2)], [2]).

i.e. start with a list of numbers from 2 to N/2, then go through the list multiplying each number by 2 and then subtracting 1.

This will produce a list of all odd numbers between 3 and N.

Dak · June 7, 2007

that still seems like a lot of superfolouse maths just to eliminate even numbers... out of interest, can't you do the equivelent of YT's STEP 2 in erlang?

or something simple like:

oddIntegers=[3]				#a list
while oddIntegers[-1] <= n:		#till the last list entry has passed n 
oddIntegers += oddIntegers+2	#append the next odd number

is there some advantage to not using steps or something like the above, or is it because the above 'isn't functional programing' (not that i really understand what that is) so isn't part of erlang?

bascule · June 8, 2007

There's one very important aspect to Erlang which doesn't allow for the above: its state is immutable, also known as "single assignment"

Once you've assigned a variable, it cannot be changed. Variables are only variables because they can change from unbound to bound. Once bound, that's it cannot be unbound ever again (for that particular stack frame)

This approach is used in conjunction with optimized tail recursion. In Erlang and many other functional languages, if the last thing to be evaluated in that function is a call to itself, it works sort of like a "goto" statement, to where the old stack frame is cleared and a new stack frame created. This means a function can call itself indefinitely without eating up stack frames. This is the typical approach used for looping in most functional languages. It's also the approach used by lambda calculus.

It would be possible to write a small tail-recursive function to build a list of odd numbers

Aeternus · June 8, 2007

  sieve :: [integer] -> [integer]
  sieve [] = []
  sieve (h:t) = h : (sieve [x | x <- t, x `mod` h /= 0])

  primes = sieve [2..]

That's the sieve of eratosthenes in haskell You can just type primes and it'll keep generating primes as long as it doesn't run out of memory/stack space. It's not amazingly fast but it should do 65000 reasonably quickly.

Here is some rather old and extremely crappy code a few years ago when I was first learning bits of C++. The style is crap and it's horribly unreadable, and probably bug ridden but it works for the most part -

prime.cpp

// Primes.cpp : Defines the entry point for the console application.
//

// Include the include 
// Basic Headers
#include "stdafx.h"

// Basic Variable Declaration
unsigned long int num;
unsigned long int temp; // Value for use as highest prime (before init)
unsigned long int maxsize; 
bool found;
char numb[20];
char tchar[2];
float sqr;
long int pos;
short int add;
bool check;
char path[500];
unsigned long int fileplace;
bool dump; // Value to check whether an eop dump is needed 
long int test; // Remember later, for increase +1



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


// Ask for File
cout << "Please enter the path of the file to save to";
cin >> path;


// Set File Pointers
ofstream ofileptr(path,ios::app);
ifstream ifileptr(path,ios::in);


// Initialise Fileplace
fileplace  = 0;
temp = 0;


// Read in the Number to go up to
cout << "Please Enter a Number to go up to \n";
cin >> num;


// Read in the memory to use
cout << "Please enter the number of primes you would like to store in ram";
cin >> maxsize;


// Create new queue
queue tlist;




// Ensure file pointer is at the beginning of the file
ifileptr.seekg(0,ios::beg);

// Go through file and load into list (from beginning)
while((!ifileptr.eof()) && (tlist.size < maxsize))
{
	// Read in the number from file
	ifileptr >> temp;


	// Check whether list is full
	if(((tlist.size+1) == maxsize)  && (temp != 0))
	{
		// Add to list
		tlist.addToList(temp);
		// Find the place in file
		fileplace = ifileptr.tellg();

	}
	else if((tlist.size < maxsize) && (temp != 0) )
	{
		// Add to list
		tlist.addToList(temp);
	}

}


// Grab last prime
ifileptr.seekg(-50,ios::end);

while(!ifileptr.eof())
{
	ifileptr >> temp;

}


// Reset File Pointer
ifileptr.seekg(0,ios::beg);


// Check for file
if(temp == 0)
{

	// Add the first two prime numbers to the list
	// 
	tlist.addToList(2);
	tlist.addToList(3);

	// Set Temp to last known value
	temp = 3;

}







// Initialise Add
add = 1;

// Go up until the number
for( unsigned long int i = temp+1;i < num;i+=add)
{


// Initialisation

	// Initialise found
	found = false;

	// Reset File
	ifileptr.clear();

	ifileptr.seekg(0,ios::beg);

	// Reset File
	ifileptr.clear();



	// Get the max divisible number
	sqr = sqrt(i);
	// Intialise temp
	temp = 0;


	// Print out every 5000
	if(((i+1) % 5000)==0)
	{

		cout << i << "\n";
	}


	// Check for divisibility by 2 and then increment by two if necessary
	if((add!=2) && ((i%2)==1))
	{
		add = 2;
	}





// Memory


	// Reset List
	tlist.setToHead();


	// Go through all know primes and check
	while((!tlist.eol) && (tlist.crecord->prime <= (sqr+1)))
	{



		// Check if its not a prime
		if((i % tlist.crecord->prime) == 0) 
		{
			// Not a Prime
			found = true;
			break;



		};

		// Move on one
		tlist.moveOnOne();

	};












// FILE

	// Check for found 
	if((found != true) && (tlist.size >= maxsize))
	{
		// Set file pointer to the end of the list
		ifileptr.seekg(fileplace,ios::beg);

		// Check file
		// Check that it can still be prime
		// Check that the number is less than the square
		// Check that it is not the end of file
		while((found != true) & (temp<=(sqr+1))  & (!ifileptr.eof()))
		{   
			// Read in the lines

				// Read from File
				ifileptr >> temp;



			   if(!ifileptr.eof())
				{
					if((i % temp) == 0) 
					{
						// Found
						found = true;
						break;

					}



				}


		}



	}// End of Check if





// Check For Prime found


	// Check for prime
	if (!found)
	{
		// Check for need of pagefile
		if((tlist.size+1) == maxsize)
		{
			// Add to list
			tlist.addToList(i);

			// Initialise temp
			temp = 0;


			// Get File size
			ifileptr.seekg(0,ios::end);

			// Put to near end of file
			if(ifileptr.tellg() <100 )
			{
				// Go to beginning
				ifileptr.seekg(0,ios::beg);

			}
			else
			{
				// Go to near end
				ifileptr.seekg(100,ios::end);

			}

			// Go to near end of file
			ifileptr.seekg(-50,ios::end);

			// Go through file
			while(!(ifileptr.eof()))
			{
				ifileptr >> temp;
			}


			// Put put pointer to correct place
			ofileptr.seekp(0,ios::end);

			// Set temprecord to the head
			tlist.setToHead();

			// Go through list
			while(tlist.crecord != NULL)
			{
				// Check for the right place to start
				if(temp < tlist.crecord->prime)
				{
					// Write to file
					ofileptr << tlist.crecord->prime << "\n";
				}

				// Read in the next in the list
				tlist.moveOnOne();
			}


			// Find the place in file
			fileplace = ifileptr.tellg();
			// Set put pointer to the end
			ofileptr.seekp(0,ios::end);
		}
		else if(tlist.size >= maxsize)
		{
			// Dump to File



			// Convert from number to letter
			sprintf(numb,"%d",i);

			// Write the number to file
			ofileptr << numb << endl; 




		}
		else
		{
			// Add the prime to the list
			tlist.addToList(i);
		}



	};




};





// End of File Dump

// Check for dumped
if(!(tlist.size >= maxsize))
{
	// File didnt fill list 

	// Set File Pointer to beginning
	ofileptr.seekp(0,ios::beg);

	// Set record to head
	tlist.setToHead();


	if(tlist.crecord != NULL)
	{

		// Finished load all primes into file
		while(!tlist.eol)
		{
			sprintf(numb,"%d",tlist.crecord->prime);

			// Write the number to file
			ofileptr << numb << "\n"; 


			// Move list on one
			tlist.moveOnOne();


		};
	}

}





ifileptr.close();
ofileptr.close();

cout << "Thank You";




}

stdafx.h

// stdafx.h : include file for standard system include files,
//  or project specific include files that are used frequently, but
//      are changed infrequently
//

#if !defined(AFX_STDAFX_H__0475B6EC_188A_4437_B811_1E42BBB8D475__INCLUDED_)
#define AFX_STDAFX_H__0475B6EC_188A_4437_B811_1E42BBB8D475__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

#include <stdio.h>

#include <iostream>

using std::ios;

#include <fstream>
#include <math.h>



// Use the io class
using std::ofstream;
using std::ifstream;
using std::cout;
using std::cin;
using std::endl;




// The list
class queue {

// Type Declarations

struct listrecord;


//  The Pointer
typedef listrecord *listptr;



// A listrecord
struct listrecord {

	unsigned long int prime;
	listptr next;
};


// Variable Declarations

public:
// Vars
listptr head;  // Head
listptr tail; // Tail
long int size; // Current Size of List
listptr crecord; // Current Record
bool eol; // Boolean value for End Of List







// Function Declarations
public :

// Add To Linked list (at tail)
void addToList(long int number)
{
		// Prime found (Add to list)
		crecord = new(listrecord);
		crecord->prime = number;

		// Check there are primes
		if(tail != NULL)
		{
			tail->next = crecord;
		}
		else 
		{
			// No primes, add head etc
			head = crecord;
		}

		// Set the Tail
		tail = crecord;
		// Set the termination value
		tail->next = NULL;
		// Increment the size
		size++;



}

// Delete from List Declaration
void deleteFromList()
{
	// Move on one
	head = head->next;
	// Delete the node from memory
	delete(crecord);
	// Set The Record to the next record on
	crecord = head;
	// Decrease the list size
	size-= 1;

}


// Move the crecord on one
void moveOnOne()
{
		// Move on one
		crecord = crecord->next;
		// Check for EOL
		if(crecord == NULL)
		{
			// Set eol to true
			eol = true;
		}
		else
		{

			// Set eol to false
			eol = false;
		}
}


// Set crecord to head
void setToHead()
{

	// Set Crecord to head
	crecord = head;

	// Check for EOL
	if(crecord == NULL)
	{
		// Set eol to true
		eol = true;
	}
	else
	{
		// Set eol to false
		eol = false;
	}


}



public :	

// Class Constructor
queue()
{
	// Initialise
	head = NULL;
	tail = NULL;
	size = 0;
	eol = false;


}


};



// TODO: reference additional headers your program requires here

//{{AFX_INSERT_LOCATION}}

// Microsoft Visual C++ will insert additional declarations immediately before the previous line.

#endif // !defined(AFX_STDAFX_H__0475B6EC_188A_4437_B811_1E42BBB8D475__INCLUDED_)

As I said it sucks rather horribly and is a bit buggy but it can calculate the first million primes in 1 or 2 seconds. I'll probably get around to rewriting it in C or something at some point but there's not really much point as there are things that calculate primes far quicker. I think it's a bit buggy with the file handling (as well as other stuff) so if you do run it, it's probably best to tell it to store around the number of primes you want in memory (or say 1/5 or so of that number), otherwise it may stall. It didn't originally as I did use it to calculate all of them up to around a billion or so but it's crappy code and I didn't really know much about portability then and was writing it with dev c++ on windows. Oh well, memories

YT2095 · June 8, 2007

10 FOR A%=1 TO 65535 STEP 2

20 FOR B%=2 TO A%/3

30 C = A%/B% : D%=C

40 IF D%=C GOTO 70

50 NEXT

60 PRINT A%;

70 NEXT A%

was my original, and running on an overclocked BBC, I get to number 3169 in 15 mins.

now using line 20 as : 20 FOR B%=2 TO SQRA%

Im in the mid 12000`s in 8 minutes! :cool:

Niiiice!

Aeternus · June 8, 2007

  #include<stdio.h>
  #include<math.h>

  int main(int argc, char* args[]) {
      /* GO GO YT CODE! */
      int A, B;
      float C;
      int value = 1000000;
      if (argc > 1) {
          sscanf(args[1], "%d", &value);
      }
      for (A = 2; A <= value; A++) {
          C = sqrt((float)A);
          for(B=2; B <= C; B++) {
              if ((A % B) == 0) goto endofloop;
          }
          printf("%d\n", A);
          endofloop : ;
      }
  }

Here's a simple C adaptation of YTs code, without 1 as a prime

Sign In

Prime Number Algorithm...

Recommended Posts

YT2095

timo

bascule

Dak

bascule

Dak

bascule

Aeternus

YT2095

Aeternus

Create an account or sign in to comment

Create an account

Sign in

Browse

Activity

Important Information