Fedora Linux Support Community & Resources Center
  #46  
Old Yesterday, 04:46 AM
RupertPupkin Offline
Registered User
 
Join Date: Nov 2006
Location: Detroit
Posts: 6,484
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 Yesterday, 07:45 AM
ocratato Online
Registered User
 
Join Date: Oct 2010
Location: Canberra
Posts: 2,484
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
  #48  
Old Yesterday, 04:46 PM
lsatenstein Offline
Registered User
 
Join Date: Jun 2005
Location: Montreal, Que, Canada
Posts: 3,816
linuxfedorafirefox
Re: Programming Challenge - list prime numbers

I like APL. Here is a one line expression for the of Eratosthenes based on variable X

For the fun of it I assigned X 199

Code:
(2=+⌿0=(⍳X)∘.|⍳X)/⍳X←199
The Response

Code:
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
Here it is again for the primes under 100

Code:
(2=+⌿0=(⍳X)∘.|⍳X)/⍳X←100
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
The logic is as follows.
"Build a matrix X by X , Eliminate all rows which are multiples of two above 2"
then scan the matrix (across the columns), eliminating rows where a modulo zero occurred.
and VOILA.
__________________
Leslie in Montreal

Interesting web sites list
http://forums.fedoraforum.org/showth...40#post1697840
Reply With Quote
  #49  
Old Today, 02:53 AM
ocratato Online
Registered User
 
Join Date: Oct 2010
Location: Canberra
Posts: 2,484
linuxfirefox
Re: Programming Challenge - list prime numbers

Quote:
Originally Posted by lsatenstein View Post
I like APL. Here is a one line expression for the of Eratosthenes based on variable X


Code:
(2=+⌿0=(⍳X)∘.|⍳X)/⍳X←100
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
The logic is as follows.
"Build a matrix X by X , Eliminate all rows which are multiples of two above 2"
then scan the matrix (across the columns), eliminating rows where a modulo zero occurred.
and VOILA.
Very concise.
However, I suspect few of us know the APL symbols, so how it works is a bit too opaque for a thread that is aiming to illustrate various programming alternatives.
Secondly, you say that your solution uses an X by X matrix. Does this mean that it is O(N^2) in space and time requirements ?
__________________
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
  #50  
Old Today, 05:19 AM
RupertPupkin Offline
Registered User
 
Join Date: Nov 2006
Location: Detroit
Posts: 6,484
linuxfedorafirefox
Re: Programming Challenge - list prime numbers

Quote:
Originally Posted by ocratato View Post
how it works is a bit too opaque for a thread that is aiming to illustrate various programming alternatives
Eh, I don't see a problem with being "too opaque." Part of the fun of these threads is trying to figure out how the heck people did what they did. To be honest, I'm not sure that lsatanstein's APL code is all that much more "opaque" than your latest 200+ line C++ program.

Your second objection, about scalability, is the more legitimate one. lsatanstein only ran it for prime numbers less than 100, so he got a 100 by 100 matrix. I'd like to see how his program runs with 1 million as the upper limit. I suspect it might have problems.

For example, I am somewhat familiar with APL, so I was able to translate his program into its equivalent in the J programming language, which is an APL-ish clone/imitator. For primes up to 100 it works fine:
Code:
$ jconsole
   X =. 100
   (2=(X$1) +/ . * 0=(1+i.X)|/(1+i.X)) # (1+i.X)
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
But after a certain point it pretty much gives up:
Code:
   X =. 1000
   echo (2=(X$1) +/ . * 0=(1+i.X)|/(1+i.X)) # (1+i.X)
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 3...
That "3..." after 353 was put there by J. But it gets worse:
Code:
   X =. 1000000
   (2=(X$1) +/ . * 0=(1+i.X)|/(1+i.X)) # (1+i.X)
|out of memory
|   (2=(X$1)+/ .*0=(1+i.X)    |/(1+i.X))#(1+i.X)
So yeah, I suspect lsatanstein's program might have similar issues.

I think from now on an upper limit of 1 million should be the minimum that an example in this thread should be able to do.
__________________
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
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: 06:10 (Tuesday, 28-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