Introduction

The Monty Hall problem is a relatively well-known puzzle based loosely on the game show “Let’s Make a Deal.” It’s known for its ability to confuse, and the winning strategy to the puzzle may not seem intuitive at first. Here, I want to describe and prove the winning strategy (with Bayes’ Theorem) and model the puzzle with Python to show that the described strategy actually wins the puzzle with the described probability.

This blog post was featured in Webucator’s Python training. Specifically, here is the video.

Problem Description

You are on a game show and there are three doors in front of you, as in the illustration below. Behind one of the doors sits a prize, but the door has been chosen at random. You’re asked to select a door, and if you select correctly you will win the prize. So, it should seem obvious at this point that, all things being equal and random, there is a 1/3 probability of you selecting the correct door.

But, there’s a catch: When you announce which door you have chosen, the game show host opens one of the other doors. That is, he opens a door that is neither the door you selected nor the door containing the prize. He then asks if you’d like to change your choice of door. There are only two doors remaining – the one you selected and another closed door. Should you stick with your original choice of door, or should you switch your choice to the only other closed door? Does it matter?

A Winning Strategy

Yes, it matters. Always change your choice of door. Always. Simply put, if you stick with your original door selection, you will only win one out of every three times you play the game. If you change your selection, you will win two out of every three times you play.

Proving The Strategy

Intuitively, it is pretty straightforward that your original choice (when there are three doors to choose from) has a 1/3 probability of winning the game. Where things become complicated is when a door is removed. Many would argue that you now have a 1/2 probability of winning by selecting either door, but this isn’t the case.

One critical aspect that you must realize is the following: The game show host can select a door to remove because he knows (a) which door you selected and (b) which door has the prize. The act of removing a door from the game is not random. In fact, in many cases, the host must remove a specific door. For instance, if you select door 1 and the prize is behind door 3, the host has no other option than to remove door 2. You must realize that the choice of which door to remove is conditional on both the prize door and the door you initially select. The only “new” information that a removed door gives you is that it tells you where the prize isn’t. How can you use that information?

Think of the situation this way: When you made your original door selection, there was an equal and random 1/3 probability that the prize was behind any particular door. That situation looks like the following picture, where we’ve selected door 1 as an example.

Put another way, there is a 1/3 probability that door 1 contains the prize and a 2/3 probability that the prize is behind one of the two other doors. That looks like this picture.

So, what happens when we open door 2 or 3? Actually, nothing. The probability doesn’t magically change – it doesn’t change at all. There’s still a 1/3 probability that the prize is behind door 1 and a 2/3 probability that the prize is not behind door 1. The only new information we have is about where the prize is not – as removing a door is completely dependent on the host knowing where the prize isn’t. So, now we have a situation that looks like this.

Of course, the fact that door 2 was opened in this example shows us that there is a 0% probability that the prize is behind door 2. The implication, as strange as it feels, is that there is a 2/3 probability that the prize is behind door 3… because the combined probability of the prize being behind doors 2 and 3 is 2/3 and we know that the prize isn’t behind door 2. Clearly, you should change your selected door to door 3 in this situation.

Rigorous Proof Using Bayes’ Theorem

Bayes’ Theorem allows us to compute conditional probabilities. That is, Bayes’ Theorem will allow us to compute a probability of the form $P(A|B)$, read as “the probability of event A given that event B holds.” We will use Bayes’ Theorem here to compute the probability of winning the game based on switching or not switching our choice in door once a door is removed.

Without loss of generality, we will select door 1 for our initial door choice. Also without loss of generality, we will say that door 2 is opened (as in the above illustrations), leaving us with doors 1 and 3 as possible prize doors. We define the following events: $D_1$, $D_2$, and $D_3$ are the events that the prize is behind doors 1, 2, and 3, respectively. $O_2$ is the event of the host opening door 2, which is dependent on two things: (a) our initial choice of door 1 and (b) the host’s knowledge that the prize isn’t behind door 2.

We wish to compute two things. The first is formally written as $P(D_1|O_2)$. That is, we wish to compute the probability of the prize being behind door 1 given that the host has removed door 2. This represents the probability of winning the game when we do not change door selection. By Bayes’ Theorem, we have
$P(D_1|O_2)=\frac{P(O_2|D_1)P(D_1)}{P(O_2)}.$
Two of the three components on the right hand side are relatively straightforward. Namely, $P(D_1)$ is the probability that the prize is behind door 1, which is 1/3. The conditional probability $P(O_2|D_1)$ is literally read as “the probability that door 2 is removed when the prize is located behind door 1.” In such a situation, the host may remove either door 2 or 3, and so this is equal to 1/2. The denominator is more challenging here. $P(O_2)$ is the probability that door 2 is removed from selection without the condition that the prize is in any particular location, and we must calculate that directly.

