I got interested in this problem several years ago, as one of several unchosen options for a project in a graduate course that my wife was taking (or maybe it was someone else from work, I don’t remember).  The problem was to implement an efficient algorithm for solving the 2x2x2 Rubik’s cube.  I like this problem because it has several nice self-contained nuggets of mathematics in it.  The motivation for this post is one of those little sub-problem nuggets (dense collision-free hashing of permutations) that makes for a particularly simple and elegant algorithm, solving not just a given scrambled cube but all possible configurations, in just 100 lines of Python code.

The 2x2x2 cube is much simpler than the more common 3x3x3 cube.  The state space is obviously smaller; there are only or 3,674,160 possible configurations, compared with over for the larger cube.  Also, there is a simple and compact representation for those configurations and moves between them.  To see this, imagine fixing one corner cubie by holding it with the fingers of one hand, and applying moves by rotating any of the three faces not involving the fixed cubie (see image below; the fixed cubie is the one not visible).  With these 6 basic moves– clockwise and counterclockwise quarter turns of 3 faces– it is possible to realize all possible cube states, consisting of permutations and orientations of the 7 non-fixed cubies.

Scrambled 2x2x2 Rubik’s cube.

From this perspective, we represent a cube state as a pair , where is a permutation in  and is a 7-vector with elements in the finite field .  The permutation indicates the positions of the cubies, and the elements of the vector indicate the 3 possible orientations of each of the cubies.  (Although we can realize every permutation, not all possible combinations of orientations are possible; the orientations of any 6 cubies determine the remaining orientation.  We can impose this constraint by requiring the sum modulo 3 of elements of to be zero.)

Each of the six basic moves may also be represented in the same way, with the calculation of a new cube state resulting from applying move to cube state given by

where the group action is the natural one permuting the coordinates of .

Ok, so where are we?  The basic idea is to compute the graph of all possible cube states, where the 6 directed edges leaving each vertex correspond to the 6 possible quarter turn moves.  If we then compute– and store– the shortest path spanning tree rooted at the vertex corresponding to the initial solved state, we effectively have a single “map” of the solution for every possible scrambled configuration of the cube.

Now for the interesting part, and the reason for this post: a paper by Myrvold and Ruskey (see the reference below) describes a very handy, linear-time algorithm for converting permutations in to and from distinct integers in .  This is handy in this case because it allows us to use simple integer array indices in the adjacency list representation of the cube graph, as well as in the predecessor vector representation of the shortest path spanning tree.  For reference, here is my Python implementation of their algorithm, excerpted from the Rubik’s cube code:

def rankperm(p):
    """Return rank in [0, n!) of given permutation of range(n).

    W. Myrvold and F. Ruskey, Ranking and Unranking Permutations in
    Linear Time. Information Processing Letters, 79 (2001):281-284.
    """
    p = np.array(p)
    q = np.array(p).argsort()
    r = 0
    for k in range(len(p) - 1, 0, -1):
        s = p[k]
        p[k], p[q[k]] = p[q[k]], p[k]
        q[k], q[s] = q[s], q[k]
        r += s * np.math.factorial(k)
    return r

def unrankperm(r, n):
    """Return permutation of range(n) with given rank.

    W. Myrvold and F. Ruskey, Ranking and Unranking Permutations in
    Linear Time. Information Processing Letters, 79 (2001):281-284.
    """
    p = list(range(n))
    for k in range(n - 1, 0, -1):
        s, r = divmod(r, np.math.factorial(k))
        p[k], p[s] = p[s], p[k]
    return p

You can download the full implementation at the usual location here, as well as on GitHub.  In addition to the “business logic” that computes and saves the cube graph, there is also a GUI that lets you play with the cube or automatically solve it.  Both pieces of code use Visual Python, the former using Numpy which is bundled with VPython, and the latter using the 3D graphics.  The image above is a screenshot; as you can see, this is the result of just a few modifications and additions to the 3x3x3 version from earlier this year.

Reference:

1. W. Myrvold and F. Ruskey, Ranking and Unranking Permutations in Linear Time.  Information Processing Letters, 79 (2001):281-284. [PDF]

Although the jacket blurb for this blog admits the possibility of the occasional rant, I do my best to confine my thoughts here to those that highlight some specific idea, either mathematical, technological, or otherwise scientific.  My usual yardstick for determining whether a topic is worthwhile is that I should have something to contribute, be it a truly new result, or perhaps just an interesting different perspective on existing ideas.  Cane-shaking is for politics and religion.

