Prime numbers

A prime number is a number that has no divisors, apart from '1' and itself. There are lots of numbers and some of them are prime. If you want to hide things, like data, you need prime numbers. I need not hide things but I still want (to play with) prime numbers.

The first prime number is '1'. The second one is '2'. Two is the only even prime number. How about '3'? It is divisible by '1'. And by '3' but not by '2'. So three is prime as well. Four has three divisors: 1, 4 and 2. So '4' is not prime.

Century?

A century is what I call a logical series of 100 (decimal) consecutive numbers starting at xxx00 and ending at xxx99. I want to see how the primes are distributed across centuries. I didn't invent the word myself. I borrowed it from "snooker" where it is used in "century break".

Finding primes

Finding primes is easy. There is just one approach: brute force. Just try to divide a number by any other number and if it is not divisible, the number is prime. So if you want to see if 11 is prime, just divide it by 2, 3, 4, 5, 6, 7, 8, 9 and 10. But why would you try to divide by '9'? If it is divisible by 9 it is also divisible by '3'. The same applies for '2', '4', '6' and '8'. So to see if '11' is prime it is enough to test divide it by 2, 3, 5 and 7. That is: all prime numbers coming before 11.

First attempt: Prim01

Below is the first version of the primefinder. It finds all prime numbers below 1000. It starts out at '3' and steps up by '2' so it only tests the odd numbers.

The actual work is done in 'isPrime'. I have setup an array of ten thousand integers to store all detected prime numbers. isPrime divides the number to be tested against each prime in the array. it does so with the MOD operator. If MOD returns '0', the number has a divisor. When not, it is the next prime number and it will be stored in the array.

MODULE prim01;

IMPORT In, Out;

CONST	maxIndex	= 10000;

VAR	primes		: ARRAY maxIndex OF INTEGER;
	n		: INTEGER;


PROCEDURE isPrime (num : INTEGER) : BOOLEAN;

VAR	i	: INTEGER;

BEGIN
  i := 1;
  WHILE  primes [i] # 0  DO
    IF  num MOD primes [i] = 0  THEN  RETURN FALSE  END;
    INC (i)
  END;
  primes [i] := num;
  RETURN  TRUE
END isPrime;
			 

PROCEDURE Init;

VAR	i	: INTEGER;

BEGIN
  primes [0] := 1;
  primes [1] := 2;
  FOR  i := 2  TO  maxIndex - 1  DO  primes [i] := 0  END;
END Init;


BEGIN
  Init;
  FOR  n := 3  TO  1000  BY  2  DO
    Out.Int (n, 4);
    Out.String (" is ");
    IF  isPrime (n) = FALSE  THEN
      Out.String ("not ")
    END;
    Out.String ("prime.");
    Out.Ln
  END;
END prim01.
   
Compile the source prim01.mod with obc prim01.mod -o prim01

See it run:

   3 is prime.
   5 is prime.
   7 is prime.
   9 is not prime.
  11 is prime.
  13 is prime.
 ...   ...
 989 is not prime.
 991 is prime.
 993 is not prime.
 995 is not prime.
 997 is prime.
 999 is not prime.
   
For a first test it is nice. The output is a bit boring though.

Improving the output format : prim02.mod

I changed prim01 such that the output became more in association with the nature of prime numbers:

      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 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997
   
Now this is more suitable for showing part of an infinite series. As you can see, the values 1 and 2 are still missing.

Prim02 : the source

Compile the source prim02.mod with obc prim02.mod -o prim02

MODULE prim02;

IMPORT Out;

CONST	maxIndex	= 10000;

VAR	primes		: ARRAY maxIndex OF INTEGER;
	n		: INTEGER;


PROCEDURE isPrime (num : INTEGER) : BOOLEAN;

VAR	i	: INTEGER;

BEGIN
  i := 1;
  WHILE  primes [i] # 0  DO
    IF  num MOD primes [i] = 0  THEN  RETURN FALSE  END;
    INC (i)
  END;
  primes [i] := num;
  RETURN  TRUE
END isPrime;
			 

PROCEDURE Init;

VAR	i	: INTEGER;

BEGIN
  primes [0] := 1;
  primes [1] := 2;
  FOR  i := 2  TO  maxIndex - 1  DO  primes [i] := 0  END;
END Init;


BEGIN
  Init;
  FOR  n := 3  TO  1000  BY  2  DO
    IF  isPrime (n)  THEN
      Out.Int (n, 4)
    END
  END;
  Out.Ln
