Fedora Linux Support Community & Resources Center
  #46  
Old Today, 04:46 AM
RupertPupkin Offline
Registered User
 
Join Date: Nov 2006
Location: Detroit
Posts: 6,483
linuxfedorafirefox
Re: Programming Challenge - list prime numbers

Here's an example in Tcl for listing the first N >= 1 primes. For this you'll need the Tcllib package from the Fedora repos (dnf install tcllib).

Save this script as getPrimes.tcl:
Code:
#!/usr/bin/tclsh
package require math::numtheory
puts -nonewline "2"
set N [lindex $argv 0]
set counter 1
set p 3
while {$counter < $N} {
   if {[math::numtheory::isprime $p] == 1} {
      puts -nonewline ", $p"
      incr counter
   }
   incr p 2
}
puts ""
Running the program to list the first 150 primes:
Code:
$ ./getPrimes.tcl 150
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73,
79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163,
 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443,
 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557,
 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647,
 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757,
 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863
__________________
OS: Fedora 25 x86_64 | Machine: HP Pavilion a6130n | CPU: AMD 64 X2 Dual-Core 5000+ 2.6GHz | RAM: 7GB PC5300 DDR2 | Disk: 400GB SATA | Video: ATI Radeon HD 4350 512MB | Sound: Realtek ALC888S | Ethernet: Realtek RTL8201N
Reply With Quote
  #47  
Old Today, 07:45 AM
ocratato Offline
Registered User
 
Join Date: Oct 2010
Location: Canberra
Posts: 2,482
linuxfirefox
Re: Programming Challenge - list prime numbers

Here is my latest version. It removes the even numbers from the sieve, eliminates a lot of modulo operations, and multi-threads the output. The downside is that you have to concatenate the output files, but there was nothing in the specifications saying the list had to be all in one file

Code:
// a multithreaded seive of Eratosthenes
//
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <bitset>
#include <thread>
#include <cmath>


using namespace std;
typedef unsigned long int ul64;

//
// Use a bitset to hold the flags indicating if we
// have found a prime. 
// The size should be small enough to fit in a CPU cache (32KB)
// given that the bitset stores 8 bits per byte.
// Each bit (at position i) represents if the number at 2*i+1 is a prime.
const ul64 BITSETSZ = 200000;
typedef bitset<BITSETSZ> Bits;

// Ideally this would be determined at run time.
const int NumThreads = 4;

// ------------------------------------------------

// mapping between integers and indexes
inline ul64 n2x( ul64 n ) { return (n-1)/2; }
inline ul64 x2n( ul64 x ) { return x*2+1; }

// Provide our own alternative to printf that is dedicated to
// printing integers and buffers its results to reduce I/O overhead.
//
inline void print_int(ul64 val, FILE *fp)
{
	// keep a buffer for each thread
	thread_local char outbuf[4096];
	thread_local int p = 0;

	char chars[16];
	int digits = 0;

	//
	// Pass zero to flush the buffer.
	if ( val == 0 )
	{
		fwrite(outbuf, 1, p, fp );
		p = 0;
		return;
	}

	do
	{
		chars[digits++] = ((val % 10) + 0x30);
		val /= 10;
	} while (val && digits < 16);


	while (digits>0)
	{
		outbuf[p++] = chars[--digits];
	}
	outbuf[p++] = '\n';
 
	if ( p > 4080 )
	{
		fwrite(outbuf, 1, p, fp );
		p = 0;
	}
}

// ------------------------------------------------

void print_page( Bits &p, ul64 pageNum, ul64 limit, FILE *fp )
{
	bool firstPage = (pageNum == 0);
	bool lastPage = ((pageNum+1)*BITSETSZ*2 >= limit);

	if( firstPage && lastPage )
	{
		print_int(2, fp);
		for( ul64 i=1, k; (k=x2n(i))<=limit; i++ )
		{
			if( p[i] ) print_int( k, fp );
		}
	}
	else if( firstPage )
	{
		print_int(2, fp);
		for(ul64 i=1; i<BITSETSZ; i++)
		{
			if( p[i] ) print_int( x2n(i), fp );
		}
	}
	else if( lastPage )
	{
		ul64 x = pageNum * BITSETSZ;
		for( ul64 i=0; x2n(x+i)<=limit; i++ )
		{
			if( p[i] ) print_int( x2n(x+i), fp );
		}
	}
	else
	{
		ul64 x = pageNum * BITSETSZ;
		for( ul64 i=0; i<BITSETSZ; i++ )
		{
			if( p[i] ) print_int( x2n(x+i), fp );
		}
	}
}

// ------------------------------------------------

