



Programming & Packaging A place to discuss programming and packaging. 
Yesterday, 04:46 AM

Registered User


Join Date: Nov 2006
Location: Detroit
Posts: 6,484


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 DualCore 5000+ 2.6GHz  RAM: 7GB PC5300 DDR2  Disk: 400GB SATA  Video: ATI Radeon HD 4350 512MB  Sound: Realtek ALC888S  Ethernet: Realtek RTL8201N

Yesterday, 07:45 AM

Registered User


Join Date: Oct 2010
Location: Canberra
Posts: 2,484


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 multithreads 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 (n1)/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(kP);
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("primes70.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[t1] = 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.

Yesterday, 04:46 PM

Registered User


Join Date: Jun 2005
Location: Montreal, Que, Canada
Posts: 3,816


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.

Today, 02:53 AM

Registered User


Join Date: Oct 2010
Location: Canberra
Posts: 2,484


Re: Programming Challenge  list prime numbers
Quote:
Originally Posted by lsatenstein
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.

Today, 05:19 AM

Registered User


Join Date: Nov 2006
Location: Detroit
Posts: 6,484


Re: Programming Challenge  list prime numbers
Quote:
Originally Posted by ocratato
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 APLish 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 DualCore 5000+ 2.6GHz  RAM: 7GB PC5300 DDR2  Disk: 400GB SATA  Video: ATI Radeon HD 4350 512MB  Sound: Realtek ALC888S  Ethernet: Realtek RTL8201N

Thread Tools 
Search this Thread 


Display Modes 
Linear Mode

Posting Rules

You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
HTML code is Off



Current GMTtime: 06:10 (Tuesday, 28032017)







