2005-02-03 C(n,r) introduces dynamic programming
A simple problem: Calculate the value of the combinatoric function, exactly.
C(n,r) is the number of different combinations of r items that may be taken from a collection of n.
Or the number of r-element subsets that exist for any n-element set.
The best-known definition is .
Let’s make it easy. This is not a problem about implementing arithmetic on giant integers. It is guaranteed that the result of the computation will fit in a normal integer. So what about this:
int factorial(int x)
{ int result=1;
for (int i=1; i<=x; i+=1)
result*=i;
return result; }
int C(int n, int r)
{ return factorial(n)/(factorial(r)*factorial(n-r)); }
Not very good.
What if n=20 and r=19? C(20,19) is exactly 20, and that fits perfectly in any int variable, but calculating C(20,19) involves calculating factorial(20), which is huge and certainly causes overflow. So the computation will fail even though the result is nice and small.
void main(void)
{ cout < "C(20,19) = " < C(20,19) < "\n"; }
$ CC comb1.cpp
$ a.out
C(20,19) = -19
We can do a simple transformation, and discover that C(n,r) = (n/r)*C(n-1,r-1) is also true.
Also not very good.
All values of C(4,3) are obviously whole numbers, but calculating C(4,3) this way will involve calculating 4/3. Done as an int, the result is 1 which is just plain wrong. Done as a float or double, the result is 1.3333333333, which can’t possibly have enough digits to be truly accurate: an error has slipped in at the very first step.
int C(int n, int r)
{ if (n==r || r==0) return 1;
return n/r * C(n-1,r-1); }
void main(void)
{ cout < "C(20,19) = " < C(20,19) < "\n"; }
$ CC comb2.cpp
$ a.out
C(20,19) = 2
Another, slightly trickier transformation reveals that
C(n,r) = C(n-1,r-1) + C(n-1, r)
#include <iostream>
#include <stdlib.h>
#include <sys/time.h>
#include <sys/resource.h>
double getcputime(void)
{ struct timeval tim;
struct rusage ru;
getrusage(RUSAGE_SELF, &ru);
tim=ru.ru_utime;
double t=(double)tim.tv_sec + (double)tim.tv_usec / 1000000.0;
tim=ru.ru_stime;
t+=(double)tim.tv_sec + (double)tim.tv_usec / 1000000.0;
return t; }
int C(int n, int r)
{ if (n==r || r==0) return 1;
return C(n-1,r-1) + C(n-1,r); }
void main(int argc, char * argv[])
{ int n=atol(argv[1]), r=atol(argv[2]);
double t0=getcputime();
int ans=C(n, r);
double t1=getcputime();
cout < "C(" < n < "," < r < ") = " < C(n,r) <
" and took " < (t1-t0) < " S to compute\n"; }
$ CC comb4.cpp
$ a.out 20 19
C(20,19) = 20 and took 1.2e-05 S to compute
$ a.out 20 10
C(20,10) = 184756 and took 0.022407 S to compute
$ a.out 24 12
C(24,12) = 2704156 and took 0.221516 S to compute
$ a.out 28 14
C(28,14) = 40116600 and took 3.30662 S to compute
Using , and of course , the following two derivations are valid:
1) ______
=
=
=
/ = /2) ______
=
=
=
=
=
=
=
/ = /______
Base cases are clearly n=r and r=0.
A memo function
$ cat comb5.cpp
#include <iostream>
#include <stdlib.h>
static int Canswers[30][30] = { -1 };
int C(int n, int r)
{ if (n==r || r==0) return 1;
if (Canswers[n][r] == -1)
Canswers[n][r] = C(n-1,r-1) + C(n-1,r);
return Canswers[n][r]; }
void main(int argc, char * argv[])
{ if (argc!=3)
{ cerr < "Need n and r on command line\n";
exit(1); }
int n=atol(argv[1]), r=atol(argv[2]);
cout < "C(" < n < "," < r < ") = " < C(n,r) < "\n"; }
$ CC comb5.cpp
$ a.out 20 10
C(20,10) = 0
What the ...?
Stupid C++, that’s what.
#include <iostream>
#include <stdlib.h>
static int Canswers[30][30];
void init(void)
{ for (int i=0; i<30; i+=1)
for (int j=0; j<30; j+=1)
Canswers[i][j] = -1; }
int C(int n, int r)
{ if (n==r || r==0) return 1;
if (Canswers[n][r] == -1)
Canswers[n][r] = C(n-1,r-1) + C(n-1,r);
return Canswers[n][r]; }
void main(int argc, char * argv[])
{ if (argc!=3)
{ cerr < "Need n and r on command line\n";
exit(1); }
int n=atol(argv[1]), r=atol(argv[2]);
init();
cout < "C(" < n < "," < r < ") = " < C(n,r) < "\n"; }
The initialisation rules in C++ are not what people usually think they are.
Why does this design help?
The original program based on C(n,r) = C(n-1,r-1) + C(n-1, r) wasted a lot of time answering the same questions over and over again. To work out C(20,10) you have to work out C(19,9) and C(19,10). T work out C(19,10) you have to work out C(18,9) and C(19,9): four computations, but two are the same. And it only gets worse:
$ cat comb7.cpp
#include <iostream>
#include <stdlib.h>
int count[25][25];
int C(int n, int r)
{ count[n][r]+=1;
if (n==r || r==0) return 1;
return C(n-1,r-1) + C(n-1,r); }
void main(int argc, char * argv[])
{ int n=atol(argv[1]), r=atol(argv[2]);
int ans=C(n, r);
cout < "C(" < n < "," < r < ") = " < C(n,r) < "\n";
for (n=1; n<20; n+=1)
{ for (int r=1; r<=n; r+=1)
cout < "C(" < n < "," < r < ")x" < count[n][r] < ", ";
cout < "\n"; } }
$ CC comb7.cpp
$ a.out 20 10
C(20,10) = 184756
C(1,1)x97240,
C(2,1)x97240, C(2,2)x48620,
C(3,1)x48620, C(3,2)x48620, C(3,3)x22880,
C(4,1)x22880, C(4,2)x25740, C(4,3)x22880, C(4,4)x10010,
C(5,1)x10010, C(5,2)x12870, C(5,3)x12870, C(5,4)x10010, C(5,5)x4004,
C(6,1)x4004, C(6,2)x6006, C(6,3)x6864, C(6,4)x6006, C(6,5)x4004, C(6,6)x1430,
...
This shows that in computing C(20,10), C(1,1) is computed 97,240 times. What a waste! It isn’t going to change from one computation to the next, so why not just work it out once, then remember it?
In fact, why not just work out the answers to all the C(...,...)-type questions before the function is ever called? Then the function can be replaced by a simple array look-up.
The recursive function works by propagating demand: you’re asked a question, and have to ask other questions before you can work out the answer. The alternative is to propagate information: we already know the answers to the simple (base) cases; from them we can produce the answers to slightly bigger questions, and from them we can produce the answers to even bigger question yet, before the questions are ever actually asked.
$ cat comb9.cpp
#include <iostream>
#include <stdlib.h>
const int max=30;
int Canswers[max][max];
void prepare(void)
{ for (int n=0; n<max; n+=1)
{ Canswers[n][0]=1;
for (int r=1; r<n; r+=1)
Canswers[n][r] = Canswers[n-1][r-1] + Canswers[n-1][r];
Canswers[n][n]=1; } }
int C(int n, int r)
{ return Canswers[n][r]; }
void main(int argc, char * argv[])
{ prepare();
int n=atol(argv[1]), r=atol(argv[2]);
cout < "C(" < n < "," < r < ") = " < C(n,r) < "\n"; }
$ CC comb9.cpp
$ a.out 20 10
C(20,10) = 184756
This technique only works if there is some good order for filling in the table. We are guaranteed that row (n-1) will be filled before we start working on row n, and all the computations for row n depend only on values in row (n-1).
This is really what dynamic programming is about: you have a naturally recursive solution to a problem, but that solution may involve really wasteful recomputation of the same sub-results over and over again. If you can tabulate the results in some rational order in advance, all recomputation can be avoided. If you can’t, then maybe a memo function is still a possibility.
The weekend brain-teaser:
How do you find the Minimum-Sum Contiguous Subsequence of a list of integers, quickly?
For example, given this list:
{ 17, -3, -4, 5, -2, 9, -5, 5, 3, 2, -4, -6, 1, -3, 2, -1, 3, -6, 8, -2, 4, -1 }
the subsequence { -3, -4 } has a sum of -7, but the subsequence { -4, -6 } has a sum of -10 which is better, and { -4, -6, 1, -3 } has a sum of -12 which is even better. Which is the absolute best?
“Contiguous” means that the sub-sequence can start and stop anywhere, but may not skip values, so the list { -3, -4, -2, -5 } is not a valid choice.
The brute force solution is to try out every possible contiguous subsequence, finding all their sums, and remembering the minimum. This is a cubic algorithm. Can you do better?