void cull_pages( vector<Bits> & primes, ul64 limit, int thread, FILE *fp )
{
	//
	// cull out the multiples of each prime on each page
	//
	ul64 numPages = primes.size();
	Bits & primeSetZero = primes[0];

	//
	// determine range of pages for this thread.
	//
	ul64 R = (numPages - 1) / NumThreads;
	ul64 startPage = (thread - 1) * R + 1;
	ul64 endPage = startPage + R;
	if ( thread == NumThreads )
	{
		endPage = numPages;		
	}

	//
	// calculate the offsets for first page
	//
	vector<pair<ul64,ul64> > offsets;  // first is prime, second is offset
	for ( ul64 i=1; i<BITSETSZ; i++ )
	{
		ul64 P = startPage * (BITSETSZ * 2);
		if ( primeSetZero[i] )
		{
			ul64 j = x2n(i);
			ul64 k = P % j;
			k = P - k + j;
			if (!( k & 0x01) ) k += j;
			k = n2x(k-P);
			
			offsets.push_back(make_pair(j,k));
		}
	}

	//
	// process the pages for this thread
	//
	for( ul64 p=startPage; p<endPage; p++ )
	{
		Bits & primeSet = primes[p];
		for( auto &pr : offsets )
		{
			ul64 j = pr.first;
			ul64 k = pr.second;
			//
			// cull out the multiples
			//
			while( k < BITSETSZ )
			{
				primeSet[k] = false;
				k += j;
			}
			//
			// update offset for next page
			pr.second = k - BITSETSZ;
		}
		print_page( primeSet, p, limit, fp );
	}
	// flush the output buffer
	print_int(0, fp);
}

// ------------------------------------------------

int main(int argc,char *argv[])
{
        ul64 limit = atoll(argv[1]);
        printf("Calculating primes up to: %lu\n",limit);
	//
	// make sure we include the end in our seive
	limit++;

	//
	// The largest factor of limit is sqrt(limit)
	ul64 range = sqrt( limit );
	if ( BITSETSZ * 2 < range )
	{
		// 
		// Limit ourselves to 40 billion
		ul64 n = ul64(4) * BITSETSZ * BITSETSZ;
		printf( "max range for program is %lu\n", n );
		return 1;
	}

	//
	// Now we need enough bitsets to hold the entire range
	// These are just the odd numbers
	vector< Bits > primes( limit / BITSETSZ / 2 +1);
	//
	// initialise to set;
	for( auto &b : primes) b.set();

	//
	// first we calc the primes for the first page
	Bits & primeSetZero = primes[0];
	ul64 rng = n2x(sqrt(x2n(BITSETSZ)));
	for( ul64 i=1; i<=rng; i++)
	{
		if( primeSetZero[i] )
		{
			ul64 j = x2n(i);
			ul64 k = n2x(j*j);
			while( k < BITSETSZ )
			{
				primeSetZero[k] = false;
				k += j;
			}
		}
	}
	FILE *f[NumThreads+1];
	f[0] = fopen("primes7-0.dat", "w");
	print_page( primeSetZero, 0, limit, f[0] );
	print_int(0, f[0]);

	//
	// Then we cull the remaining pages.
	// start the threads
	std::thread threads[NumThreads];
	for (int t=1; t<=NumThreads; t++ )
	{
		char fname[20];
		sprintf(fname, "primes7-%d.dat", t );
		f[t] = fopen(fname, "w");
		// std::ref is our promise that primes will hang around
		// until the thread finishes.
		threads[t-1] = std::thread( cull_pages, std::ref(primes), limit, t, f[t] );
	}

	//
	// Wait until they are done
	for (int t=0; t<NumThreads; t++ )
	{
		threads[t].join();
	}

	//
	// tidy up the files
	//
	for( int i=0; i<=NumThreads; i++ )
	{
		fclose(f[i]);
	}

	return 0;
}
It is quite quick:
Code:
[brenton@acer primes]$ time ./primes7 1000000000
Calculating primes up to: 1000000000

real	0m2.918s
user	0m8.508s
sys	0m0.981s
__________________
Has anyone seriously considered that it might be turtles all the way down?
That's very old fashioned thinking.
The current model is that it's holographic nested virtualities of turtles, all the way down.
Reply With Quote
Reply

Tags
challenge, list, numbers, prime, programming

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

Forum Jump

Similar Threads
Thread Thread Starter Forum Replies Last Post
Programming Challenge: Create a list of X random numbers between Y and Z sea Programming & Packaging 6 14th December 2015 06:28 AM
A programming challenge. lsatenstein Programming & Packaging 3 4th January 2015 04:04 AM
Programming Challenge: Comparing Numbers in a String Flounder Programming & Packaging 23 10th January 2014 03:53 PM
[BASH]: Prime Numbers sea Guides & Solutions (Not For Questions) 14 21st November 2012 10:18 AM
Program to find prime numbers renkinjutsu Programming & Packaging 9 4th September 2011 07:40 AM


Current GMT-time: 14:30 (Monday, 27-03-2017)

TopSubscribe to XML RSS for all Threads in all ForumsFedoraForumDotOrg Archive
logo

All trademarks, and forum posts in this site are property of their respective owner(s).
FedoraForum.org is privately owned and is not directly sponsored by the Fedora Project or Red Hat, Inc.

Privacy Policy | Term of Use | Posting Guidelines | Archive | Contact Us | Founding Members

Powered by vBulletin® Copyright ©2000 - 2012, vBulletin Solutions, Inc.

FedoraForum is Powered by RedHat