If we consider all of the possible locations where the prize may be (behind each of the three doors) along with their corresponding probabilities, we see that
$P(O_2)=P(D_1)P(O_2|D_1)+P(D_2)P(O_2|D_2)+P(D_3)P(O_2|D_3).$
We know that $P(D_1)=P(D_2)=P(D_3)=1/3$, and so we consider the other values in question. The first of these values, $P(O_2|D_1)$ is the probability that door 2 is removed given that the prize is behind door 1. We already know that this is 1/2. The second value, $P(O_2|D_2)$ is zero since we’d never have door 2 opened if the prize is behind door 2. The third value, $P(O_2|D_3)$ is 1/1 (or 100%) because the prize being behind door 3 and our selection of door 1 only leaves door 2 to be opened. Thus, we have computed:
$P(D_1|O_2)=\frac{P(O_2|D_1)P(D_1)}{P(O_2)}=\frac{\frac{1}{2}\cdot\frac{1}{3}}{\frac{1}{3}\left(\frac{1}{2}+0+1\right)}=\frac{1}{3}.$
If I just did your probability or stats homework, you’re welcome. (Please understand the result if this is the case.) What we have just shown is that keeping our selection of door at door 1 results in a 1/3 chance of us winning the prize. What about the probability of winning when we switch choices to door 3? We have
$P(D_3|O_2)=\frac{P(O_2|D_3)P(D_3)}{P(O_2)}.$
The denominator here is exactly the same (thank goodness for that). The first value in the numerator, $P(O_2|D_3)$ is 1 (or 100%) since it is the probability of opening door 2 given that the prize is behind door 3, in which case door 2 is the only possible door to open (remember that we have initially selected door 1). The second value, $P(D_3)$ is 1/3 since is represents the unconditional probability of the prize being behind door 3. Thus, the probability of winning if we change our choice to door 3 is
$P(D_3|O_2)=\frac{P(O_2|D_3)P(D_3)}{P(O_2)}=\frac{1\cdot\frac{1}{3}}{\frac{1}{2}}=\frac{2}{3}$
and the result has been proved. QED.

A Python Model

We’ll create a Python class that will model the Monty Hall problem. Each time the class is initialized, it will select a door at random (from the 3 available) representing the contestant’s selection. It will then have the host remove a door and present the contestant with the option to change doors. We’ll run the game one million times with both situations: one million times where the contestant switches doors and one million times where the contestant sticks with the original door choice. We’ll present winning percentages in both cases. The code is a bit verbose. (I’ve used this recently as an introductory example to object oriented programming in Python.) This is saved in a file named MontyHall.py.

#!/usr/bin/python2   # MontyHall.py - Jason B. Hill 2014   # This script models the Monty Hall problem: You are on a gameshow and are # presented with three doors. Behind one door (chosen at random) sits a prize. # You must select (guess) a door. If you guess correctly, you get the prize. # But, when you select your door, the show's host throws in a bit of a twist, # opening one of the other two doors ... specifically a door that doesn't # contain the prize. He asks you if you'd like to change your door selection, # to the only other remaining door (not the one you picked originally and not # the now open door). Do you change your selection? Why?   # Answer: Yes. ALWAYS change your selection. By sticking with your original # door selection, you have a 1/3 chance of getting the prize. By switching, you # end up with a 2/3 chance of getting the prize.   # I explain why on my code blog at: code.jasonbhill.com   # I'm going to write this as a tutorial example for object oriented programming # in Python. The idea is that we'll have a class object representing the game. # We'll have variables inside that class that hold information about the game # state, and methods (functions) that advance the game state.       import random # We need the random module to select integers in [1,3] randomly     def pick_door(): """ Function to pick a door. Returns an integer 1, 2, or 3 at random. """ return random.randint(1,3)     class MontyHall: """ Class to model the Monty Hall problem. """ def __init__(self): """ Creates an instance of the Monty Hall problem. This will be called automatically when the class is formed. """ # Pick a prize door randomly and store as a variable. self.prize_door = pick_door() # We'll create variables for the selected door and removed door as # well, but we won't set them just yet. self.selected_door = None self.removed_door = None   def select_door(self): """ Randomly selects a door for the contestant. """ self.selected_door = pick_door()   def remove_door(self): """ This is how the host removes a (non-prize/non-selected) door.. """ # Pick a door at random. d = pick_door() # If that door is the prize door or the contestant's door, re-pick. while d == self.selected_door or d == self.prize_door: d = pick_door() # set the removed door to d self.removed_door = d   def switch_choice(self): """ Switches the selected door once a non-prize door is removed. """ # 1+2+3=6. There's only one choice of door if we switch our selection. self.selected_door = 6 - self.selected_door - self.removed_door   def user_wins(self): """ Determine if the user wins. Return true on win, false on lose. """ if self.selected_door == self.prize_door: return True else: return False   def run_game(self, switch=True): """ Once a problem is initialized, run the game.   'switch' determines if the user changes selection during the game. """ # The user selects a door. self.select_door() # The host removes a door. self.remove_door() # The user can switch selection of doors. if switch: self.switch_choice() # Determine if the user wins. return self.user_wins()       # Now, we'll run the game. When asked if we want to switch door selection when # the door is removed by the host, we'll say 'yes' and switch our choice of # door. We'll run that experiment one million times, always switching choice # when given the chance. Here's what that looks like:   wins, losses = 0, 0 for i in range(1000000): # make an instance of the game, call it 'm' m = MontyHall() # run the game and switch choice of door. if m.run_game(switch=True): # a return value of True means we've won wins += 1 else: # a return value of False means we've lost losses += 1   # Now that the game has been run one million times, compute/display stats. perc_win = 100.0*wins / (wins+losses)   print "\nOne million Monty Hall games (with switching):" print " won:", wins, "games" print " lost:", losses, "games" print " odds: %.2f%% winning percentage" % perc_win     # Now, we'll run the game one million times and always stick to our original # door selection every single time.   wins, losses = 0, 0 for i in range(1000000): # make an instance of the game, call it 'm' m = MontyHall() # run the game and switch choice of door. if m.run_game(switch=False): # a return value of True means we've won wins += 1 else: # a return value of False means we've lost losses += 1   # Now that the game has been run one million times, compute/display stats. perc_win = 100.0*wins / (wins+losses)   print "\nOne million Monty Hall games (staying with original choice):" print " won:", wins, "games" print " lost:", losses, "games" print " odds: %.2f%% winning percentage\n" % perc_win

