/*---------------------------------------------------------------------- File : perfect.c Contents: find all perfect numbers within a given interval Author : Christian Borgelt History : 05.11.1997 file created 17.11.1997 search interval made a program argument prime factor version added 18.11.1997 comments on prime factor version added ----------------------------------------------------------------------*/ #include #include #ifndef VERSION #define VERSION 1 /* default: compile first version */ #endif /*---------------------------------------------------------------------- In this program we use the conditional compilation facility provided by the preprocessor in a slightly different way compared to the program `dow.c' (in which conditional compilation was introduced). In that program one could change the version to be compiled only by changing the source code, i.e. by changing the definition of the preprocessor constant USESWITCH in the program itself. In this program we provide a way to change the version to be compiled without changing the source code. This is done by a conditional definition of the constant VERSION. We only define it, if it is not already defined. Whether some constant is already defined, can be checked with the preprocessor directives `#ifdef' and `#ifndef', which are short forms of `#if defined(x)' and `#if !defined(x)', where x is an arbitrary preprocessor constant. We define VERSION conditionally in order to have a fallback position, if no definition of VERSION is provided. How can we change the definition of VERSION without changing the source code? This is done with a compiler option. For the Gnu C-compiler it is the -D option (D for `define'). If you invoke the compiler on this program with 'gcc -D VERSION=2 ...', the constant VERSION will already be defined (and will have value 2), when the preprocessor processes the above directives. Hence the #define will not be executed and the second version will be compiled. ----------------------------------------------------------------------*/ /*---------------------------------------------------------------------- Functions ----------------------------------------------------------------------*/ #if (VERSION == 1) /* if to compile the first version */ int perfect (int n) { /* --- test whether n is perfect */ int d; /* divisor to try */ int sum; /* sum of divisors */ if (n <= 1) return 0; /* 1 is not a perfect number */ sum = 1; /* any number is divisible by one */ for (d = 2; d*d < n; d++) { /* try nums. up to square root of n */ if (n % d == 0) { /* if n is divisible by d, add d and */ sum += d +n/d; /* n/d to the sum (if d is a divisor, */ if (sum > n) return 0; /* n/d must also be a divisor) */ } /* if the divisor sum already exceeds */ } /* the number n, abort the function */ if (d*d == n) /* if n is a square number, */ sum += d; /* add the square root to the sum */ if (sum == n) return 1; /* only if the sum of the divisors */ return 0; /* equals the number, it is perfect */ } /* perfect() */ #endif /*---------------------------------------------------------------------- This function approaches the perfect number problem directly, i.e. it tests for divisors of the argument n and sums them. The computation is made efficient by testing for divisors only up to the square root of n. This is possible, since for any divisor d greater than the square root, there must be a corresponding divisor d' = n/d smaller than the square root. Hence if we have one, we can compute the other and need not test further numbers until we encounter it (cf. the program `prime.c' where the same trick is used to speed up a prime number test). With this approach, we only have to make sure that we do not add twice the square root of a square number. This is done by testing only numbers d whose square is less than n and dealing with square numbers separately. In analogy to the function prime4 (see file prime.c) it may be better to write d <= n/d to avoid problems with the range of values for an integer number. But just as in `prime.c' we deliberately ignore such subtleties. Exercise: There are some improvements other than testing divisors only up to the square root. These improvements can be derived from those used to speed up the prime number test, like testing only odd number as divisors, if the number itself is odd, etc. (cf. program prime.c). Implement these improvements. ----------------------------------------------------------------------*/ #if (VERSION == 2) /* if to compile the second version */ int perfect (int n) { /* --- test whether n is perfect */ int pfcnt; /* number of prime factors */ int pfs [9]; /* vector of prime factors */ int exps[9]; /* exponents of prime factors */ int des [9]; /* exponents in divisor */ int m = n; /* copy of the number to test */ int d; /* divisor to try */ int e; /* exponent of a prime factor */ int i; /* increment of divisor, index */ int sum; /* sum of divisors */ if (n <= 1) return 0; /* 1 is not a perfect number */ /* --- determine prime factors and their exponents --- */ pfcnt = 0; /* initialize prime factor counter */ if (m % 2 == 0) { /* if m is even, */ pfs[pfcnt] = 2; e = 0; /* then 2 is a prime factor */ do { e++; m /= 2; } while (m % 2 == 0); exps[pfcnt++] = e; /* divide out prime factor 2 */ } /* and note its exponent */ if (m % 3 == 0) { /* if m is divisible by 3, */ pfs[pfcnt] = 3; e = 0; /* then 3 is a prime factor */ do { e++; m /= 3; } while (m % 3 == 0); exps[pfcnt++] = e; /* divide out prime factor 3 */ } /* and note its exponent */ for (d = 5, i = 2; d*d <= m; d += i, i ^= 6) { if (m % d == 0) { /* if m is divisible by d, */ pfs[pfcnt] = d; e = 0; /* then d is a prime factor */ do { e++; m /= d; } while (m % d == 0); exps[pfcnt++] = e; /* divide out prime factor d */ } /* and note its exponent */ } if (m > 1) { /* store the last prime factor */ pfs[pfcnt] = m; exps[pfcnt++] = 1; } if (pfcnt <= 1) return 0; /* prime numbers are not perfect */ /* --- construct all divisors and sum them --- */ for (i = pfcnt; --i >= 0; ) /* copy all prime factor exponents */ des[i] = exps[i]; /* to the buffers for the divisor */ des[0]--; d = n /pfs[0]; /* compute first divisor and exps. */ sum = 1; /* any number is divisible by one */ do { /* -- divisor construction loop */ sum += d; /* sum divisor; if the sum exceeds */ if (sum > n) return 0; /* the number, abort the function */ for (i = 0; i < pfcnt; i++) { /* traverse divisor exponents */ if (des[i] > 0) { /* if an exponent can be reduced, */ des[i]--; d /= pfs[i]; /* reduce this exponent by one, */ break; /* compute the new divisor and */ } /* abort the loop */ e = des[i] = exps[i]; /* if an exponent is zero, */ while (--e >= 0) /* reset it to the original value */ d *= pfs[i]; /* and compute the new divisor by */ } /* multiplying the old by pf^exp */ } while (d > 1); /* repeat if there is another divisor */ if (sum == n) return 1; /* only if the sum of the divisors */ return 0; /* equals the number, it is perfect */ } /* perfect() */ #endif /*---------------------------------------------------------------------- In this function the problem is approached in a completely different way. The idea is to determine all prime factors of the number n that is to be tested and to construct all divisors of n from these prime factors. Hence the function has two clearly separated parts. In the first part the prime factors and their exponents are determined and stored in the vectors `pfs' and `exps', respectively. In the second part divisors are constructed by modifying the exponents of the prime factors. All possible combinations of exponents are generated (except one: the exponents as determined for the number n, which corresponds to divisor n). For each generated combination of exponents, the corresponding divisor is added to the sum (except one: all exponents zero, which corresponds to divisor 1 and thus terminates the loop --- divisor 1 is taken care of before starting the loop). The first part does not need much comment, since the trick used in the loop to test for divisors is already known from the program `prime.c'. We test for divisors 2 and 3 before entering the loop and in the loop we skip multiples of two and three. Every time we find a divisor, we determine its exponent by dividing out this divisor and counting how often it is present in the number. To determine the prime factors of n we use a copy of n (stored in m) since we need n to determine later on whether the sum of the divisors equals n. Hence we should not modify n. There are only three more questions that may arise about the prime factor computation. The first is: Why do we have nine fields in the vectors to store the prime factors and their exponents in? Well, the smallest number with more than nine different prime factors is 6469693230 = 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 29, and this number is already much larger than the largest number representable as an `int'. Hence any number representable as an `int' has less than 10 different prime factors (we need not worry about identical prime factors, since these are taken care of by the exponents). The second question is: In the loop, we only skip multiples of two and three. Hence we sometimes try a non-prime number as a divisor. Why do we end up with only prime numbers in the vectors? This is due to the fact that we divide out all prime factors found. For example, let n = 21025 = 5^2 * 29^2. This number is divisible by 25. 25 is not divisible by 2 or by 3 and is thus tested as a divisor in the loop. Hence, at first glance, one may think that 25 should appear in the vector `pfs'. But this is not the case, since before we reach 25 and test it, we already found that m (the copy of n) is divisible by 5 and divided out this prime factor. That is, when we try d = 25, it is m = 841 = 29^2, which is not divisible by 25 and thus 25 is skipped. It is easily verified that the factors of any non-prime number are divided out before this number is tried as a divisor, and thus no non-prime number is written into the vector `pfs'. The third question is: Why do we need this extra test on m > 1 after the loop? This is due to the fact that we stop the loop, if d*d > m. In this case the last prime factor, if it has exponent 1, is not divided out from m. Take n = 14 as an example. The loop is terminated with d = 5 and m = 7. If we did not check for m > 1, we would fail to notice the prime factor 7. Now for the second part of the above function, which is much more difficult to understand than the first. In this part, as indicated above, all combinations of exponents of the prime factors are generated and the corresponding divisors are summed. We start by copying the exponents of the prime factors of the number n to a buffer `des', which is intended to hold the prime factor exponents of the generated divisor. We then compute the first divisor as n / (smallest prime factor) and adapt the exponents in the vector `des' accordingly (we reduce the exponent of the smallest prime factor). In the loop, we maintain this relation. It is always d = \prod_{i=0}^{pfcnt-1} pfs[i]^des[i]. This is achieved by adapting d every time an exponent in the vector `des' is changed. The method to generate all possible combinations of exponents is based on the following scheme: We start with the vector of exponents for the first divisor. Then we reduce the exponent of the smallest prime factor, which yields a second divisor. The third divisor is obtained by reducing the same exponent a second time. This reduction is carried out until the exponent is zero. In this case it cannot be reduced any further. So we restore the original value (the exponent as found for the number n) and reduce the exponent of the second smallest prime factor. After this we start reducing the exponent of the smallest prime factor again. We have to to this, since we have not yet used all the values of the exponent of the smallest prime factor in combination with a reduced exponent of the second smallest prime factor. When the exponent of the smallest prime factor reaches zero again, we restore the original value (just as before) and reduce the exponent of the second smallest prime factor a second time. If the exponent of the second smallest prime factor reaches zero, we restore its original value and reduce the exponent of the third smallest prime factor. We then restart looping through the reductions of the exponents of the smallest and the second smallest prime factor. This is necessary, since all the combinations of exponents of the smallest and the second smallest prime factor have not yet been used in combination with a reduced exponent of the third smallest prime factor. In the same way we handle the other exponents until we reach the situation in which all exponents are zero. This corresponds to d = 1. d = 1 indicates that we generated all possible combinations and thus can abort the loop. Maybe the method is best understood, if you carry it out on a piece of paper for a simple example, e.g. for n = 180 = 2^2 * 3^2 * 5^1. Then the first divisor is d = 90 = 2^1 * 3^2 * 5^1, the second is d = 45 = 2^0 * 3^2 * 5^1, the third is d = 60 = 2^2 * 3^1 * 5^1, and so on. Note that with the above scheme the divisors generated first are larger than those generated later on (although they are not sorted in descending order, this statement holds on average). This leads to the favourable situation that the sum increases fast in the beginning and slowly later on. With this behaviour numbers n with a lot of divisors, i.e. whose sum of divisors exceeds their own value by far, can be filtered out very early, since every time we add a divisor, we check whether the sum of the divisors already exceeds n and abort the function, if this is the case. The only question that remains is: Why should one use such a complicated function, if there is a much simpler solution (i.e. the first version)? The reason is that the prime factor version is much faster than the other. To test the numbers up to 1000000 whether they are perfect, it takes only one sixth of the time with the above function than it does with the first (simpler) version. And this ratio increases for larger numbers. ----------------------------------------------------------------------*/ int main (int argc, char *argv[]) { /* --- find perfect numbers */ int min, max; /* search interval */ int n; /* number to try */ if (argc != 3) { /* if wrong number of arguments given */ printf("usage: %s min max\n", argv[0]); printf("find perfect numbers in the interval [min, max]\n"); return 0; /* print a usage message */ } /* and abort the function */ min = atoi(argv[1]); /* convert lower limit to a number */ max = atoi(argv[2]); /* convert upper limit to a number */ if ((min <= 0) || (max < min)) { /* check for correct arguments */ printf("%s: illegal interval\n", argv[0]); return 1; } for (n = min; n <= max; n++) /* test the numbers in the interval */ if (perfect(n)) /* if a number n is perfect, */ printf("%i\n", n); /* print it */ return 0; /* return 'ok' */ } /* main() */ /*---------------------------------------------------------------------- You should not need much comment on the main function. The command line evaluation is done as usual (cf. e.g. program heron.c). In the loop the function 'perfect' is called for each n and if it returns 'true', the number is printed. (Note that in C any number not equal to zero is interpreted as 'true' and zero as 'false'. Hence you need not write 'if (perfect(n) != 0)', 'if (perfect(n))' is sufficient.) ----------------------------------------------------------------------*/