END prim02.
   
The main differences are in the MAIN routine which is now a lot simpler.

Add some statistics: Prim03.mod

Now would't it be nice to see in one glance how many primes were detected in a given domain? Yes of course. So I made prim03.mod. What it adds is:

For very big domains (above 1 - 400000) the 32 bit INTEGER holding the iterations count overflows.... So either don't go too far or don't pay attention to a negative count.

Prim03.mod : the source

MODULE prim03;

IMPORT	Out;

CONST	maxIndex	= 100000;
	maxCnt		= 16;

VAR	primes		: ARRAY maxIndex OF INTEGER;
	m, n, total, t	: INTEGER;


PROCEDURE isPrime (num : INTEGER) : BOOLEAN;

VAR	i	: INTEGER;

BEGIN
  i := 1;
  WHILE  primes [i] # 0  DO
    IF  num MOD primes [i] = 0  THEN  RETURN FALSE  END;
    INC (t);
    INC (i)
  END;
  primes [i] := num;
  INC (total);
  RETURN  TRUE
END isPrime;
			 

PROCEDURE Init;

VAR	i	: INTEGER;

BEGIN
  m := 2;
  t := 0;
  primes [0] := 1;
  primes [1] := 2;
  total := 2;
  FOR  i := 2  TO  maxIndex - 1  DO  primes [i] := 0  END;
  Out.String ("     1     2")  
END Init;


BEGIN
  Init;
  n := 3;
  WHILE  n < 5000  DO
    IF  isPrime (n)  THEN
      Out.Int (n, 8);
      INC (m);
      IF  m = maxCnt  THEN
        Out.Ln;
	m := 0
      END
    END;
    INC (n);
    IF  n MOD 100 = 0  THEN		(* Century break?	*)  
      Out.Ln;
      m := 0
    END
  END;
  Out.Ln;
  Out.String ("Total primes = ");	Out.Int (total, 6);	Out.Ln;
  Out.String ("Iterations   = ");	Out.Int (t, 9);		Out.Ln
END prim03.
   
As you can see, the main loop now has a WHILE instead of a FOR loop. This is easier to control. Also, the STEP 2 has been abandoned. The number-variable is incremented by 1 so also all even numbers are tested for primeness. Of course this is ridiculous but this makes detecting century breaks easier.

Compile the source prim03.mod with obc prim03.mod -o prim03

       1       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
     ...
    4703    4721    4723    4729    4733    4751    4759    4783    4787    4789    4793    4799
    4801    4813    4817    4831    4861    4871    4877    4889
    4903    4909    4919    4931    4933    4937    4943    4951    4957    4967    4969    4973    4987    4993    4999

Total primes =    670
Iterations   =    229.092
   
670 primes in the range 0 - 5000. And it took 230 thousands iterations through the isPrime loop. Execution was fast: less than a second. Of course I tested with bigger numbers as well:
  268123  268133  268153  268171  268189  268199
  268207  268211  268237  268253  268267  268271  268283  268291  268297
  268343
  268403  268439  268459  268487  268493
  413353
  461707  461717
  499801  499819  499853  499879  499883  499897
  499903  499927  499943  499957  499969  499973  499979

Total primes =  41539
Iterations   = 864337884
   
In the range 0 - 500.000 there are 41.539 primes and the isPrime loop is used 864.337.884 times (almost a billion times). It took the program close to 10 seconds to execute. Notice centuries 2683 and 4133. These centuries contain just one prime number. There must be centuries that hold no primes at all I guess.
Also notice the iterations count. 0-5000 is 230k iterations and 0-500000 is 864337k. So multiplying the range by 100 multiplies the iteration count by 4000.

Speeding things up

By now it is evident that most of the time is spent in the isPrime WHILE loop. For moderate numbers, say upto 10.000, it is no big deal, but for testing '499979' for primeness it takes 41k iterations. That's a lot. Do we really need all these tests? No we don't. It is enough to test num for divisors upto SQRT (num). So at first I wanted to have an integer square root routine. It has been made. And it works. Find it in the navigator frame on the right.

But there is a better and easier way. Let's have a look at the divisors of the number 105:

   1 x 105		105 x 1
   3 x 35		 35 x 3
   5 x 21		 21 x 5
   7 x 15		 15 x 7
   