Running the program, we get the following output.

$python MontyHall.py One million Monty Hall games (with switching): won: 666675 games lost: 333325 games odds: 66.67% winning percentage One million Monty Hall games (staying with original choice): won: 333779 games lost: 666221 games odds: 33.38% winning percentage As of writing this (March 31, 2014), the version of Boost included in Debian Wheezy is 1.49. This is a bit problematic for me, as some very useful things are included in more recent versions (e.g., better syslog event handlers in the logging library). One of the problems I sometimes face with Boost is in remembering how to compile/install it. The Boost “Getting Started” pages basically say that most of Boost is header-only and doesn’t need to be compiled, but that there are 16 libraries that need to be compiled. Unfortunately, it stops there. So, this is for future reference and for others that may find it useful. How does one compile/install all of Boost? The first thing we need to do is make sure all the dependencies for Boost are available. Using your favorite package manager, make sure you have the following. (Depending on your distro, they may be named slightly differently. But, make sure you have the dev version when required.) - build-essential - g++ - python-dev - autotools-dev - libicu-dev - libbz2-dev Next, download and decompress the Boost version that you’d like. I’m using 1.55. $ wget http://downloads.sourceforge.net/project/boost/boost/1.55.0/boost_1_55_0.tar.gz $gunzip boost_1_55_0.tar.gz$ tar xvf boost_1_55_0.tar

Finally, as root, build and install. Since I’m placing the installation in /usr/local (as is the default), I’m not issuing any prefix declaration here. There are really just two commands:

(root) $./bootstrap.sh (root)$ ./b2 --with=all -j 2 install

The -j n option will build the libraries using n threads, so this can be changed based on your system. After completion, you should be able to include and link against the required headers and libraries.h

I use GNOME 3 in OpenSUSE relatively frequently at work. I’m familiar with this “shell” flavor of GNOME and have watched it evolve over recent years. But, I’ve stuck with XFCE on my personal machines for several reasons:

• XFCE is very customizable.
• It looks professional and doesn’t throw useless flashy buttons or movements in front of you.
• While there are some pretty horrible themes, it generally hasn’t been hit with the ugly stick. This was my main reason for moving away from GNOME 2, which was my first Linux desktop environment after starting with Ubuntu years ago.

I revisited GNOME 3 today in Debian for my personal laptop. I’m sad to say that I still can’t use it. While there are great aspects to it, there are still bits that infuriate me.

Icon Sizes

Believe it or not, the issue of icon size was the original thing that pushed me from GNOME 2 to XFCE. Here’s a sample of what toolbar icons look like in GNOME 2.

This is the icon equivalent of bad kerning. Icons are simply stuffed full-size into a toolbar, with no space between them and no vertical padding in the toolbar. This drives me nuts and it looks horrendous. At least in XFCE, I get reliable spacing and padding around icons to make them look reasonable. Yep… that’s the original reason why I switched to XFCE. GNOME2 was just hit with the ugly stick in a big way.

Moving on to GNOME 3, it gets worse. While the toolbars are basically gone, I now have the currently open application showing in the top toolbar as follows.     So, Iceweasel (the Debian equivalent of Firefox) is open here and the icon is so big that it won’t even fit in the toolbar. Plus, there’s text to show me the application name associated with the icon that is now unrecognizable because it’s so huge. This is simply GNOME 3 trying to add some design appeal, while forgetting that this Iceweasel icon and text do absolutely nothing, except for possibly closing Iceweasel. It gets worse. If I click on this monstrosity, it pops up a menu with a single option: “Quit,” which stays open until I click the icon/text again. It’s useless, beyond letting you know which program you’re using, in case you forget. And, you can’t get rid of this ugliness.

Sure, there are add-on extensions to remove this little useless graphic, but I haven’t been able to get them to work. So, I’m stuck with ugly, space-sucking, useless graphics on my desktop. Good job.

While we’re on the topic of icon sizes, I should mention that I was able to find an extension that would allow me to change the notification icons on the right-hand size of the top toolbar in GNOME 3. They now, after being modified from their original form, look as follows.

That actually looks good. I like that. It’s clean and functional.

Activities

I like the mouse gestures associated with the “Activities” button in the upper-left side of my screen, but I don’t understand why it is called “Activities.” It would be more appropriately called “Stuff” or “Dashboard” or even “Oberoth,” the name of my laptop.

If you name it “Activities,” then it should actually be a route for accessing items that are actually activities: “write an email,” “browse a website,” etc. Instead, “activities” are things like “removable devices” and “windows.” I don’t even see why “activities” is there in the first place.

If I go in to “activities” and then click on “applications,” I get an array of huge icons associated with programs. I do like that I can start typing the name of a program and have this list be filtered. But, the array of icons and text isn’t incredibly easy to use. I was able to reduce the icon size a bit, but now we find ourselves in the following situation.

Just to recap, the six application options in these six icons are as follows: “gedit”, “GNU Im…”, “GParted,” “Hugin …”, “Hugin …”, “Hugin …”, “LibreO…”, “LibreO…”, and “LibreO…”. Only two of these nine options remove any doubt as to what they access. This is just horribly designed.

Conclusion

