Prime Number Sieves

Black Hole

Black Hole

May contain traces of nut
Black Hole said:
Did you ever use MS-DOS? The Linux command line isn't all that different.
Click to expand...
oijonesey said:
Yeah a little bit but most of the time only to change disks and directories and run games!! I have found some not too scary looking sites with Linux for Beginners I might start to go through. I think at face value it struck me that when OS's had developed into things you could operate with the tip of your finger suddenly if felt like Linux wanted to drag me right back to weekends with C&VG typing loads of lines of BASIC in!! BUT the Pi was all about teaching kids (and their Dads) that you can do more with a computer than get a high score so that's what I'm gonna learn meself!! My utopian daydream is that one day I may actually come up with something that ends up on the Hummy - rather than my current 30 odd year old personal best of
10 PRINT "FART"
20 GOTO 10
Click to expand...
That beats me - I once left the Physics department's PDP8 printing prime numbers on the teletype all weekend - but a roll of paper reading "FART!" all the way along would have been much more effective!
 
I remember writing something similar to oijonesey in a Dixons store in Manchester many, many years ago! Of course, I have 'matured' since then.....not. :)
 
It wasn't very big, only 6 digits as far as I can remember (I probably still have the paper roll somewhere!). It was a simple program written in good old BASIC (not the fancy Visual stuff) when a non-programmer could still understand programming, hand-typed on a TTY complete with line numbers. Something like (off the top of my head, and ignoring output formatting to print in multiple columns rather than just one):
Code: 
10 DIM A[1000]
20 LET P=3, A[1]=2, N=2
30 PRINT "2"
40 LET I=1
50 IF A[I]*A[I]>P GOTO 85
55 LET X=P/A[I]
60 IF X-INT(X)=0 GOTO 90
70 LET I=I+1
80 GOTO 50
85 PRINT P
86 LET A[N]=P
87 LET N=N+1
88 IF N>1000 GOTO 110
90 LET P=P+2
100 GOTO 40
110 END
I make no apology for the details not being quite right, not having written any BASIC since the '80s. Stick this in a modern computer and it would run to the end in a flash - with a PDP8 running interpreted BASIC it chuntered along like the Saturday football results.
 
I have rewritten it in 'C'. It could probably be cleaned up some more (to lose the goto) but it takes about 20ms to run on the Humax.
Code: 
#include <stdio.h>
#include <stdlib.h>

#define NUM_PRIMES 1000

int main()
{
unsigned a[NUM_PRIMES+1];
unsigned p=1;
unsigned n, i;
a[1]=2;
printf("2\n");
for (n=2; n<=NUM_PRIMES; n++) {
p2:  p+=2;
  i=1;
  while (1) {
    if (a[i]*a[i] > p) {
      printf("%d\n", p);
      a[n]=p;
      break;
      }
    if (p % a[i] == 0) goto p2;
    i++;
    }
  }
}
 
There's an error in my logic - the program breaks after 1000 primes, not the square of the 1000th prime as I intended.

I loaded a CBM Basic emulator onto my iPad and debugged the code (not in a very elegant way) and this runs:
Code: 
10 DIM A(1000)
20 LET P=3
21 LET A(1)=2
22 LET N=2
30 PRINT "2"
40 LET I=1
50 IF A(I)*A(I)>P GOTO 85
55 LET X=P/A(I)
60 IF X-INT(X)=0 GOTO 90
70 LET I=I+1
80 GOTO 50
85 PRINT P
86 LET A(N)=P
87 LET N=N+1
88 IF N>1000 GOTO 110
90 LET P=P+2
100 GOTO 40
110 LET I=1
120 IF A(I)*A(I)>P GOTO 180
130 LET X=P/A(I)
140 IF X-INT(X)=0 GOTO 190
150 LET I=I+1
160 GOTO 120
180 PRINT P
190 LET P=P+2
200 IF P<A(1000)*A(1000) GOTO 110
210 END
At the current rate I estimate it will take about 50 hours to complete (primes less than about 80,000,000).
 
9.5 hours, the Commodor is up to 10.5 million (I guess the PDP8 might have been limited by the speed of the teletype).

B*ll*x I just closed the app by mistake.
 