The same number combinations appear twice: once below root (num) and once above. So we can stop testing when the divisor exceeds the root of num. And we can simply flag this condition when the quotient of num / prime gets less than the prime itself. So, 7 x 15 is followed by 15 x 7 and this is detected. Below is an example:
num = 499979		sqrt (499979) = 707

prime numbers around tipping point : 677   683   691   701   709   719

 num	 div   quo      quotient > divisor
------	 ---   ---
499979 / 677 = 738		yes
499979 / 683 = 732		yes
499979 / 691 = 723		yes
499979 / 701 = 713		yes
499979 / 709 = 705		no
499979 / 719 = 695		no
   
This way it is easy to detect the flipping point (around sq root of the number) without the need to determine that root. Notice that the number 499979 need only be tested until number 707 instead of until number 499973 (prime number 42000) so for the higher numbers this saves around 40000 tests per number.

Prim04 : the source

Below is the source code of the new Prim04 program. The changes are in the isPrime routine. In the old version, adding new primes was just adding it at the current array index. But since we now only go through a versy small part of the array, it is needed to add it at the end and to uodate the end position.

MODULE prim04;

IMPORT	Out;

CONST	maxIndex	= 1000000;
	maxCnt		= 16;

VAR	primes		: ARRAY maxIndex OF INTEGER;
	m, n, total, t	: INTEGER;


PROCEDURE isPrime (num : INTEGER) : BOOLEAN;

VAR	i, quo, pri	: INTEGER;

BEGIN
  i := 1;
  LOOP
    pri := primes [i];
    IF  pri = 0  THEN  EXIT  END;
    IF  num MOD pri = 0  THEN  RETURN FALSE  END;
    quo := num DIV pri;
    IF  quo < pri  THEN  EXIT  END;
    INC (t);
    INC (i)
  END;
  i := total;
  primes [i] := num;
  INC (total);
  RETURN  TRUE
END isPrime;
			 

PROCEDURE Init;

VAR	i	: INTEGER;

BEGIN
  m := 2;
  t := 0;
  primes [0] := 1;
  primes [1] := 2;
  total := 2;
  FOR  i := 2  TO  maxIndex - 1  DO  primes [i] := 0  END;
  Out.String ("     1     2")  
END Init;


BEGIN
  Init;
  n := 3;
  WHILE  n < 500000  DO
    IF  isPrime (n)  THEN
      Out.Int (n, 8);
      INC (m);
      IF  m = maxCnt  THEN
        Out.Ln;
	m := 0
      END
    END;
    INC (n);
    IF  n MOD 100 = 0  THEN  
      Out.Ln;
      m := 0
    END
  END;
  Out.Ln;
  Out.String ("Total primes = ");	Out.Int (total, 6);	Out.Ln;
  Out.String ("Iterations   = ");	Out.Int (t, 9);		Out.Ln
END prim04.
   
Now when we run this program for 0 - 500.000 it takes less than a second to run and produces:
prim03.mod
   
       1       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

    4703    4721    4723    4729    4733    4751    4759    4783    4787    4789    4793    4799
    4801    4813    4817    4831    4861    4871    4877    4889
    4903    4909    4919    4931    4933    4937    4943    4951    4957    4967    4969    4973    4987    4993    4999

Total primes =    670
Iterations   =    229.092

-------------

prim04.mod

       1       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
     
    4703    4721    4723    4729    4733    4751    4759    4783    4787    4789    4793    4799
    4801    4813    4817    4831    4861    4871    4877    4889
    4903    4909    4919    4931    4933    4937    4943    4951    4957    4967    4969    4973    4987    4993    4999

Total primes =    670
Iterations   =     14.600
   
230k iterations reduced to less than 15k. That is 93% less. This explains why prim04 needs close to half a second where prim03 still required around 10.

Push it to the limits

I changed prim04 so that it determines the prime numbers upto 10 million. Prim03 would probably take hours to run, but the new prim04 completes in 6 seconds producing:

       1       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
    ...
 9548909 9548921 9548927 9548947 9548953 9548993 9548999

 9549101 9549103 9549107 9549109 9549119 9549157 9549167
 9707911 9707933 9707947 9707977 9707987 9707989

 9708113 9708121 9708137 9708143 9708187 9708191
 9953507 9953513 9953521 9953533 9953539 9953557 9953569 9953579
 9953621
 9953711 9953717 9953719 9953737 9953743 9953747 9953753
 9967801 9967813 9967819 9967823 9967883 9967891
 9967987
 9968009 9968029 9968041 9968053
 9999601 9999637 9999653 9999659 9999667 9999677
 9999713 9999739 9999749 9999761
 9999823 9999863 9999877 9999883 9999889
 9999901 9999907 9999929 9999931 9999937 9999943 9999971 9999973 9999991