I could go on, but won’t. It’s this sort of nonsense that I don’t miss when using XFCE. I think the look and feel overall is nice, but it’s infuriating during actual use. I can’t find what I’m looking for, and thing I’m not looking for take up space. The clock, for instance, takes up about 3x the space on the GNOME 3 toolbar as my XFCE clock does on my XFCE toolbar.

I’ll try again in another year or so. Until then, XFCE still has my vote. By the way, below is a screenshot of what my XFCE environment currently looks like.

Problem

I have been using XFCE as a desktop environment in Linux ever since I became frustrated with GNOME 2 (it was hit with the ugly stick). I’ve used a few others (GNOME 3, LXDE, etc.) and some window managers as well (OpenBox), but I keep coming back to XFCE for its simplicity and elegance. What annoyed me for some time though was that the display manager of XFCE didn’t appear to support multiple monitors. In reality, you just have to know how to do it … which isn’t terribly complicated.

XRANDR

Xrandr is the “primitive command line interface to the RandR extension.” This is what will allow us to enable and configure multiple monitors via the terminal. The first thing we’ll do is see what sort of displays are available. I’m running a Lenovo ThinkPad on a dock, so I’m going to have a load of display ports at my disposal. Xrandr is almost certainly installed if you’re running a distro such as Debian, Ubuntu/Xubuntu, etc.

We can get a list of available display ports by simply typing ‘xrandr’ at the prompt.

jason@oberoth:~$xrandr Screen 0: minimum 320 x 200, current 1600 x 900, maximum 8192 x 8192 LVDS1 connected 1600x900+0+0 (normal left inverted right x axis y axis) 309mm x 174mm 1600x900 60.0*+ 40.0 1440x900 59.9 1360x768 59.8 60.0 1152x864 60.0 1024x768 60.0 800x600 60.3 56.2 640x480 59.9 VGA1 disconnected (normal left inverted right x axis y axis) HDMI1 disconnected (normal left inverted right x axis y axis) DP1 disconnected (normal left inverted right x axis y axis) HDMI2 disconnected (normal left inverted right x axis y axis) HDMI3 disconnected (normal left inverted right x axis y axis) DP2 disconnected (normal left inverted right x axis y axis) DP3 disconnected (normal left inverted right x axis y axis) As xrandr shows, I’m currently running with only my laptop screen, which is named LVDS1. This screen has a resolution of 1600×900, and it’s currently running at that resolution. If I connect a VGA cable from my dock to my monitor, I’ll get a slightly different output. jason@oberoth:~$ xrandr Screen 0: minimum 320 x 200, current 1600 x 900, maximum 8192 x 8192 LVDS1 connected 1600x900+0+0 (normal left inverted right x axis y axis) 309mm x 174mm 1600x900 60.0*+ 40.0 1440x900 59.9 1360x768 59.8 60.0 1152x864 60.0 1024x768 60.0 800x600 60.3 56.2 640x480 59.9 VGA1 connected (normal left inverted right x axis y axis) 1920x1080 60.0 + 1600x1200 60.0 1680x1050 60.0 1280x1024 75.0 60.0 1440x900 75.0 59.9 1280x960 60.0 1280x800 59.8 1152x864 75.0 1024x768 75.1 70.1 60.0 832x624 74.6 800x600 72.2 75.0 60.3 56.2 640x480 75.0 72.8 66.7 60.0 720x400 70.1 HDMI1 disconnected (normal left inverted right x axis y axis) DP1 disconnected (normal left inverted right x axis y axis) HDMI2 disconnected (normal left inverted right x axis y axis) HDMI3 disconnected (normal left inverted right x axis y axis) DP2 disconnected (normal left inverted right x axis y axis) DP3 disconnected (normal left inverted right x axis y axis)

Now, VGA1 is connected and the data shows that it is capable of running at 1920×1080. There are multiple ways to tell xrandr that we’d like to use both LVDS1 (laptop monitor) and the external monitor connected at VGA1. First, we need to determine the orientation of how we’d like to set things up. My preference is to have my laptop monitor maintain my XFCE desktop toolbar and be viewed as the “primary” screen. My external monitor will sit off to the right size of the laptop screen, and I need to tell xrandr where to place this monitor in reference to my laptop screen. I do that like so.

xrandr --output LVDS1 --auto --primary --output VGA1 --auto --right-of LVDS1

This tells xrandr that the first output is LVDS1, which is set to automatically register resolution and also be set as the primary output. The second output is VGA1, which also has resolution set automatically and resides to the right of LVDS1. Once we issue that command, the second monitor becomes active and I can move my mouse pointer and drag windows between monitors. If we issue an xrandr command now, we see this.

