This article is inspired by the idea of “literate programming,” a technique pioneered by an ancient sorcerer named Donald Knuth. In this paradigm, a program is written as a document that explains the problem being solved. I accomplished this in a fairly naïve way by using a script that walks through the code in-order and generates a markdown file from the comments. Usually, literate programming is done with more sophisticated tools that assemble the code out-of-order so that it can be presented in a manner optimized for human comprehension. I didn’t feel the need to do this, because the program here is short and already pretty coherent as-is.¹ Thus, you could copy the code here directly into a python session and have a working result.

The program in question is the solution to a problem from Project Euler, a website that offers over 600 math problems designed to be solved through programming. I thought long and hard before publishing this article, because when you solve a Project Euler problem, you’re given the message:

We hope that you enjoyed solving this problem. Please do not deprive others of going through the same process by publishing your solution outside of Project Euler

However, this problem has been out for a long time now, and many others have posted their solutions. Since it is already very easy to cheat, I don’t think that posting one additional solution is likely to deprive any additional people of the chance to solve the problem. And if someone is going to cheat, they might as well do it with a solution that carefully explains the problem in detail.

Nevertheless, if you have not solved the problem yourself, I strongly encourage you to try and do it yourself first. Even if you get stuck, I would recommend coming back to it another day before reading my answer.

And just to keep my conscience clear, I have included a mistake in this program. It will give the wrong answer unless you make a small change.

Now, without further ado:

Problem 151

A printing shop runs 16 batches (jobs) every week and each batch requires a sheet of special colour-proofing paper of size A5.

Every Monday morning, the foreman opens a new envelope, containing a large sheet of the special paper with size A1.

He proceeds to cut it in half, thus getting two sheets of size A2. Then he cuts one of them in half to get two sheets of size A3 and so on until he obtains the A5-size sheet needed for the first batch of the week.

All the unused sheets are placed back in the envelope.

At the beginning of each subsequent batch, he takes from the envelope one sheet of paper at random. If it is of size A5, he uses it. If it is larger, he repeats the ‘cut-in-half’ procedure until he has what he needs and any remaining sheets are always placed back in the envelope.

Excluding the first and last batch of the week, find the expected number of times (during each week) that the foreman finds a single sheet of paper in the envelope.

Give your answer rounded to six decimal places using the format x.xxxxxx.

Solution

To solve this problem, we’ll need some way to represent a bag. We could write a Bag class, but that would probably just over-complicate the situation. The bag is really just 5 integers: one for each paper size. Python is already pretty great at working with integers, so we would just be creating more work for ourselves by writing a class.

(Since this is just a math/programming exercise, there’s no need to worry about maintainability, sharability, and such.)

We could instead use an array of integers. For instance, bag[0] would be the count of a1 sheets, bag[1] would be the count of a2 sheets, etc. A bag that contained 1 A2 sheet, 2 A4s and 2 A5s would be written thus:

array_bag = [0, 1, 0, 2, 2]

But there’s actually an even simpler way, which is to just use a single integer. Each digit will contain the number of sheets of paper of the corresponding size. The above bag would be:

int_bag = 1022

(0) - 10000's place - A1
 1  - 1000's place  - A2 
 0  - 100's place   - A3
 2  - 10's place    - A4
 2  - 1's place     - A5

We’ll need a few utility functions to work with numbers in this format.

First, we’ll want to be able to convert a paper size to the digit place in which it is counted, as in the table above. To get the place of the nth digit in a number, you just take 10^n. Since the size numbering is “backwards,” we need to take 10^(5 - size)

def place(size):          # Python uses the ** operator for exponentiation,
    return 10**(5 - size) # because ^ is used for bitwise XOR, as in C.

To get the number of sheets of a particular size, we need to strip away the other digits. For instance, sheets_of_size(3, 1211) should return 2.

def sheets_of_size(n, bag):
    r = bag       # Start with the number of the bag.         eg. 1211
    p = place(n)  # Get the place of the number                    100
    r = r / p     # Divide the bag number by the place number       12.11
    r = trunc(r)  # Truncate the digits after the decimal point     12
    return r % 10 # Remove the preceding digits by taking % 10       2

One of the main advantages of representing bags as integers is that it’s easy to add and subtract sheets: we just need to add an integer. To add an a3, and a4, and an a5, instead of writing:

array_bag[2] += 1
array_bag[3] += 1
array_bag[4] += 1

we can just write:

int_bag += 111

Indeed, we’ll have to do exactly this to represent the foreman cutting a piece of paper from the bag. Cutting a paper always yields one of each smaller paper size, including A5, since the foreman uses up one of them.

Let’s write a function to do this. It will take a bag and the paper size to be cut, and return the resulting bag.

eg. cut(1, 10000) → 1111
    cut(2,  1011) →  122
    cut(3,  1210) → 1121
    cut(5,  1101) → 1100