Total primes = 664.580
Iterations   = 276.809.518
   
In the first 10 million numbers, there are close to 665 thousand primes and the software required a quarter of a billion iterations to find this out. And all of it in less than 10 seconds using an interpreted language. I wonder how long this routine would take when it was written in Python.

Take a look at centuries 95490 and 97080. These centuries contain zero primes. Other centuries contain just one or two primes. Which is perfectly understandable since there are around 3000 primes (and hence possible divisors) at that stage. And this raises the question: are there empty millenia? Yes of course there are. But where?

With the new program 'factor' I checked one of these centuries (no need to chwck even numbers and those ending at '5') and indeed all numbers in that century were prime.

Is '1' a prime?

A prime is a number divisible by '1' and by itself. So in this sense, '1' is a prime since it is divisible by '1' and by '1'. But it's a silly prime. So '1' will be skipped from the primes list as of now (in this blog).
Next dilemma; is '2' a prime? It meets the criteria, both divisors are not the same, but there are no values between '1' and '2' so, it's a bit cheating. '2' is a silly prime as well. It is prime, but it is divisible by any number smaller than itself (i.e '1'). Let's call it a borderline condition. It gets the benefit of the doubt.

Getting an argument

The current version, prim04, needs the source to be changed and recompiled when a new range of primes is required. That's not neat in programming. I want ti be able to use the command line for specifying my range of numbers. So prim05 was created. It uses the same 'engine' as prim04 but now with another bodywork.

Prim05 : the source code

MODULE prim05;

IMPORT	Args, Conv, Out;

CONST	maxIndex	= 1000000;
	maxCnt		= 16;

VAR	primes		: ARRAY maxIndex OF INTEGER;
	m, n, total, t,
	maxNum		: INTEGER;


PROCEDURE isPrime (num : INTEGER) : BOOLEAN;

VAR	i, quo, pri	: INTEGER;

BEGIN
  i := 0;
  LOOP
    pri := primes [i];
    IF  pri = 0  THEN  EXIT  END;
    IF  num MOD pri = 0  THEN  RETURN FALSE  END;
    quo := num DIV pri;
    IF  quo < pri  THEN  EXIT  END;
    INC (t);
    INC (i)
  END;
  i := total;
  primes [i] := num;
  INC (total);
  RETURN  TRUE
END isPrime;


PROCEDURE GetParams;

VAR	n, count	: INTEGER;
	option		: ARRAY 33 OF CHAR;

BEGIN
  count := Args.argc;
  IF  count = 1  THEN  
    maxNum := 1000
  ELSE
    Args.GetArg (1, option);
    maxNum := Conv.IntVal (option)
  END;
END GetParams;


PROCEDURE Init;

VAR	i	: INTEGER;

BEGIN
  m := 1;
  t := 0;
  primes [0] := 2;
  total := 1;
  GetParams;
  FOR  i := 2  TO  maxIndex - 1  DO  primes [i] := 0  END;
  Out.String ("      2")
END Init;


BEGIN
  Init;
  n := 3;
  WHILE  n < maxNum  DO
    IF  isPrime (n)  THEN
      Out.Int (n, 8);
      INC (m);
      IF  m = maxCnt  THEN
        Out.Ln;
	m := 0
      END
    END;
    INC (n);
    IF  (n MOD 100 = 0)  THEN  
      Out.Ln;
      m := 0
    END
  END;
  Out.Ln;
  Out.String ("Total primes = ");	Out.Int (total, 6);	Out.Ln;
  Out.String ("Iterations   = ");	Out.Int (t, 9);		Out.Ln
END prim05.
   
As you can see there are two options: Since the processing engine is the same for prim04 and prim05 it doens't make sense to show sample outputs.

Bigger is better

I want it all and I want it now.... I want to be able to process bigger number. I want to find bigger prime numbers. What would be the easiest way to accomplish this?
I could make one easy tune-up: use LONGINT's (64 bit integers, ranging from roughly 0 to 10^19 and the negative numbers) That would open the gates to very large prime numbers. There's only one drawback. Which speed penalty will I have to pay for this since I am running on a 32 bit operating system, so LONGINT's will have to be processed in software (an INTEGER is held in one register, a LONGINT will occupy two registers).