jason@oberoth:~xrandr Screen 0: minimum 320 x 200, current 3520 x 1080, maximum 8192 x 8192 LVDS1 connected 1600x900+0+0 (normal left inverted right x axis y axis) 309mm x 174mm 1600x900 60.0*+ 40.0 1440x900 59.9 1360x768 59.8 60.0 1152x864 60.0 1024x768 60.0 800x600 60.3 56.2 640x480 59.9 VGA1 connected 1920x1080+1600+0 (normal left inverted right x axis y axis) 477mm x 268mm 1920x1080 60.0*+ 1600x1200 60.0 1680x1050 60.0 1280x1024 75.0 60.0 1440x900 75.0 59.9 1280x960 60.0 1280x800 59.8 1152x864 75.0 1024x768 75.1 70.1 60.0 832x624 74.6 800x600 72.2 75.0 60.3 56.2 640x480 75.0 72.8 66.7 60.0 720x400 70.1 HDMI1 disconnected (normal left inverted right x axis y axis) DP1 disconnected (normal left inverted right x axis y axis) HDMI2 disconnected (normal left inverted right x axis y axis) HDMI3 disconnected (normal left inverted right x axis y axis) DP2 disconnected (normal left inverted right x axis y axis) DP3 disconnected (normal left inverted right x axis y axis) We see that screen 0 now registers as being 3520×1080, and it holds two outputs: LVDS1 with resolution 1600×900 at 0×0 and VGA1 with resolution 1920×1080 at location 1600×0. The only awkward issue here is that my external monitor is taller than my laptop. Right now, the monitors are aligned vertically at the top. What if I want them aligned vertically by the bottom of the screens? Keeping in mind that the external monitor is 180 pixels taller than the laptop, I could position the external monitor manually as so. xrandr --output LVDS1 --auto --primary --output VGA1 --auto --pos 1600x-180 Rotating the External Monitor One thing I like to do while coding, or even reading webpages, is to rotate my external monitor for a taller reading profile. I do that as follows. xrandr --output LVDS1 --auto --primary --output VGA1 --auto --right-of LVDS1 --rotate right Maybe some time in the future I’ll come up with a script to automate this. There’s this annoying thing that’s happening at work. It started with one or two sites being blocked from access while on the company network. It has grown to the point where standard, run-of-the-mill stuff is blocked. For instance, the Google image search for “Linux Mint XFCE” (I’ve been considering switching my work desktop to Mint) turns up results that are blocked by the network. Here’s a screenshot of the warning that a Linux Mint screenshot image link produces. It would be one thing if I only ran into this “access denied” warning while browsing during lunch. But, this keeps popping up for work-related searches, including the recent issue of not being able to access a page about Python’s Twisted module because it … apparently … contained information about weapons. Forwarding HTTP Traffic over SSH So, what I’m going to do is to reroute all of my browser traffic over SSH and through a server that I administrate. Come to think of it, this is a safe solution for many situations (travels, working on sensitive data over an open wifi connection, etc.). First, you need to have SSH access to a server that can act as a proxy for your web browser connection. I’m using a server that I’ll call proxy.myserver.com. We’ll set up the forwarding here. There are a few command options that we use for the ssh command here. Namely, the D option specifies a local application level port forward and the N tells proxy.myserver.com that we don’t want to interact… we’ll just open a connection in the terminal and let it sit there without a prompt. The real data transfer work is going to be done through the web browser. Here’s how we set it up, forwarding port 9999 from my desktop to proxy.myserver.com, I enter the following in the terminal on my desktop.  ssh -ND 9999 jason@proxy.myserver.com

After entering the password, it just sits there, which is exactly what should happen. Now, let’s set up the browser. I’m using Iceweasel (Firefox) on Debian. By entering the options menu at Edit->Preferences and finding our way to Advanced->Connection->Settings, we get the following menu.

By entering a manual proxy configuration with a SOCKS host, we will be able to utilize our SSH port 9999 connection.

Now, Firefox/Iceweasel will use the proxy connection through proxy.myserver.com to fetch any web traffic. If for some reason that connection dies, we’ll get a warning from the browser stating that the proxy isn’t accepting any connections.