But this week, after finishing a couple of great books, I plan to do nothing more than gush a bit about reading interesting words that are not my own.  I suppose I could hide behind the guise of a book review as a flimsy attempt at justifying not having much more to say than, “What are you doing still reading here?  You should be reading over there.”

The first book is A Stubbornly Persistent Illusion, edited by Stephen Hawking.  It is a collection of “The Essential Scientific Works” of Albert Einstein, containing not just his technical papers, but some fascinating public addresses and writings for a more general audience.  I found particularly interesting the rather intimate “self-obituary” in which he describes his progression from youth, through the past 67 years of his life, to his then-current perspective on science as a mode of thought:

“As the first way out there was religion, which is implanted into every child by way of the traditional education-machine.  Thus I came– despite the fact that I was the son of entirely irreligious (Jewish) parents– to a deep religiosity, which, however, found an abrupt ending at the age of 12.  Through the reading of popular scientific books I soon reached the conviction that much in the stories of the Bible could not be true…  Out yonder there was this huge world, which exists independently of us human beings and which stands before us like a great, eternal riddle, at least partially accessible to our inspection and thinking.”

I am always fascinated, and inspired, by these “human” views into the life of a man known to many people only by his equations.  Even the people who fundamentally change our understanding of the world– and Einstein was certainly one of them– are not cold, unfeeling machines or superhumans churning out equations.  They are simply men and women, whose very human and common desire to understand the world, and awe at the prospect of doing so, is just as inspiring to me as their very unique ability to actually realize that understanding.

The second book is the motivation for the title of this post, How to Solve It, by George Polya.  This is actually my second time reading this book, but I suppose I feel more influence from it now that I spend more time trying to teach.

This book is striking to me for two reasons.  First, it is inspiring in the same way that the previous book is: Polya knows the “Aha!” moment, knows how it is shared by both the student and the teacher.  But he suggests that those moments need not always be produced (i.e., proved) by the teacher to be consumed (i.e., verified) by the student.  It is possible for the student not just to learn solutions presented by others, but to learn how to solve it himself:

“The author remembers the time when he was a student himself, a somewhat ambitious student, eager to understand a little mathematics and physics.  He listened to lectures, read books, tried to take in the solutions and facts presented, but there was a question that disturbed him again and again: ‘Yes, the solution seems to work, it appears to be correct; but how is it possible to invent such a solution? … And how could I invent or discover such things by myself?'”

This phenomenon is what makes teaching such an exciting thing to do.  It is certainly rewarding when a student says, “I understand what you did there.”  But it is infinitely more rewarding when a student says, “I wondered what would happen if I did this, and look what I found…”

But the second striking feature of this book is the extent to which it succeeds at its goal.  Polya accomplishes something that at first glance seems oxymoronic: he systematizes heuristic.  More simply, he describes, with precision and clarifying examples, the common approaches to solving mathematical problems that may otherwise seem to have nothing in common.  In hindsight, much of the guidance seems like common sense.  But there seems to be something uncannily useful in making that common sense explicit.  As I read it, I find myself thinking that this is a book that everyone should read, not just teachers and students of mathematics.

Finally, it’s simply a fun book to read.  Polya knows how to write with some wit, such as when describing an example of how not to do it:

“The traditional mathematics professor of popular legend is absentminded.  He usually appears in public with a lost umbrella in each hand.  He prefers to face the blackboard and to turn his back on the class.  He writes a, he says b, he means c, but it should be d.”

The Boston Globe reported last week on an elderly couple who were making– and spending– huge amounts of money on the Massachusetts Cash WinFall state lottery.  In my usual habit, I started poking around in the mathematics involved, and also as usual, found more interesting effects than I expected.

The game is pretty simple: for each $2 ticket, you select 6 numbers from 1 to 46 (without replacement).  The parimutuel semi-weekly jackpot for matching all 6 numbers varies from $500,000 to about $2.5 million, with additional fixed payouts of $2, $5, $150, and $4,000 for matching 2 to 5 numbers, respectively.  As is usually the case, this game has negative expected return, with state lottery management taking about $1.10 from each $2 ticket sold.

What makes this game interesting is that when the jackpot reaches about $2 million without a winner– which happens several times a year– the excess jackpot is “rolled down” into larger fixed payouts for matching fewer than all 6 numbers.  This last happened on the 14 July drawing, with payouts increased to $26, $802, $19,507, and $2,392,699 for matching 3 through 6 numbers, respectively.  The result is a game with positive expected return for the ticket buyer… even if you never expect to win the jackpot!