def cut(size, bag):
    p = place(size) 
    return int(
          bag            # Take the original bag         eg.  1210
        - p              # Remove the sheet                 - 100
        + ((p - 1) / 9)  # Add one of each smaller sheet    + (99 / 9)
    ) 

To solve the problem, we’ll need to test whether a bag is has only one sheet remaining. We could calculate this mathematically fairly easily, but it would be slower and require more code then just listing each bag.

def only_one(bag):
    return (bag == 1   or bag == 10   or 
            bag == 100 or bag == 1000 or bag == 10000)

We’ll also want to be able to get the total number of sheets in the bag, for reasons that will soon become clear.

def total_sheets(bag):
    return (
          sheets_of_size(1, bag)
        + sheets_of_size(2, bag)
        + sheets_of_size(3, bag)
        + sheets_of_size(4, bag)
        + sheets_of_size(5, bag)
    )

Now that we’ve laid the groundwork, let’s get into exactly how we’ll solve the problem. The obvious solution is to check every possible outcome of each batch, count the number where the bag has only one sheet in it, and find the average.

This is the right path, but we also have to remember that not every possibility is equally likely. For instance, if the foreman starts the batch with the bag 0120 (1 A3, 2 A2s), he’s twice as likely to cut an A2 as an A3. Consequently, we need to get every possibility paired with its chance of occurring. However, the chance of getting a bag depends on the chance of getting the bag that came before it. Consequently, we need to descend a tree of possibilities. Four batches down, it will look like this:

10000 (100%)
    1111 (100%)
        222 (25%)
            133  (8.33%)
            213  (8.33%)
            221  (8.33%)
        1022 (25%)
            233  (5%)
            1013 (10%)
            1021 (10%)
        1101 (25%)
            212  (8.33%)
            1012 (8.33%)
            1100 (8.33%)
        1110 (25%)
            221  (8.33%)
            1021 (8.33%)
            1101 (8.33%)

At each batch, there are more and more possibilities, and the probability of one of them occurring falls each time. On each level, the probabilities add up to 100%, because some bag is guaranteed to result from each batch the foreman runs.

To get all the probabilities, we’ll need write a function that will take a bag and its probability of occurring, and return a list of the possible bags after the foreman runs the batch, and their probabilities of occurring.

eg. possibilities([1022, .25]) → [[ 233, .05], 
                            
   [1013, .1 ], 
                            
   [1021, .1 ]]

We’ll use proportions instead of percents to represent probabilities. A proportion is the same as a percent chance divided by 100. Basically, instead of “per 100” it’s “per 1.” If this is confusing, one way to think of it is that if something happens one half of the time, it happens .5 of the time.

def possibilities(poss):
    bag = poss[0]  # We're using an array instead of two arguments so that
    prob = poss[1] # we can pass items in the list returned by this function
                   # as arguments to the function.
        
    p = [] # The possibilities will be stored in a list.
    for size in range(1, 6): # For each paper size
        num_of_size = sheets_of_size(size, bag)
        if num_of_size > 0:
            p.append([
                # The `cut` function that we wrote earlier returns the
                # resulting bag if we cut a sheet of `size`.
                cut(size, bag),
                
                # Multiplying the probabilities of two events gives the
                # probability that both will occur. Since a result can only
                # occur if the result that produced it also occurred, we need 
                # to multiply the probability of both events together.
                # its probability 
                (num_of_size / total_sheets(bag)) * prob
            ])

    return p

Now, using the above function, we need to recursively traverse the probability tree and add up the probabilities of the cases where there is only one sheet in the bag.

We’ll use a function that will take two parameters. One, poss will be a possibility, in the same format as above: [bag, probability]. The next, n, will be the number of recursions that should be used (or the depth that the possibility tree should be descended.

def expected(poss, n):
    bag = poss[0]
    prob = poss[1]

    # If this result contains only one sheet, augment the expected value 
    # by the probability of getting this result.
    p_sum = int(only_one(bag)) * prob

    # If we've reached the maximum number of recursions, don't explore
    # the possibilities which result from this one.
    if n == 0: return p_sum
    
    # For each resulting possibility, augment the expected value by 
    # the output of this function applies to it, reducing `n` to 
    # prevent infinite recursion.
    for e in possibilities(poss):
        p_sum += expected(e, n - 1)

    return p_sum

Now we’ve done everything we need to solve the problem! We just need to be sure to pass the right parameters to the functions we wrote.

Because the foreman starts every week with an A5 sheet (10000), the probability is 100% (1).

start = [10000, 1]

The answer has to exclude the first and last of the 16 batches, so we will go 14 iterations deep rather than 15.

solution = expected(start, 14)

And finally, we output our answer rounded to 6 decimal places!

print(round(solution, 6))