I was person #100 to solve Project Euler problem 451 a while back. The problem and solution are detailed in the source code below. The solution was found in 28 seconds on a 20-core Xeon machine. As noted in the code, using the Chinese Remainder Theorem could help efficiency, but I don’t have that coded currently and so I iterated over potential C.R.T. candidates instead, which provided a quickly coded and self-contained solution.

  /* * Jason B. Hill / jason@jasonbhill.com * pe451.c / solves Project Euler problem 451 * build: gcc pe451.c -fopenmp -O3 -o pe451 * execute: ./pe451 * * I logged in to see that problem 451 had around 90 solutions, and decided to * attempt to be one of the first 100 solvers to enter the proper solution. As * such, I coded this a bit quickly. A "clean" solution would make use of the * Chinese Remainder Theorem, and the timing could probably be brought down * a tiny bit. But, all things considered, this runs very fast itself. * * Problem: Consider the number 15. There are eight positive integers less than * 15 which are coprime to 15: 1, 2, 4, 7, 8, 11, 13, 14. The modular inverses * of these numbers modulo 15 are: 1, 8, 4, 13, 2, 11, 7, 14 because * 1*1 mod 15=1 * 2*8=16 mod 15=1 * 4*4=16 mod 15=1 * 7*13=91 mod 15=1 * 11*11=121 mod 15=1 * 14*14=196 mod 15=1 * Let I(n) be the largest positive number m smaller than n-1 such that the * modular inverse of m modulo n equals m itself. So I(15)=11. Also I(100)=51 * and I(7)=1. Find the sum of I(n) for 3<=n<=2·10**7. * * The solution is documented in the code below. The code is in C with OpenMP. * * The result is found in 28 seconds on a 20-core Xeon and and 2:48 on a 2-core * core i7. Here is some of the output: * prime sieve created - 0.110000 seconds * list of 1270607 primes created - 0.050000 seconds * list of prime factors created for integers <= 20000000 - 5.350000 seconds * prime factor exponents computed for all integers - 0.980000 seconds */   #include <omp.h> #include <stdio.h> #include <stdlib.h> #include <stdint.h> #include <time.h>     /*****************************************************************************/ /* The following structs and routines are for integer factorization */ /*****************************************************************************/   /* * Notes on factoring: * * We need to consider the prime factorization of positive integers less than * 2*10^7. How many distinct prime factors can such numbers have? That's easy. * Since the product of the 8 smallest primes (2*3*5*7*11*13*17*19 = 9,699,690) * is less than 2*10^7 and the 9 smallest primes (9,699,690*23 = 223,092,870) * is greater than 2*10^7, we know that each positive integer less than 2*10^7 * can only have 8 distinct prime factors. Let's create a C-struct that can * keep track of this data for each factored integer. */   /* * A struct to hold factor information for an integer. * p holds a list of distinct prime factors. (at most 8) * e holds a list of exponents for the prime factors in p. * num_factors says how long the lists p and e are. */ struct factors_t { uint64_t p[8]; // list of distinct primes uint8_t e[8]; // list of exponents for the primes p[] uint64_t pp[8]; // prime power factor (p[i]**e[i]) uint8_t num_factors; // how many distinct prime factors };     /* * More notes on factoring: * * Every positive integer n has a prime factor <= sqrt(n). We use this fact to * build a prime sieve. First, we construct a function to return the square * root of a uint64_t. Then, we'll use that function to create a sieve (which * is returned as a pointer/list of uint8_ts... effectively Booleans). */   /* * Return the square root of a uint64_t */ uint64_t sqrt_uint64(uint64_t n) { uint8_t shift = 1; uint64_t res, s;   while ((1<<shift) < n) shift += 1; res = (1<<((shift>>1) + 1)); while (1) { s = (n/res + res)/2; if (s >= res) return res; res = s; } }     /* * Return a prime sieve identifying all primes <= limit. * This is just a list of uint8_t's where 0 means the index is composite and * 1 means the index is prime. */ uint8_t * prime_sieve(uint64_t limit) { uint8_t *sieve; // this will be the pointer that is returned uint64_t i,j; uint64_t s = sqrt_uint64(limit);   // allocate memory for sieve sieve = malloc((limit + 1) * sizeof(uint8_t));   // set initial values in the sieve sieve[0] = 0; sieve[1] = 0; sieve[2] = 1;   // set other initial odd values in sieve to 1, even values to 0. for (i=3;i<=limit;i++) { if (i%2==0) sieve[i] = 0; else sieve[i] = 1; }   // unset composite numbers (evens are already unset) for (i=3;i<=(s+1);i=i+2) { if (sieve[i]==1) { j = i*3; // sieve[i] prime, sieve[2*i] already 0 while (j<=limit) { sieve[j] = 0; j += 2*i; } } } return sieve; }     /* * Determine if a value is prime using a provided sieve. */ _Bool is_prime(uint64_t n, uint8_t * sieve) { if ((1 & sieve[n]) == 1) return 1; else return 0; }     /*****************************************************************************/ /* The following functions are specific to Project Euler problem 451 */ /*****************************************************************************/   /* * Notes: * * Given a number n, we're looking for the largest m < n-1 such that m^2 is 1 * modulo n. That's a mouthful. First, field theory tells us that we're only * interested in integers that are relatively prime with n (which is why the * problem asks for m < n-1 ... otherwise it would be a lot easier). * * Now, if n=p^k is a prime power and we have x^2 = 1 mod p^k, we consider: * a) 1 is a solution for x since 1*1=1 mod p^k = 1. * b) p^k-1 is a solution since (p^k-1)^2=p^(2k)-2p^k+1 (mod p^k) = 1. But, * that's too big since we need m < n=p^k-1. :-( * c) Any other solutions? Well... that's sort of complicated. Ugh. Let's first * consider the case when p=2. For n=p^1=2^1, we only have x=1 as a * solution. When n=p^2=2^2, we only have x=1 and x=3. For n=p^k=2^k, we * also find that 2^(k-1)+1 and 2^(k-1)-1 square to 1 modulo 2^k. For other * primes, this doesn't happen and we only have 1 and p^k-1 as solutions. * * Thus, what we're going to do is as follows: * 1) Factor each integer within the given range (factoring a range of numbers * is much faster than factoring each one individually). We'll store that * factorization information using the struct defined above. * 2) Because of step 1, determining if a number is prime or a power of a prime * is easy. In case of primes and prime powers, we consider if the prime is * even or odd, returning the appropriate maximum square root of 1 directly. * 3) When the integer is composite having more than a single distinct prime * factor, we use the Chinese Remainder Theorem "in spirit" and iterate * over potential candidates (instead of computing the C.R.T. result * directly). We're only considering possible m that are both relatively * prime with n AND such that m^2 is +1 or -1 modulo the prime power * factors of n. When a candidate does not satisfy these properties, we * simply move to the next candidate. */   /* * Compute the next candidate (the next largest number that may satisfy the * given equivalence relations). This is done relative to the largest non-2 * prime power in the factorization of n. * * n is the integer for which we're computing the largest square root. * m is the largest odd prime power factor of n. * c is the current candidate (assume uninitialized when set to 0). */ uint64_t next_candidate(uint64_t n, uint64_t m, uint64_t c) { uint64_t r;   // check if c is initialized if (c==0) c = n-1; // we can only consider values of \pm1 mod m r = c % m; if (r==m-1) return c-(m-2); else return c-2; }     /* * Determine if the current candidate is (1) relatively prime to n and, if it * is, (2) a square root of unity modulo n. We can tweak this later for timing * as the test for being relatively prime may cut performance a bit. We'll see. * * n is the integer for which we're computing the largest square root. * cnd is the current candidate. * factors is a list of factors_t structs (prime, exponent info) */ _Bool verify_candidate(uint64_t n, uint64_t cnd, struct factors_t * factors) { uint8_t i,j; uint64_t pp;   // verify that cnd is not divisible by any prime in factors[n].p[] for (i=0;i<factors[n].num_factors;i++) { if (cnd % factors[n].p[i] == 0) return 0; } // verify cnd modulo 2 when exponent of 2 in n is > 2 if (factors[n].p[0] == 2 && factors[n].e[0] > 2) { pp = factors[n].pp[0]>>1; if (cnd % pp != 1 && cnd % pp != pp-1) return 0; } // verify other primes for (i=0;i<factors[n].num_factors;i++) { if (factors[n].p[i] == 2) continue; pp = factors[n].pp[i]; if (cnd % pp != 1 && cnd % pp != pp-1) return 0; } return 1; }     /* * Given a positive integer n, find the largest positive m < n-1 such that * gcd(n,m)=1 and m**2 (mod n) = 1. */ uint64_t largest_sqrt(uint64_t n, struct factors_t * factors) { uint8_t i; uint32_t j;   // if n is an odd prime, or a power of an odd prime, return 1 if (factors[n].num_factors == 1 && factors[n].p[0] != 2) return 1;   // if n is a power of 2 if (factors[n].num_factors == 1) { // if n is 2 or 4 if (factors[n].e[0] < 3) return 1; // if n is 2**e for e >= 3 return (factors[n].pp[0]>>1)+1; }   // find the maximum odd prime power factor of n uint64_t pp = 1; uint64_t cnd;   for (i=0;i<factors[n].num_factors;i++) { if (factors[n].p[i] != 2 && factors[n].pp[i] > pp) { pp = factors[n].pp[i]; } }   // get the first candidate w.r.t. moppf cnd = next_candidate(n, pp, 0);   while (!verify_candidate(n, cnd, factors)) cnd = next_candidate(n, pp, cnd);   return cnd; }     /*****************************************************************************/ /* Execute */ /*****************************************************************************/ int main() { uint64_t s = 0; // final output value uint64_t i,j,k; // iterators uint64_t limit = 20000000; // upper bound for computations uint64_t num_primes = 0; // count for primes <= limit uint8_t e; // exponents for prime factoring int tid; // thread id for openmp double time_count; // timer clock_t start, start_w; // time variables   uint8_t *sieve; // prime sieve uint64_t *primes; // list of primes struct factors_t *factors; // prime factors and exponents     /* start outer timer */ start_w = clock();     /* make the prime sieve */ start = clock(); sieve = prime_sieve(limit); time_count = (double)(clock() - start) / CLOCKS_PER_SEC; printf("prime sieve created - %f seconds\n", time_count);     /* form list of primes from sieve */ start = clock(); // compute the number of primes in sieve for (i=2;i<=limit;i++) { if (is_prime(i, sieve)) { num_primes = num_primes + 1; } } primes = malloc(num_primes * sizeof(uint64_t)); j=0; for (i=2;i<=limit;i++) { if (is_prime(i, sieve)) { primes[j] = i; j++; } } time_count = (double)(clock() - start) / CLOCKS_PER_SEC; printf("list of %llu primes created - %f seconds\n", num_primes, time_count);     /* fill out a factors_t struct for each integer below limit */ start = clock(); // allocate memory for factors_t factors = malloc(sizeof(struct factors_t) * (limit + 1)); // set the initial number of factors for each number to 0 for (i=1;i<=limit;i++) factors[i].num_factors=0; // for each prime, add that prime as a factor to each of its multiples for (i=0;i<num_primes;i++) { j = 1; // start at 1*p for each prime p while (j*primes[i]<=limit) { k = factors[j*primes[i]].num_factors; // get proper index for p factors[j*primes[i]].p[k] = primes[i]; // add prime to p factors[j*primes[i]].num_factors++; // increase index j++; // increase multiple of prime } } time_count = (double)(clock() - start) / CLOCKS_PER_SEC; printf("list of prime factors created for integers <= %llu - %f seconds\n", limit, time_count);     /* compute exponents for each prime in factor lists */ start = clock(); for (i=2;i<=limit;i++) { for (j=0;j<factors[i].num_factors;j++) { e=1; k=factors[i].p[j]; while (i % (k*factors[i].p[j]) == 0) { e++; k*=factors[i].p[j]; } factors[i].e[j] = e; factors[i].pp[j] = k; } } time_count = (double)(clock() - start) / CLOCKS_PER_SEC; printf("prime factor exponents computed for all integers - %f seconds\n", time_count);     /* find largest square root of unity for each integer under limit; add to s */ #pragma omp parallel for reduction(+:s) shared(factors) schedule(dynamic,1000) for (i=3;i<=limit;i++) { s = s+largest_sqrt(i, factors); } time_count = (double)(clock() - start) / CLOCKS_PER_SEC; printf("result for 3<=i<=%llu: %llu - %f seconds\n", limit, s, time_count);     time_count = (double)(clock() - start_w) / CLOCKS_PER_SEC; printf("\ntotal time: %f seconds\n\n", time_count);   free(sieve); free(primes); free(factors);   return 0; }

