Consider the following general class of urn and marble problems. Given 1, 2, ... m numbered immovable urns and a collection of n identical marbles, how many different arrangements of the marbles in the urns are possible by allowing up to a maximum of p marbles in each urn? Note that any number of the n marbles can be used from 0 up to all n of them.
Specifically, if you have 64 different urns and 24 identical marbles and the rule that you can have from 0 to a maximum of 4 marbles in an urn, how many unique arrangements are possible? It appears that no easy direct computational formula exists to solve this problem. If only up to one marble is allowed in each urn, a sum of the binomial coefficients _{64}C_{x} from x = 0 to 24 would yield one possible formula. Allowing two or more marbles in each urn creates more computational difficulty. So, this problem may benefit from an appropriately programmed computer. In the phrases of Nicolas Wirth, this problem may benefit from appropriate algorithms and data structures.
How should the computer be programmed? Using C allows great speed in execution. As we see below, speed and an effective algorithm are not enough to solve this problem. Other languages, which are descendants of C, have built in data structures that allow a problem solver to focus on development of a more efficient algorithm.
There are 5^{m} arrangements to check, where m is the number of urns. Tony's original algorithm decreased the number of arrangements to be counted by skipping all the arrangements where more than 24 marbles were used. Consider the following program that has some changes to streamline the inner counting loop. For modularity during development, add_marble and count_marbles were separate functions in the original program. Code from these functions was placed in-line in the main function for speed considerations. This C program keeps Tony's original arrangement counter algorithm, which handled quite long integers nicely by creating a 100 digit, base 1 billion counter. The running_count variable was introduced to reduce the number of times the number of marbles in use was tallied.
// urn3.c /* A few revisions for speed: 1. add_marble and count_marbles both count the marbles: discard one. 2. Passing constants in a subroutine call is a waste of time. 3. In-line subroutines when memory is not an issue and particularly when you only call them once. 4. Use break; and forget structured programming. 5. Don't use dynamic memory allocation unless memory is an issue. 6. Try counting marbles less often. */ #include<stdio.h> #include<stdlib.h> #include<time.h> int main (int argc,char **argv) { int bctr[100]; // base billion counter with 100 "digits". int temp,running_count,i,b,t; int m,n,p; // m is number of urns. n is number of marbles. p is max in urn. int urns[101]; // array to represent the urn/marble states possible. time_t start,finish; if (argc != 4) { printf("\nUsage: urn3 urns marbles max-per-urn\n"); exit(1); } m = atoi(argv[1]); n = atoi(argv[2]); p = atoi(argv[3]); m = m < 101 && m > 0 ? m : 64; n = n > 0 ? n : 24; p = p > 0 && p >= n/m ? p : n/m; // don't bother with the easy cases for (i=0;i < m+1;i++) urns[i]=0; // Initialize urns to zero. for (b=0;b < 100;b++) bctr[b]=0; // Initialize the base billion counter. start = time(0); running_count = 0; while (urns[m] == 0) { urns[0]++; // Always add 1 marble to rightmost urn. running_count++; for(t=0;t < m;t++) { if(urns[t] > p) { // If an urn exceeds the max running_count -= p - 1; // adjust running count of marbles urns[t]=0; // empty urn urns[t+1]++; // and add 1 marble to urn on left. } else break; // exit when stops "carrying". } if (running_count > n) { // Add elements in urns from the left to the right. When sum exceeds n+1, // fill the rest of the urns to the right with p marbles each. temp=0; t = m - 1; // Start at leftmost urn. while (temp < n+1 && t >= 0) { temp+=urns[t]; t--; } running_count = temp; for (;t >= 0;t--) { // fill any uncounted urns running_count += p - urns[t]; urns[t]=p; } } if(running_count <= n) { bctr[0]++; for(b = 0;b < 100;b++) { // carry if(bctr[b] >= 1000000000) { bctr[b]=0; bctr[b+1]++; } else { break; } } } } finish = time(0); printf("For %d urns and %d marbles with a max of %d in each urn,\n", m,n,p); printf("there were "); b=99; while (bctr[b] == 0 && b)b--; // skip leading billions printf("%ld", bctr[b--]); while (b >= 0) printf("%09ld", bctr[b--]); printf(" arrangements counted in %ld seconds.\n", finish - start); return 0; }produced
Starting with 18 urns: Sun Feb 27 11:54:00 EST 2005 Results with 18 urns: Sun Feb 27 23:36:37 EST 2005 For 18 urns and 24 marbles with a max of 4 in each urn, there were 103917328350 arrangements counted in 42157 seconds. Starting with 19 urns: Mon Feb 28 08:41:00 EST 2005 Results with 19 urns: Tue Mar 1 18:59:02 EST 2005 For 19 urns and 24 marbles with a max of 4 in each urn, there were 262270328730 arrangements counted in 123482 seconds.
Even with the changes, which improved the speed by about a factor of 6, this C program took almost 36 hours to run for just m=19 urns. So, how far away was the finish line of 64 urns using this C program? So far we have:
Urns | Arrangements | Seconds to calculate | Arrangements per second |
---|---|---|---|
13 | 470597480 | 69 | 6820253 |
14 | 1562441800 | 298 | 5243093 |
15 | 4855374080 | 1149 | 4225738 |
16 | 14208711350 | 4127 | 3442866 |
17 | 39381411950 | 13668 | 2881285 |
18 | 103917328350 | 42157 | 2465007 |
19 | 262270328730 | 123482 | 2123955 |
>>> m = 20 >>> seconds = 123482 >>> while m < 65: ... seconds *= 3 ... m += 1 ... >>> seconds 364804441630310046775834926L >>> years = seconds / (3600*24*365) >>> years 11567872958850521523LIt looks like we need a better algorithm! Even to get to 24 urns it would take over 4 months of computing with the hardware at hand. (In retrospect, when we notice the rate of increase goes down after there are more urns than marbles, 10^{19} years is a gross overestimate. However, even 30 million years, to get to 64 urns this way, is still prohibitive.
What next? Can we use anything from those classes in Discrete
Math, Theory of Computation and Data Structures? If we make
the observation that the set of all arrangements can be
partitioned into subsets depending on the number of marbles in
urn_{x}, then the cardinality of the whole set is the
sum of the cardinalities of the 5 subsets.
arrangements(m urns, p marbles) = | arrangements(m-1 urns, p-0 marbles) + |
arrangements(m-1 urns, p-1 marbles) + | |
arrangements(m-1 urns, p-2 marbles) + | |
arrangements(m-1 urns, p-3 marbles) + | |
arrangements(m-1 urns, p-4 marbles) |
Now we have a tree of height m enumerating our arrangements. That still leaves 5^{m-1} nodes to count. A little smaller number (and a bad pun). To avoid calculating any node more than once, an array of intermediate results would be helpful: arrangements_{m,p} The arrays in Python, Perl and Ruby are not the rectangular arrays of C, but rather nested arrays. The access is still array[row][column], but the structure is more similar to char **array in C. These arrays can contain heterogenous objects like integers of widely varying length.
The facility to handle really long integers would also be helpful. In Perl it is necessary to include the Math::BigInt module. Before about Python version 2.3, you need to specifically use long integers. In Ruby the transition from regular to long integers is automatic.
The following programs in Python, Perl and Ruby show that each has some advantages for making this problem more manageable. Perl is the easiest language in which to set up our doubly subscripted array of intermediate values. Perl ran all 64 urns in about 1.2 seconds. Python is the most cumbersome language in which to set up the array, but its long integer implementation was about 15 times faster than Perl's. Setting up the array in Ruby is just slightly harder than in Perl and easier than in Python. Ruby's long integers are just as fast as Python's.
m, the number of urns, n, the number of marbles and p, the maximum number of marbles per urn, are entered as arguments on the command line. p is a global constant. Our formula has be implemented as a recursively called function. Notice that m has been called x in that function and n has been called y. This is to draw attention to the somewhat subtle differences in the rules for variable scope in Perl, Python and Ruby.
Finally, to make our long integers more readable, we add some commas. Perl has the oldest and still most complete support for regular expressions, including zero width assertions. Python allows look behind assertions of fixed length. Ruby only allows look ahead assertions. The effect remains the same: Scan the integer from left to right, squeezing in a comma everywhere there is digit to the left and groups of 3 digits to the right. This is not the most efficient algorithm possible. Since displaying the results is not part of the inner loops, an effective algorithm is fine.
#!/usr/bin/python # urn.py import sys,time,re # Main subroutine calculates a tree from its root, out to the leaves. def node(x,y): global a sum = 0L for i in range(min(y,p),-1,-1): if a[x - 1][y - i] > 0: # If a sub-tree was already calculated then sum += a[x - 1][y - i] # do not recalculate it. elif x == 1: # Level 1 is always a leaf node. sum += 1 else: sum += node(x - 1,y - i) a[x][y] = sum # Save each newly calculated sub-tree. return sum # Get the parameters from the command line. if len(sys.argv) != 4: # Give up if no parameters. print "Usage: urn.py urns marbles max-per-urn" sys.exit(1) m = int(sys.argv[1]) n = int(sys.argv[2]) p = int(sys.argv[3]) # Initialize the intermediate results array (0 marbles -> 1 arrangement) a = [] for i in range(m + 1): # In Python range(5) is range(0,5) or 0,1,2,3,4 b = [] # (like C++ STL iterators). b.append(1L) for j in range(1,n + 1): b.append(0L) a.append(b[:]) # b[:] is a new copy of b, not another reference to b # Run start = time.time() for i in range(1,m + 1): ans = re.sub(r'(?<=\d)(?=(\d\d\d)+$)',',',str(node(i,n))) # Add commas print "%2d %d %d %30s" % (i, n, p, ans) # for readability. print "Calculation took %d seconds" % (time.time() - start)produced
jm@cbw:~/ltu/mcs$ time ~/berard/urn.py 64 24 4 1 24 4 5 2 24 4 25 3 24 4 125 4 24 4 625 5 24 4 3,125 6 24 4 15,625 7 24 4 78,005 ... 18 24 4 103,917,328,350 19 24 4 262,270,328,730 ... 63 24 4 1,590,937,630,234,684,879,080 64 24 4 2,194,493,026,399,270,064,325 Calculation took 0 seconds real 0m0.442s user 0m0.070s sys 0m0.010s
#!/usr/bin/perl -w # urn.pl use Math::BigInt; # Main subroutine calculates a tree from its root, out to the leaves. sub node { my ($x,$y) = @_; my $z = $p <= $y ? $p : $y; my $sum = Math::BigInt->new(0); for my $i (reverse (0..$z)) { if ($a[$x - 1][$y - $i] > 0){ # If a sub-tree was already calculated then $sum += $a[$x - 1][$y - $i]; # do not recalculate it. } elsif ($x == 1) { # Level 1 is always a leaf node. $sum += 1; } else { $sum += node($x - 1,$y - $i); } } $a[$x][$y] = $sum; # Save each newly calculated sub-tree. return $sum; } # Get the parameters from the command line. die "Usage: urn.pl urns marbles max-per-urn" unless (scalar @ARGV == 3); # Give up if no parameters. $m = $ARGV[0]; $n = $ARGV[1]; $p = $ARGV[2]; # Initialize the intermediate results array (0 marbles -> 1 arrangement) for my $i (0..$m) { $a[$i][0] = Math::BigInt->new(1); for my $j (1..$n) { $a[$i][$j] = Math::BigInt->new(0); } } # Run. $start = time(); for my $i (1..$m) { ($ans = substr(node($i,$n),1)) =~ s/(?<=\d)(?=(\d\d\d)+$)/,/g; # Add commas printf "%2d %d %d %30s\n",$i,$n,$p,$ans; # and remove + for readability. } $dif = time() - $start; $plural = $dif != 1 ? "s" : ""; print "Calculation took $dif second$plural\n"; # Display the intermediate results table for small n if ($n < 10) { for my $i (0..$m) { for my $j (0..$n) { printf "%4d ",$a[$i][$j]; } print "\n"; } }produced
jm@cbw:~/ltu/mcs$ time ~/berard/urn.pl 64 24 4 1 24 4 5 2 24 4 25 ... 63 24 4 1,590,937,630,234,684,879,080 64 24 4 2,194,493,026,399,270,064,325 Calculation took 1 second real 0m1.224s user 0m1.130s sys 0m0.000s
#!/usr/bin/ruby # urn.rb # Main subroutine calculates a tree from its root, out to the leaves. def node(x,y) sum = 0 Array(0..Array[y,$p].min).reverse_each do |i| if ($a[x - 1][y - i] > 0) # If a sub-tree was already calculated then sum += $a[x - 1][y - i] # do not recalculate it. elsif (x == 1) # Level 1 is always a leaf node. sum += 1 else sum += node(x - 1,y - i) end end $a[x][y] = sum # Save each newly calculated sub-tree. return sum; end # Get the parameters from the command line. unless (ARGV.length == 3) # Give up if no parameters. print "Usage: urn.rb urns marbles max-per-urn\n" exit(1) end m = ARGV[0].to_i n = ARGV[1].to_i $p = ARGV[2].to_i # Initialize the intermediate results array (0 marbles -> 1 arrangement) $a = Array.new Array(0..m).each do |i| $a[i] = [1] + [0] * n end # Run. start = Time.now Array(1..m).each do |i| ans = node(i,n).to_s ans.gsub!(/(\d)(?=(\d\d\d)+$)/,'\1,') # Add commas printf "%2d %d %d %30s\n",i,n,$p,ans # for readability. end printf "Calculation took %d seconds\n",Time.now - startproduced
jm@cbw:~/ltu/mcs$ time ~/berard/urn.rb 64 24 4 1 24 4 5 2 24 4 25 ... 63 24 4 1,590,937,630,234,684,879,080 64 24 4 2,194,493,026,399,270,064,325 Calculation took 0 seconds real 0m0.356s user 0m0.080s sys 0m0.000s
The developers of Python, Perl and Ruby all have tried to make, "the difficult possible." Hopefully these examples will help foster an appreciation of the work of Guido van Rossum, Larry Wall, Yukihiro Matsumoto and their many friends.
Our improved algorithm could certainly have been written in a number of languages, including the language of the original algorithm: C. To paraphrase the ideas of Matz, the creator of Ruby, about the feel of a perfect language: Many languages have similar power. The "perfect" language for an individual programmer depends on the way that programmer approaches the task of writing code. The likelihood of the existence of a perfect language for every programmer and every task is vanishingly small.
These examples are too short to allow you to pick a favorite programming language. You will have to try them out yourself. You might find you like them all.
Revised March 7, 2005