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:
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 = 864337884In 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.
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 7The 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 = 705This 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.no 499979 / 719 = 695no
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:
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.023sI 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.926The 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.
The python version
This week, my friend Patrick Brunier volunteered to make a version of prim06 in python. He did a fine job since the outputs of the python version and the obc version were identical. The python version just took longer to execute. Where the obc version complete in 12 minutes, the python version took
real 39m51.943s user 39m50.588s sys 0m0.813s40 minutes. Now that's a big penalty to pay for having to use a language that's in fashion just because it was named after Monty Python. Still, a big hand for Patrick. And here is his source.
# PRYM06 - a program to calculate prime numbers.
# Copyright (C) 2017 Patrick Brunier
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, ONLY version 3 of the License.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
## PRYM06 is a Python port of the Prim06.mod that was originally made by Jan Verhoeven using Oberon
## Have a look at his fantastic "Fruttenboel" site for the original Prim06.mod and many other interesting topics.
## https://fruttenboel.verhoeven272.nl/obc/obc10.html
## Goal of PRYM06 was to compare interpreted Python executing time to compiled Oberon execution time.
##
## The program uses some bad Python practices, like the use of globals. But the intention was, to make a 1:1 Python
## version of the original.
##
## Finding all primes up to 100.000.000 (maximum_number variable) took Prim06 (Oberon 3.0.7) 7 minutes compared to
## 32 minutes for Prym06 (Python 3.6.2).
from __future__ import print_function
# USER VARIABLES
maximum_number = 100000000
# CONSTANTS
maximum_index = 10000000
maximum_values_per_line = 4
# RUNTIME VARIABLES
found_primes = [0] * maximum_index
found_primes[0] = 2
number_of_found_primes = 1
number_of_iterations = 0
number_printed_on_line = 0
def is_prime(number):
global number_of_found_primes
global found_primes
global number_of_iterations
i = 0
while True:
prime = found_primes[i]
if prime == 0:
exit()
if (number % prime) == 0:
return False
quo = number / prime
if quo < prime:
break
i += 1
number_of_iterations += 1
i = number_of_found_primes
found_primes[i] = number
number_of_found_primes += 1
return True
print("%24d" % found_primes[0], end='')
number_printed_on_line += 1
n = 3
while n < maximum_number:
if is_prime(n):
print('%24d' % n, end='')
number_printed_on_line += 1
if number_printed_on_line == maximum_values_per_line:
print('')
number_printed_on_line = 0
n += 1
if n % 100 == 0:
print('')
number_printed_on_line = 0
print('')
print("Total primes = ", number_of_found_primes)
print("Iterations = ", number_of_iterations)
I get a headache going through this source. It's no wonder why Guido called his 'language' after a satirical comedy
show.
Speed tests
Since prim06 relies on brute force arithmetics, I used it as a single core benchmark. With some surprising results. First I compiled a statically linked version of primo06.mod
obc -C prim06.mod -o prim06staticand used it as a test.
time prim06static 100000000 >kkYes that is 100000000 or 100.000.000 or 100 million
Downloads
Page created 01 Feb 2017,