Problem

The fraction $\frac{49}{98}$ is a curious fraction, as an inexperienced mathematician in attempting to simplify it may incorrectly believe that $\frac{49}{98}=\frac{4}{8}$, which is correct, is obtained by cancelling the 9s.

We shall consider fractions like $\frac{30}{50}=\frac{3}{5}$ to be trivial examples.

There are exactly four non-trivial examples of this type of fraction, less than one in value, and containing two digits in the numerator and denominator.

If the product of these four fractions is given in its lowest common terms, find the value of the denominator.

Solution in Python

We really only need to iterate from 10 to 99 for the numberator and from the numerator to 99 for the denominator, since we require that the fraction be less than one. Also, we can ignore any denominators that have zeros (since the only other possible integer to match in the numerator is nonzero). We form a list from the numerator and denominator, check for redundancy, remove the redundancies when found, then cross multiply to find if the fractions are equivalent to their redundant-removed siblings.

import time import fractions   start = time.time()   p = fractions.Fraction(1,1)   for a in range(10, 100, 1): for b in range(a+1, 100, 1): if b % 10 == 0 or a == b: continue La, Lb = [a/10, a%10], [b/10, b%10] if any(i in Lb for i in La) and not all(i in Lb for i in La): if La[0] in Lb: x = La[0] else: x = La[1] La.remove(x) Lb.remove(x) if a*Lb[0] == b*La[0]: p *= fractions.Fraction(La[0],Lb[0])   elapsed = time.time() - start   print "result %s found in %s seconds" % (p, elapsed)

