/*---------------------------------------------------------------------- File : prime.c Contents: test whether a given number is prime Author : Christian Borgelt History : 05.11.1997 file created 09.11.1997 remark on 'i = 6-i;' versus 'i ^= 6;' added 12.11.1997 exercise on conditional compilation added 14.11.1997 check for (n < 2) moved to 'prime' functions ----------------------------------------------------------------------*/ #include #include /*---------------------------------------------------------------------- Functions ----------------------------------------------------------------------*/ int prime1 (int n) { /* --- the simplest approach */ int d; /* divisor to test */ if (n < 2) return 0; /* 2 is the smallest prime number */ for (d = 2; d < n; d++) /* test all numbers less than n */ if (n % d == 0) return 0; /* if a divisor found, n is not prime */ return 1; /* if no divisor found, n is prime */ } /* prime1() */ /*---------------------------------------------------------------------- This approach directly implements the definition of a prime number: it is divisible only by one and itself. That is, it is not divisible by any number in {2,...,n-1}, where n is the prime number. We first check whether n is less than 2 to avoid a wrong result for n = 1. Then we test all numbers in {2,...,n-1} whether they divide n. If we find any such number, we know immediately that n cannot be prime, and hence abort the function with a negative result (i.e. 0 = false). Only if all numbers in {2,...,n-1} do not divide n without a remainder, we return a positive result (i.e. 1 = true). ----------------------------------------------------------------------*/ int prime2 (int n) { /* --- skip even numbers */ int d; /* divisor to test */ if (n < 2) return 0; /* 2 is the smallest prime number */ if (n <= 2) return 1; /* 2 is a prime number */ if (n % 2 == 0) return 0; /* if n is even, it is not prime */ for (d = 3; d < n; d += 2) /* test all odd numbers less than n */ if (n % d == 0) return 0; /* if a divisor found, n is not prime */ return 1; /* if no divisor found, n is prime */ } /* prime2() */ /*---------------------------------------------------------------------- First improvement on the algorithm used in prime1. It is obvious that any number which is divisible by an even number must also be divisible by two, i.e. must be even itself. Hence any test whether n is divisible by d, where d is even and not equal to two is pointless. We would have already found out that n is not prime when trying to divide by two, i.e. when testing whether n is even. Hence we perform only one test whether n is even before the loop and then in the loop only look at odd numbers d. This is easily achieved by starting with d = 3 and increasing d by 2 in each step. ----------------------------------------------------------------------*/ int prime3 (int n) { /* --- skip multiples of 2 and 3 */ int d; /* divisor to test */ int i; /* increment for divisor */ if (n < 2) return 0; /* 2 is the smallest prime number */ if (n <= 3) return 1; /* 2 and 3 are prime numbers */ if (n % 2 == 0) return 0; /* if n is even or if it is */ if (n % 3 == 0) return 0; /* divisible by 3, it is not prime */ for (d = 5, i = 2; d < n; d += i, i ^= 6) if (n % d == 0) return 0; /* if a divisor found, n is not prime */ return 1; /* if no divisor found, n is prime */ } /* prime3() */ /*---------------------------------------------------------------------- Second improvement on the algorithm used in prime. The idea is the same as in prime2, only taken one step further. In addition to multiples of two (even numbers) we also skip multiples of three in the loop and perform a single test, whether n is divisible by three before starting the loop. However, the loop itself gets much more complicated now, since we cannot use a fixed increment on the divisor any longer (in prime2, we could increase the divisor by two in each step to traverse the odd numbers). It turns out that to skip multiples of two and three, we have to use an alternating scheme, increasing d by two in one step, four in the next, then two again, and so forth. If you do not understand this directly, recall that any multiple of six is a multiple of two as well as of three. Hence, the scheme must repeat itself in cycles of length 6 --- which it does: 2 + 4 = 6. Recall also that between any multiple of 6 and the next, there are two numbers that are divisible by two and one that is divisible by three. Think about how these numbers are placed between the first multiple of six and the next and the scheme will become obvious. Now for the implementation of this alternating scheme. A direct approach would use an increment i set to 2 in the beginning (if starting with d = 5), test in each step whether i is 2, set it to 4, if this is the case, and to 2 otherwise. This approach works fine, but we can do better. One possibility is to compute in each step 'i = 6 -i'. As you can easily verify, this produces the alternating scheme. Another approach, which is used in the above function, exploits the way in which numbers are represented in the computer. They are stored as binary numbers (cf. exercise 11) with 2 having bit 1 set (since 2 = 2^1) and 4 having bit 2 set (since 4 = 2^2). This special representation can be used to produce the alternating scheme of increments of two and four by applying an operation that changes both bits to their other possible state. An operation to switch the value of a bit is the logical xor function (exclusive or, represented in C by '^') with a bit that is set to one. Here is the truth table of the exclusive or operation: a^b 0 1 Let a be an arbitrary bit (0 or 1) and let b = 1. 0 0 1 Then computing a = a^b just flips the bit in a. 1 1 0 Doing an xor with b = 0 leaves the bit in a unchanged. When applying the xor operation to numbers, it is applied to each bit separately. So we only need to find a value for b that has a bit set for each bit we want to flip in a and that is zero for any other bit. For our increment we need to flip bits 1 and 2, so we need a number with bits 1 and 2 set and all other bits zero. This number is 6 = 2^2 + 2^1 = 4 + 2. And that's why we have this strange 'i ^= 6' in the for-statement. (Remark: You may think that it is always better to use the statement 'i = 6-i;' instead of 'i ^= 6;', because it is much easier to understand. In general you are right, since your first concern should be to write programs that are easy to understand and to maintain. However, 'i ^= 6;' may be preferable, if you are coding for speed. Although on a Sun Sparc processor, both commands take the same time to execute, on an Intel Pentium 'i ^= 6;' takes only two thirds of the time it takes to execute 'i = 6-i;'. Of course, another reason why 'i ^= 6;' is used here is to teach you the bitwise xor operation.) It is clear that we could take the idea we exploited in prime2 and prime3 even further and skip multiples of 5, 7, 11 etc. in the loop. In the limit we thus arrive at the sieve of Eratosthenes (see program eratos.c, to be found here later). ----------------------------------------------------------------------*/ int prime4 (int n) { /* --- test only up to n/2 */ int d; /* divisor to test */ if (n < 2) return 0; /* 2 is the smallest prime number */ for (d = 2; d+d <= n; d++) /* test numbers up to n/2 */ if (n % d == 0) return 0; /* if a divisor found, n is not prime */ return 1; /* if no divisor found, n is prime */ } /* prime4() */ /*---------------------------------------------------------------------- In this function a different idea is exploited to speed up the prime number test. It is based on the question whether it is necessary to test numbers d up to n-1 or whether it is possible to stop earlier. To find an answer to this question, consider the number n-1 for an arbitrary n. Is n-1 a divisor? Obviously, it is a divisor only for n = 2 (and then it is 1, so it is not tested). For any number n > 2, it is 1 < n/(n-1) < 2 and thus n-1 cannot be a divisor. Now look at n-2. As you can easily verify, it cannot be a divisor of any number n > 4, since for n > 4, it is 1 < n/(n-2) < 2. In the same way, we derive that n-3 cannot be a divisor of any number n > 6. Generalizing, we arrive at: n-k is not a divisor of n, if n > 2k. Hence we need to test numbers d only up to n/2. We test this using d+d <= n, since additions are often faster than multiplications or divisions. Although most compilers recognize automatically that for certain constant values (e.g. 2, 3 or 4) a division or a multiplication can be replaced by a more efficient command, you cannot be sure of this. So it is a good idea do replace at least 2*d by d+d as above. Actually, it may be even better to write 'd <= n-d' to avoid problems with the range of values for an integer number. We deliberately ignore such subtleties (just as we did, without comment, in prime2 and prime3). ----------------------------------------------------------------------*/ int prime5 (int n) { /* --- test only up to square root */ int d; /* divisor to test */ if (n < 2) return 0; /* 2 is the smallest prime number */ for (d = 2; d*d <= n; d++) /* test nums. up to square root of n */ if (n % d == 0) return 0; /* if a divisor found, n is not prime */ return 1; /* if no divisor found, n is prime */ } /* prime5() */ /*---------------------------------------------------------------------- This function carries the idea exploited in prime4 one step further. To understand the method, assume that you have an arbitrary natural number n, its square root r, and a divisor d of n with d >= r. Since d is a divisor, the quotient n/d must also be a divisor. Now d >= r and r*r = n, hence n/d <= r. That is, before we could reach the divisor d and test it, we would have tested n/d, found that it divides n, and aborted the function. It follows that we need not test numbers larger than r. If there is a divisor d >= r, there must also be divisor d' <= r. Another way to understand this method is to start from function prime4. We derived in the comment on prime4 that no number d > n/2 can be a divisor. But to humbly assert such a weak condition is necessary only as long as we do not know whether n is even. If we already know that n is not divisible by 2, we can strengthen the assertion to: There cannot be a divisor d > n/3. This is so, because for any such d we can easily derive that 1 < n/d < 3. Hence, to be a divisor, n/d must be 2. But then 2 must also be a divisor, which it is not. And then, when we tested whether n is divisible by 3 and found that it is not, we can again reduce the limit for d, namely to n/5 (since 4 cannot be a divisor because 2 is no divisor). In the limit we arrive again at the condition d*d <= n used above. To implement this idea we simply change the condition in the for-statement to d*d <= n (compared to prime1 and prime4). In analogy to prime4 it may be better to write d <= n/d to avoid problems with the range of values for an integer number. Again we deliberately ignore such subtleties. ----------------------------------------------------------------------*/ int prime6 (int n) { /* --- prime3 and prime5 combined */ int d; /* divisor to test */ int i; /* increment of divisor */ if (n < 2) return 0; /* 2 is the smallest prime number */ if (n <= 3) return 1; /* 2 and 3 are prime numbers */ if (n % 2 == 0) return 0; /* if n is even or if it is */ if (n % 3 == 0) return 0; /* divisible by 3, it is not prime */ for (d = 5, i = 2; d*d <= n; d += i, i ^= 6) if (n % d == 0) return 0; /* if a divisor found, n is not prime */ return 1; /* if no divisor found, n is prime */ } /* prime6() */ /*---------------------------------------------------------------------- Of course the ideas of prime3 and prime5 can be combined. This is done in function prime6. It tests the numbers up to the square root of n skipping multiples of 2 and 3. ----------------------------------------------------------------------*/ int main (int argc, char *argv[]) { /* --- main function */ int n; /* number to test */ if (argc != 2) { /* if wrong number of arguments given */ printf("usage: %s n\n", argv[0]); printf("test whether n is a prime number\n"); return 0; /* print a usage message */ } /* and abort the program */ n = atoi(argv[1]); /* convert argument to a number */ if (n < 1) { /* test for a correct argument */ printf("%s: illegal number\n", argv[0]); return 1; } printf("prime1(%i) = %i\n", n, prime1(n)); printf("prime2(%i) = %i\n", n, prime2(n)); printf("prime3(%i) = %i\n", n, prime3(n)); printf("prime4(%i) = %i\n", n, prime4(n)); printf("prime5(%i) = %i\n", n, prime5(n)); printf("prime6(%i) = %i\n", n, prime6(n)); return 0; /* return 'ok' */ } /* main() */ /*---------------------------------------------------------------------- Note that the number of arguments must be 2, since the first argument is always the program's name (cf. program heron.c). This explains also, why we ignore argv[0] (which contains the name) and convert only argv[1] to a number. Note also that this program does not detect misspelled numbers. If it is called with 'prime 7xc8', then n will simply be 7. To see the difference in running time of the different functions try 'prime 3537151'. Exercise: Use conditional compilation to compile programs each of which runs only one of the six versions. (Conditional compilation is explained in the file `dow.c'. Hint: Use a preprocessor constant VERSION set to the version number of the prime number test and test that number with '#if (VERSION == x) ... #endif', where x is a version number.) ----------------------------------------------------------------------*/