xyz321 said:
I have rewritten it in 'C'. It could probably be cleaned up some more (to lose the goto) but it takes about 20ms to run on the Humax.
Click to expand...
I have delved into my recollection of C (via a flowchart and a C language reference); can you tell me how to run this please?
Code: 
#include <stdio.h>
#include <stdlib.h>
#define NUM_PRIMES 1000
int main()
{
unsigned a[NUM_PRIMES+1];
unsigned p=3;
unsigned n=2;
unsigned i;
a[1]=2;
printf("2\n");
/* Seed the array of NUM_PRIMES prime numbers (2 is already done), outputting as we go */
do {
  for (i=1; (a[i]*a[i])<=p; i++) {
    if (p % a[i] == 0) break;
    }
  if (a[i]*a[i]>p) {
    printf("%d\n", p);
    a[n]=p;
    n++;
    }
  p+=2;
  } while (n <= NUM_PRIMES);
/* Optimised code to test for primes once the array is seeded, good up to the square of the largest seed */
do {
  for (i=2; (a[i]*a[i])<=p; i++) {
    if (p % a[i] == 0) break;
    }
  if (a[i]*a[i]>p) printf("%d\n", p);
  p+=2;
  } while (p < a[NUM_PRIMES]*a[NUM_PRIMES]);
}
 
If you were to call it 'primes.c':
Code: 
gcc primes.c -o primes
./primes
Note that the operator '^' is the bitwise xor operator which is probably not what you want. It is probably best to use a * a as in your original basic version. Otherwise you would need the 'pow(x, y)' function which is in the math library.
 
Thanks. I had imagined there was a raise to power operator which was more streamlined than a general purpose a*b. Source in post 8 adjusted.
 
Well what do you know. I had to install the gcc package, and correct a couple of typos from eliminating the "^", but it runs! It got up to 1,000,000 in seconds.
 
Yup. The execution might have been throttled by the output routine. I'm gonna try this:
Code: 
#include <stdio.h>
#include <stdlib.h>
#define NUM_PRIMES 1000
int main()
{
unsigned a[NUM_PRIMES+1];
unsigned p=3;
unsigned n=2;
unsigned i, c;
a[1]=2;
/* Seed the array of NUM_PRIMES prime numbers (2 is already done), outputting as we go */
do {
  for (i=1; (a[i]*a[i])<=p; i++) {
    if (p % a[i] == 0) break;
    }
  if (a[i]*a[i]>p) {
    a[n]=p;
    n++;
    }
  p+=2;
  } while (n <= NUM_PRIMES);
/* Optimised code to test for primes once the array is seeded, good up to the square of the largest seed */
do {
  for (i=2; (a[i]*a[i])<=p; i++) {
    if (p % a[i] == 0) break;
    }
  if (a[i]*a[i]>p) {
    c=p;
    n++;
    }
  p+=2;
  } while (p < a[NUM_PRIMES]*a[NUM_PRIMES]);
n--;
printf("%d\n", n);
printf("%d\n", c);
}
 
Apparently 62,710,559 is the 3,713,161th prime, and taking out the print statements saved 9 minutes.

I'm baffled by the printf formats though, I took the easy way out.
 
Instead of removing the printf statements, you could try:
Code: 
time ./primes > /dev/null
The first parameter in the printf statements is the output format. %d is integer output (it should be %u since the variables are declared as unsigned), '\n' is newline. If you wish to print n and c on the same line:
Code: 
printf("%u %u\n", n, c);
or
printf("%u ",n);
printf("%u\n",c);
To tidy up the output by allocating 10 characters to 'n':
Code: 
printf("%10u %u\n", n, c);
 
Black Hole said:
Apparently 62,710,559 is the 3,713,161th prime, and taking out the print statements saved 9 minutes.
Click to expand...

Update: I have confirmed 62,710,559 is a prime number, and according to the data from THIS WEB SITE (click) it is the 3,713,160th prime. n got post-incremented, so I should have subtracted 1 before printing. I have adjusted the source code in post 14.
 
Mike0001 said:
Prime Number Sieves are so pre-second century!
Click to expand...
In the late '80's (as opposed to the late '70's) I took to torturing an Amstrad 1640 with Mandelbrot Set images that required days to generate. Now they render in seconds to full iPad resolution.
 
You must log in or register to reply here.
Back
Top