We’ve used Python’s fractions module to make things a bit easier. When executed, we get the following.

result 1/100 found in 0.00531482696533 seconds

Problem

Consider all integer combinations of $ab$ for $2\le a\le 5$ and $2\le b\le 5$:

22=4, 23=8, 24=16, 25=32
32=9, 33=27, 34=81, 35=243
42=16, 43=64, 44=256, 45=1024
52=25, 53=125, 54=625, 55=3125

If they are then placed in numerical order, with any repeats removed, we get the following sequence of 15 distinct terms:

4, 8, 9, 16, 25, 27, 32, 64, 81, 125, 243, 256, 625, 1024, 3125

How many distinct terms are in the sequence generated by $ab$ for $2\le a\le 100$ and $2\le b\le 100$?

Solution in Python

Update: After a conversation with my friend and co-worker Larry, I’ve changed just two lines of this code – using a Python set instead of a list. That has taken the runtime down considerably. The set version runs in about 2% the time of the list version. This is because checking if an element is in a set is a constant time operation, while checking membership in a list is a linear time operation.

This is one of those problems that can be over-engineered. If you simply code the basic happy-path Python version, you find that it runs very quickly and returns the result. You could consider unique factorizations and so on … and that was my initial idea, but coding the “brute-force” Python version is all you really need here.

#!/usr/bin/python   import time   # L = [] # old version with a list L = set() # new version with a set   limit = 101   start = time.time()   for a in range(2,limit): for b in range(2,limit): c = a**b if c in L: pass # else: L.append(c) # old list version L.add(c) # new version with set   elapsed = time.time() - start   print "found %s ints in %s seconds" % (len(L), elapsed)

When executed, we get the correct result in under a second.

$python 29.py # found 9183 ints in 0.66696190834 seconds # old list version found 9183 ints in 0.0150249004364 seconds Problem The$n^\text{th}$term of the sequence of triangle numbers is given by,$t_n=\frac{1}{2}n(n+1)$; so the first ten triangle numbers are: 1, 3, 6, 10, 15, 21, 28, 36, 45, 55, … By converting each letter in a word to a number corresponding to its alphabetical position and adding these values we form a word value. For example, the word value for SKY is 19 + 11 + 25 = 55 =$t_{10}$. If the word value is a triangle number then we shall call the word a triangle word. Using words.txt (right click and ‘Save Link/Target As…’), a 16K text file containing nearly two-thousand common English words, how many are triangle words? Solution in Python #!/usr/bin/python import time # start the clock start = time.time() # turn the string of words into a list of python strings with open('words.txt','r') as f: words = f.read().split(',') words = [list(word.strip('\"')) for word in words] f.close() # we should have an idea of how long the longest word is, # giving us an idea of the magnitude of the triangle words m = max([len(word) for word in words]) # form triangle numbers up to the given range triangles = [n*(n + 1)/2 for n in range(1,2*m)] # make a dictionary map for character values vals = {} s = list("ABCDEFGHIJKLMNOPQRSTUVWXYZ") for c in s: vals[c] = s.index(c) + 1 # count triangle words triangle_words = 0 for word in words: if sum([vals[c] for c in word]) in triangles: triangle_words += 1 # terminate the clock elapsed = time.time() - start print "found %s triangle words in %s seconds" % (triangle_words, elapsed) When executed, this returns as follows. found 162 triangle words in 0.00315809249878 seconds Problem An irrational decimal fraction is created by concatenating the positive integers: $0.123456789101112131415161718192021...$ It can be seen that the 12th digit of the fractional part is 1. If$d_n\$ represents the nth digit of the fractional part, find the value of the following expression.

$d_1\times d_{10}\times d_{100}\times d_{1000}\times d_{10000}\times d_{100000}\times d_{1000000}$

Python Solution

This was coded in Python using the Sage Notebook (for speed). The idea is to concatenate the string representations of increasing integers into a single string, turn that string into a list with the appropriate index range, and take the appropriate product from within that list. I use the range 400000 to form the string, guessing (correctly) that this would provide enough digits in the resulting string/list.

import time   start = time.time()   s = "" for i in range(400000): s += str(i) d = list(s) n = prod([ZZ(d[10**i]) for i in range(7)])   elapsed = time.time() - start   print "result %s found in %s seconds" % (n, elapsed)

When executed, we get the following result.

result 210 found in 0.369875907898 seconds