But as Steve Jacobs used to say (or still says, for all I know), “Expectation isn’t everything.”  You still have a lousy 2% probability of actually making any money at all with just one ticket.  Expected value is only realized in the long run, considering the average result of a large number of repeated trials.

So the couple in the Globe article simply realized the long run, and actually executed that large number of repeated trials, buying about $614,000 worth of tickets.  This strategy not only has a nice expected return of about $184,000, but it is extremely low risk; the probability of making money jumps from 2% with one ticket to over 98% with 307,000 tickets.  The couple have already won nearly $1 million this year alone.

What got my attention and motivated this post was the following paragraph from the article:

Mark Kon, a professor of math and statistics at Boston University, calculated that a bettor buying even $10,000 worth of tickets would run a significant risk of losing more than they won during the July rolldown week. But someone who invested $100,000 in Cash WinFall tickets had a 72 percent chance of winning. Bettors like the Selbees, who spent at least $500,000 on the game, had almost no risk of losing money, Kon said.

The interesting question is, how to calculate that 72% probability figure?  As far as I can tell, this problem does not have a nice analytical solution.  The “big hammer” approach is to estimate by simply simulating the drawing many times.  But I wondered if there might be a more efficient way, either analytically or via a simpler form of estimation.

My first idea was to replace the actual game with a more tractable one: consider a lottery where each $2 ticket yields only two possible outcomes: either it wins a payout with probability , or it loses with probability .  We can compute the necessary values of and so that the single ticket outcome has the same (positive) expected value and variance as the actual game… but let us also conservatively assume that we never win the jackpot, so that the distribution is not so skewed.

In this simplified game, each ticket wins $4,510.95 with probability about 0.00052.  The expected return on a single ticket is the same as the actual (jackpot-less) game, about $.34, and the variance is the same as well.  The actual distribution is different, of course, but we should see the same general behavior as we buy more and more tickets: the overall expected return will increase linearly, and more importantly, the probability of making money will also increase, from the lowly 0.00052 with one ticket, approaching 1 for a sufficiently large number of tickets.

Right?

Well, not exactly.  The problem turned out to be more interesting than that.  Yes, we can efficiently compute the probability of winning money from tickets in our simplified game (I will leave it to the interested reader to work out the details).  But I was surprised to find that that probability is not monotonic as a function of , as the following plot shows.

Probability of winning vs. number of tickets.

The red curve corresponds to the actual lottery, and was generated via simulation; see the source code at the end of this post.  The behavior is what I expected: the more tickets you buy, the higher the probability that you win money.

The blue curve, corresponding to the simplified game, is more interesting.  First, just fixing the first and second moments of the distribution did indeed approximate the real behavior as well as I had hoped.  We can efficiently estimate the probabilities in the article: even buying $10,000 worth of tickets still only wins money about half the time; buying $100,000 wins about 75% of the time; and buying $614,000 wins money about 97% of the time.

But the unexpected and interesting behavior is the jaggedness of the blue curve.  There are many large jumps, at times nearly 20%, where the probability of an overall win decreases with the purchase of a single additional ticket.

ObPuzzle: what is going on here?

For reference, following is the source code for simulating the actual game, with commented parameters for the simplified game as well.

#include "math_Random.h"
#include <iostream>

int main()
{
    const int num_samples = 10000;

    // The simplified single-payout game.
    //const int num_outcomes = 2;
    //const double probability[] = {0.9994806467243879, 1.};
    //const double payoff[] = {-2, 4508.953581737101};

    // The actual game during the 14 July rolldown week.
    const int num_outcomes = 6;
    const double probability[] = {0.8312777261949867, 0.9776294385532591,
        0.9987251808751723, 0.999974270881075, 0.9999998932401705, 1.};
    const double payoff[] = {-2, 0, 24, 800, 19505, 2392697};

    math::Random rng;

    // Evaluate buying increasing numbers of tickets.
    for (int num_tickets = 25000; num_tickets <= 300000; num_tickets += 25000)
    {
        int count = 0;
        for (int i = 0; i < num_samples; ++i)
        {
            // Buy and cash in tickets.
            double win = 0;
            for (int j = 0; j < num_tickets; ++j)
            {
                double p = rng.nextDouble();
                for (int k = 0; k < num_outcomes; ++k)
                {
                    if (p <= probability[k])
                    {
                        win += payoff[k];
                        break;
                    }
                }
            }

            // Record whether we make or lose money.
            if (win >= 0)
            {
                ++count;
            }
        }

        // Display probability of winning (i.e., not losing money).
        std::cout << num_tickets << "\t" <<
            static_cast<double>(count) / num_samples << std::endl;
    }
}