/*---------------------------------------------------------------------- File : choose.c Contents: compute n choose k for given n and k Author : Christian Borgelt History : 29.04.1996 file created 06.11.1997 some comments added 19.11.1997 comments completed ----------------------------------------------------------------------*/ #include #include #include /*---------------------------------------------------------------------- Functions ----------------------------------------------------------------------*/ int choose (int n, int k) { /* --- compute n choose k */ int i; /* loop variable */ int r = 1; /* result of n choose k */ if (k > n) return 0; /* check range of k */ for (i = 1; i <= k; i++, n--) { if (INT_MAX /n < r) /* if r * n is out of range, */ return -1; /* abort the function */ r = (r *n) /i; /* compute \prod_{i=1}^k (n-i+1)/i */ } /* = n!/(k!(n-k)!) */ return r; /* return result */ } /* choose() */ /*---------------------------------------------------------------------- n choose k, also known as `binomial coefficient', is defined as n!/(k!(n-k)!). The above function computes the value of n choose k without computing the factorials. The idea underlying this is that the fraction n!/(k!(n-k)!) can be reduced to (\prod_{i = n-k+1}^n i)/k!. With this reduction the loop should be clear. What may cause problems is the following: You already know from the program `euler.c' that a quotient whose nominator and denominator are both integer numbers is evaluated as an integer number, i.e. any fraction is discarded from the result. Hence the question may arise, why we can do all computations using only integer numbers. The reason is that, strange enough, none of the quotients computed above has a fraction, so there is nothing to discard. To understand this, look at the first quotients that are computed. In the first execution of the loop it is i = r = 1 and thus n/1 is computed, which obviously does not have a fractional part. The second step (i = 2) yields n(n-1)/2 (where n is the original n, not the one modified by the loop). This cannot have a fractional part, because either n or n-1 must be even. Not both of two consecutive number can be odd. Hence n(n-1) must be divisible by two. Now look at the next step (i = 3). The quotient computed is then n(n-1)(n-2)/(2*3). We already know that the division by 2 cannot lead to a fractional part, since either n or (n-1) must be even. For an analogous reason, either n or n-1 or n-2 must be divisible by 3. To understand this consider the remainders of n, if it is divided by three. This remainder can be 0, 1 or 2. If it is 0, n is divisible by 3. If it is 1, then n-1 must be divisible by 3, since n = m*3 +1 (for some m) because of the remainder 1. If the remainder is 2, then n-2 must be divisible by 3, since n = m*3 +2 (for some m) because of the remainder 2. The next step (i = 4, which will be the last we consider here) is to compute the quotient n(n-1)(n-2)(n-3)/(2*3*4). Among the four number n, n-1, n-2, and n-3 there must be one that is divisible by 4 and two that are divisible by two. Hence n(n-1)(n-2)(n-3) must be divisible by 8. And from the above argument concerning n(n-1)(n-2) it follows that it must also be divisible by 3. Hence n(n-1)(n-2)(n-3)/(2*3*4) does not have a fractional part. The scheme underlying the reasoning should be clear by now und thus we can establish the result that no quotient computed in the function above has a fractional part. Hence integer computations are sufficient. What remains to be explained is the test `if (INT_MAX /n < r)'. INT_MAX is a preprocessor constant that is defined in the header file `limits.h'. Its value is the largest number that can be represented as an `int'. In the same manner there is a constant `FLT_MAX', which is the largest number representable as a 'float', and so on. Take a look at `/usr/include/limits.h', if you are interested. With the above test we make sure that the result of `r*n' does not exceed the largest number that can be represented as an `int'. If we did not check this and `r *n' exceeded INT_MAX, the function would yield a wrong result. It is important to note that you cannot write the above condition as `if (n *r > INT_MAX)', since INT_MAX is the largest number that can be represented as an `int' and even if `n *r' were greater than INT_MAX, the condition `(n *r > INT_MAX)' would be false, since `n *r' is computed as an integer number and no integer number can be greater than INT_MAX. The condition used in the function above is safe, since for n >= 1, it is INT_MAX / n <= INT_MAX and thus must be representable. If the test fails, that is, if `n *r' exceeds the largest integer number, a negative number is returned to indicate that no result can be computed. ----------------------------------------------------------------------*/ int main (int argc, char *argv[]) { /* --- main function */ int n, k; /* arguments */ int r; /* result (n choose k) */ if (argc != 3) { /* if wrong number of arguments given */ printf("usage: %s n k\n", argv[0]); printf("compute n choose k, i.e. n!/(k!(n-k)!)\n"); return 0; /* print a usage message */ } /* and abort the program */ n = atoi(argv[1]); /* convert n to a number */ k = atoi(argv[2]); /* convert k to a number */ if ((n < 0) || (k < 0)) { /* check range of n and k */ printf("%s: argument out of range\n", argv[0]); return 1; } r = choose(n, k); /* compute n choose k */ if (r < 0) { /* if the computation failed, abort */ printf("%s: result cannot be computed\n", argv[0]); return 1; } printf("%i\n", r); /* print result */ return 0; /* return 'ok' */ } /* main() */ /*---------------------------------------------------------------------- Since the function `choose' can fail (see above), we have to check the result before we print it. In the function `choose' we return a negative number to indicate that the computation failed. This is safe, since n choose k is always greater than or equal to 0. Hence, if the function `choose' returns a negative number, we have to print an error message. ----------------------------------------------------------------------*/