Umeå universitet
#include "stdio.h"
#include "sys/time.h"
#include "stdlib.h"
long dclock() {
/* Returns the time in 10^-3 seconds */
struct timeval tp;
struct timezone tzp;
gettimeofday(&tp, &tzp);
return 1000 * (tp.tv_sec % 1000) + tp.tv_usec / 1000; }
void sieve(int n) {
unsigned int i, j, primes = 0, ops = 0;
long time;
char* mark = (char*) malloc(n * sizeof(char));
time = - dclock();
memset(mark, 1, n; /* Set all marks to 1 */
for (i = 2; i < n; i++)
if (mark[i]) { /* i is prime */
primes++;
for (j = 2 * i; j < n; j += i) {
ops++;
mark[j] = 0; } }
time += dclock();
printf("n = %10d, n / primes = %8.3f, ops / n = %8.3f, time = %8.3f\n",
n, (float) n / primes, (float) ops / n, 0.001 * time);
free(mark); }
int main() {
int n;
for (n = 4; 1; n *= 2)
sieve(n);
return 1; }
Running the program on a Pentium III with a 1 GHz processor gives the following output:
n = 4, n / primes = 2.000, ops / n = 0.000, time = 0.000 n = 8, n / primes = 2.000, ops / n = 0.375, time = 0.000 n = 16, n / primes = 2.667, ops / n = 0.812, time = 0.000 n = 32, n / primes = 2.909, ops / n = 1.031, time = 0.000 n = 64, n / primes = 3.556, ops / n = 1.297, time = 0.000 n = 128, n / primes = 4.129, ops / n = 1.500, time = 0.000 n = 256, n / primes = 4.741, ops / n = 1.680, time = 0.000 n = 512, n / primes = 5.278, ops / n = 1.826, time = 0.000 n = 1024, n / primes = 5.953, ops / n = 1.961, time = 0.000 n = 2048, n / primes = 6.628, ops / n = 2.079, time = 0.000 n = 4096, n / primes = 7.262, ops / n = 2.185, time = 0.000 n = 8192, n / primes = 7.969, ops / n = 2.281, time = 0.000 n = 16384, n / primes = 8.623, ops / n = 2.369, time = 0.001 n = 32768, n / primes = 9.330, ops / n = 2.450, time = 0.001 n = 65536, n / primes = 10.018, ops / n = 2.525, time = 0.001 n = 131072, n / primes = 10.699, ops / n = 2.595, time = 0.003 n = 262144, n / primes = 11.398, ops / n = 2.660, time = 0.012 n = 524288, n / primes = 12.083, ops / n = 2.721, time = 0.083 n = 1048576, n / primes = 12.784, ops / n = 2.779, time = 0.240 n = 2097152, n / primes = 13.477, ops / n = 2.834, time = 0.523 n = 4194304, n / primes = 14.172, ops / n = 2.885, time = 1.135 n = 8388608, n / primes = 14.869, ops / n = 2.934, time = 2.401 n = 16777216, n / primes = 15.565, ops / n = 2.981, time = 4.984 n = 33554432, n / primes = 16.259, ops / n = 3.026, time = 10.360 n = 67108864, n / primes = 16.956, ops / n = 3.069, time = 21.349 n = 134217728, n / primes = 17.652, ops / n = 3.110, time = 44.735 n = 268435456, n / primes = 18.347, ops / n = 3.149, time = 93.482One clearly can observe how n / primes increases in a logarithmic way with n. The number of primes apparantly is something like 0.70 * n / log_2 n ~= n / ln n, which might be the actual value. The time increases strongly, but this is also due to the decreasing effective speed of the processor when working on larger sets of data. This is not what interests us most. ops / n is more interesting, this is what determines the time-order of the sieve method. It is not constant, but grows slowly, much slower than logarithmic. Which function will it be?
Theorem: Erathostenes algorithm runs in O(n * loglog n) time.
Proof: If a number i is a prime, then we are processing about n / i numbers, costing O(n / i) time. If i is a non-pirme, it can be handled in O(1) time. Thus, the costs caused by all non-primes can be estimated on O(n). Now we will try to bound the costs due to the primes. Denote the number of primes in the interval [m, 2 * m] by P(m). Mathematics gives us P(m) = Theta(m / log m). Let us assume this as a fact. Particularly, there are constants m' and c so that P(m) <= c * m / log m, for all m >= m'. The contribution to the total cost of the primes smaller than m' is bounded by sum_{i = 2}^{m'} n / i = O(log m' * n) = O(n), because m' is a constant. So, we may assume that P(m) <= c * m / log m for all m. The contribution to the total cost by the primes in an interval [m, 2 * m] are given by sum_{m <= i <= 2 * m| i is prime} n / i <= 2 * n / m * P(m) <= 2 * c * n / log m. Now, we have succeeded in eliminating the uncertainty about the number and distribution of the primes, and the contribution to the cost by them can be estimated with a normal sum: sum_{k = 0}^{log n} sum_{m = 2^k <= i <= 2 * m| i is prime} n / i <= 2 * c * n sum_{k = 0}^{log n} 1 / k ~= 2 * c * n * loglog n.
loglog n is a very slowly growing function. loglog 2^4 = 2, loglog 2^8 = 3, loglog 2^16 = 4, loglog 2^32 = 5. The above sequence of values is quite well described by 0.65 * log_2 log_2 (n). For all practical puproses one might work with ops <= 4 * n. The efficiency of this method (and the speed of modern computers) is exemplified by the fact that all primes up to 2^28 could be computed in less than 100 seconds. The memory consumption can be reduced by working with individual bits but even then, we are running out of memory soon. Using the secondary memory, we might go up to 2^35, but then it is really out and over. Clearly for larger numbers one must return either to division-based primality testing or to the randomized tests presented hereafter.
For the sake of completeness, we here give the program that works with bits instead of chars and the corresponding table.
#include "stdio.h"
#include "sys/time.h"
#include "stdlib.h"
#include "string.h"
char mask[8];
char invt[8];
long dclock() {
/* Returns the time in 10^-3 seconds */
struct timeval tp;
struct timezone tzp;
gettimeofday(&tp, &tzp);
return 1000 * (tp.tv_sec % 1000) + tp.tv_usec / 1000; }
void sieve(int n) {
unsigned int i, j, primes = 0, ops = 0;
long time;
char* mark = (char*) malloc((n + 7) / 8);
time = - dclock();
memset(mark, 255, (n + 7) / 8); /* Set all marks to 1 */
for (i = 2; i < n; i++)
if (mark[i >> 3] & mask[i & 7]) { /* i is prime */
primes++;
for (j = 2 * i; j < n; j += i) {
ops++;
mark[j >> 3] = mark[j >> 3] & invt[j & 7]; } }
time += dclock();
printf("n = %10d, n / primes = %8.3f, ops / n = %8.3f, time = %8.3f\n",
n, (float) n / primes, (float) ops / n, 0.001 * time);
free(mark); }
int main() {
int i, n;
for (i = 0, n = 1; i < 8; i++, n *= 2) {
mask[i] = n;
invt[i] = 255 - n; }
for (n = 4; 1; n *= 2)
sieve(n);
return 1; }
n = 4, n / primes = 2.000, ops / n = 0.000, time = 0.000 n = 8, n / primes = 2.000, ops / n = 0.375, time = 0.000 n = 16, n / primes = 2.667, ops / n = 0.812, time = 0.000 n = 32, n / primes = 2.909, ops / n = 1.031, time = 0.000 n = 64, n / primes = 3.556, ops / n = 1.297, time = 0.000 n = 128, n / primes = 4.129, ops / n = 1.500, time = 0.000 n = 256, n / primes = 4.741, ops / n = 1.680, time = 0.000 n = 512, n / primes = 5.278, ops / n = 1.826, time = 0.000 n = 1024, n / primes = 5.953, ops / n = 1.961, time = 0.000 n = 2048, n / primes = 6.628, ops / n = 2.079, time = 0.000 n = 4096, n / primes = 7.262, ops / n = 2.185, time = 0.000 n = 8192, n / primes = 7.969, ops / n = 2.281, time = 0.000 n = 16384, n / primes = 8.623, ops / n = 2.369, time = 0.001 n = 32768, n / primes = 9.330, ops / n = 2.450, time = 0.001 n = 65536, n / primes = 10.018, ops / n = 2.525, time = 0.002 n = 131072, n / primes = 10.699, ops / n = 2.595, time = 0.004 n = 262144, n / primes = 11.398, ops / n = 2.660, time = 0.007 n = 524288, n / primes = 12.083, ops / n = 2.721, time = 0.017 n = 1048576, n / primes = 12.784, ops / n = 2.779, time = 0.044 n = 2097152, n / primes = 13.477, ops / n = 2.834, time = 0.164 n = 4194304, n / primes = 14.172, ops / n = 2.885, time = 0.561 n = 8388608, n / primes = 14.869, ops / n = 2.934, time = 1.450 n = 16777216, n / primes = 15.565, ops / n = 2.981, time = 3.250 n = 33554432, n / primes = 16.259, ops / n = 3.026, time = 6.993 n = 67108864, n / primes = 16.956, ops / n = 3.069, time = 14.798 n = 134217728, n / primes = 17.652, ops / n = 3.110, time = 31.066 n = 268435456, n / primes = 18.347, ops / n = 3.149, time = 64.809 n = 536870912, n / primes = 19.043, ops / n = 3.187, time = 135.690 n = 1073741824, n / primes = 19.738, ops / n = 3.224, time = 282.884All important numbers are the same, as they should be. The changes in the program are actually very small. At first i am always afraid of addressing the individual bits, but once you start it is almost trivial and it really saves a lot. It is to be noticed that in this way one does not only save memory, but even time. The reason is that the time consumption is mainly due to accessing the memory, saving memory access time therefore outweighs the additional cost due to the more complex address computation. In less then 5 minutes we now find all primes up to 2^30, with a program that takes less than one hour to write down. Quite impressive!
Proving the theorem was hard because we do not precisely know how the primes are distributed. In many problems of this kind the uneven distribution is coped with by considering a number of experiments (as in the analysis of quicksort) for a range of values (as here), and then arguing like: "among these experiments / in this range there are at most so many cases". Each of the cases is subsequently assumed to contribute as much as the maximum possible contribution in the range.
For any prime number n and integer a, 1 <= a < n, a^{n - 1} mod n = 1.
We are not going to prove it, but will check a few cases. Let n = 7, then
1^6 = 1 == 0 * 7 + 1 2^6 = 64 == 9 * 7 + 1 3^6 = 729 == 104 * 7 + 1 4^6 = 4096 == 585 * 7 + 1 5^6 = 15625 == 2232 * 7 + 1 6^6 = 46656 == 6665 * 7 + 1For n = 9, we find
1^8 = 1 == 0 * 9 + 1 2^8 = 256 == 28 * 9 + 4 !!! 3^8 = 6561 == 729 * 9 + 0 !!! ...So, this looks like quite a good test for determining that numbers are NOT prime. A number a so that a^{n - 1} mod n != 1, is called a witness for the non-primality of n. For a non-prime number n, a number a with a^{n - 1} mod n == 1, is called a false witness for the primality of n: it suggests that n is a prime, though n is not. Many numbers have such false witnesses. For example, 4^14 mod 15 == 1: 4 is a false witness for the primality of 15. Most numbers have a few false witnesses, but there are numbers for which almost all a are false fitnesses. So, on one hand it is too expensive to test all a, and on the other hand, it does not give sufficient certainty to test only a few of them: it might happen that we falsely conclude that n is prime.
We have here a randomized algorithm with one-sided error: if we conclude that the number is not-prime, this conclusion is always correct, but the opposite conclusion might be wrong. This is no problem, and in the domain of randomized algorithms it is common to accept such a situation provided the probability on error can be made arbitrarily small. It is then up to the user to determine how much certainty he wants to have.
Theorem: For an arbitrary odd number n > 4, if n is a prime, then the improved Fremat test returns true for all 1 < a < n, else it return true for less than one quarter of the a.
Of course when testing a number for primality, and getting true, we are not satisfied to know that it is a prime with probability more than 1 - 1/4 = 3/4. But, the solution is simple: repeatedly perform the test for randomly chosen a. Then the probability that we get k times true for a number that is not a prime is less than (1/4)^k. So, the error probability (that is, concluding prime when it is not prime), can be made arbitrarily small.
Testing for primality in this way is a nice example of a so-called Monte Carlo algorithm: an algorithm that may give a wrong answer. We have learned here:
Public key cryptosystems are based on the fact that it is quite easy to find large prime numbers, while it is almost impossible to factorize the product of two large prime numbers. In these protocols, the receiver, in the book he is called Bob, takes the initiative. The sender, in the book called Alice, uses a key computed and sent by Bob to encode her message. Actually, because Bob makes this key public anyone else could also send encoded messages to Bob (quite handy: if we want to organize a longdistance group-meeting with k participants, then we need only k keys instead of Theta(k^2)).
Bob does the following:
For encoding a message, Alice does the following:
For decoding a message from Alice, Bob does the following:
All this is quite simple, the only issue is actually handling large numbers and finding large prime numbers. The time for encoding an n-digit message is determined by the expomod operation. Expomod performs a logarithmic number of multiplications and divisions of O(n)-bit sized numbers. Depending on how one implements these operations, each takes O(n^2) or O(n^1.58). So, with n-bit encoding the cost per transfered bit is O(n * log n) or O(n^1.58 * log n).
The key issue in a cryptosystem is of course not the time it takes but its security. The security in this case comes, at least we hope, from the impossibility for the ennemy to perform the decoding even though he knows c, n and z. Without knowing s, it appears very hard to compute c^s mod z. Let us believe this. s can be easily computed when one knows phi, but it appears hard otherwise. Let us believe this. How hard is it to compute phi when one knows z? phi = (p - 1) * (q - 1) = z + 1 - (p + q). So, one somehow must figure out the sum of the prime factors of z. If we believe that the sum of the prime factors of a number cannot be computed without knowing these factors, then knacking the code is as hard as factoring a Theta(n)-digit number. This is believed to be very hard.
The good thing about this idea is that it so far resisted attack even though many good scientists were thinking about it for a long time (but we cannot be sure, maybe we only believe it has not be knacked, while someone is getting very rich). The weakness is that in the line of reasoning there are many "believes". Therefore, it is still imaginable that some time in the future one of the made assumptions proves not to have been justified.