It's probably best to just try it out.

The LONGINTed source

MODULE prim06;

IMPORT	Args, Conv, Out;

CONST	maxIndex	= 10000000;
	maxCnt		= 4;

VAR	primes		: ARRAY maxIndex OF LONGINT;
	m, n, t,
	maxNum		: LONGINT;
	total		: INTEGER;


PROCEDURE isPrime (num : LONGINT) : BOOLEAN;

VAR	i		: INTEGER; 
	quo, pri	: LONGINT;

BEGIN
  i := 0;
  LOOP
    pri := primes [i];
    IF  pri = 0  THEN  EXIT  END;
    IF  num MOD pri = 0  THEN  RETURN FALSE  END;
    quo := num DIV pri;
    IF  quo < pri  THEN  EXIT  END;
    INC (t);
    INC (i)
  END;
  i := total;
  primes [i] := num;
  INC (total);
  RETURN  TRUE
END isPrime;


PROCEDURE GetParams;

VAR	n, count	: INTEGER;
	option		: ARRAY 33 OF CHAR;

BEGIN
  count := Args.argc;
  IF  count = 1  THEN  
    maxNum := 1000
  ELSE
    Args.GetArg (1, option);
    maxNum := Conv.IntVal (option)
  END;
END GetParams;


PROCEDURE Init;

VAR	i	: INTEGER;

BEGIN
  m := 1;
  t := 0;
  primes [0] := 2;
  total := 1;
  GetParams;
  FOR  i := 2  TO  maxIndex - 1  DO  primes [i] := 0  END;
  Out.String ("      2")
END Init;


BEGIN
  Init;
  n := 3;
  WHILE  n < maxNum  DO
    IF  isPrime (n)  THEN
      Out.LongInt (n, 24);
      INC (m);
      IF  m = maxCnt  THEN
        Out.Ln;
	m := 0
      END
    END;
    INC (n);
    IF  (n MOD 100 = 0)  THEN  
      Out.Ln;
      m := 0
    END
  END;
  Out.Ln;
  Out.String ("Total primes = ");	Out.LongInt (total, 6);		Out.Ln;
  Out.String ("Iterations   = ");	Out.LongInt (t, 9);		Out.Ln
END prim06.
   

See it run

I did several cross comparisons. Both prim04 and prim06 produce the same results so that is good (and expected since both programs use the same engine). But how about the speed? For small prime ranges there is virtually no speed penalty. But for 1 - 10 miilion range the speed penalty is considerable:

oxygen:~/Oberon/prime$ time ./prime5 100000 >jj		100k, INTEGERs
	real    0m0.092s
	user    0m0.029s
	sys     0m0.005s

oxygen:~/Oberon/prime$ time ./prime6 100000 >jj		100k, LONGINTs
	real    0m0.147s
	user    0m0.093s
	sys     0m0.005s

oxygen:~/Oberon/prime$ time ./prime5 1000000 >jj	1 million, INTEGERs
	real    0m0.352s
	user    0m0.345s
	sys     0m0.007s

oxygen:~/Oberon/prime$ time ./prime6 1000000 >jj	1 million, LONGINTs
	real    0m1.568s
	user    0m1.561s
	sys     0m0.007s

oxygen:~/Oberon/prime$ time ./prime5 10000000 >jj	10 million, INTEGERs
	real    0m6.421s
	user    0m6.402s
	sys     0m0.021s

oxygen:~/Oberon/prime$ time ./prime6 10000000 >jj	10 million, LONGINTs
	real    0m31.173s
	user    0m31.162s
	sys     0m0.023s
   
I also issued the command prime5 100000000 >kk to look for primes upto 100 million. The program took 12 minutes to complete. The last section of the output was:
	99999509                99999517                99999539                99999541
	99999547                99999551                99999563                99999587		99999589
	99999611                99999617		99999623                99999643		99999677
	99999703                99999721                99999773                99999787
	99999821                99999827                99999839                99999847
	99999931                99999941		99999959                99999971		99999989
	
   Total primes =     5.761.455
   Iterations   = 6.226.689.926
   
The file 'kk' was 134 MB big. It would probably not fit in your scrollback buffer. So if you want to remember a big prime number, just use 99999989 6 nines an 8 and a 9.

Downloads

When downloading, seamonkey will try to interpret these files as MOD files, a sound format from the Amiga, AFAIK.

Page created 01 Feb 2017,