Friday, March 30, 2012

Mega Millions

With the Mega Millions jackpot over $640 million, you might be asking yourself:

How do I use Factor to win the lottery???

You can find a good article on Forbes about calculating your odds of winning. We can also calculate it with Factor. The lottery works like this: you have to choose 5 "main" numbers (1 through 56) and then pick a "mega" number (1 through 46).

We can use the "n choose k" formula in math.combinatorics to compute the number of ways of picking "5 unordered outcomes from 56 numbers and then 1 of 46 possible mega numbers":

IN: scratchpad 56 5 nCk 46 * .
175711536

Sure enough, thats 175,711,536 possibilities. So if you buy one random lottery ticket, you have a 1 in 175,711,536 chance of winning, or a 0.00000000569114597006% chance. Not that great! What are the other odds?

Since the jackpot is so large, it's got to be worth playing, right? In fact, if you take the total jackpot and the odds of each winning category, you can find the expected value of a ticket:

IN: scratchpad {
                   640000000/175711536
                   250000/3904701
                   10000/689065
                   150/15313
                   150/13781
                   7/306
                   10/844
                   3/141
                   2/75
               } sum "$%.2f" printf
$3.82

Wow, $3.82 of expected value (not counting sharing the jackpot with someone who picked the same numbers, or the fact you can only win the highest prize you qualify for)! But, how do we pick our numbers... well, Factor to the rescue! Let's randomly sample our 5 main numbers and then pick a random mega number and output the result:

IN: scratchpad : mega ( -- seq )
                   56 [1,b] 5 sample natural-sort
                   46 [1,b] random suffix ;

In the spirit of fiver, we can generate these numbers to play:

IN: scratchpad 5 [ mega ] replicate .
{
    { 2 8 25 30 43 29 }
    { 6 15 31 41 47 33 }
    { 11 20 26 29 33 15 }
    { 3 6 19 23 32 15 }
    { 16 25 31 44 50 24 }
}

And tonight we can check the winning numbers and see if it worked!

Friday, March 2, 2012

Next Permutation

Yesterday, I noticed a programming challenge to find the "next greater permutation of digits" for a given number:

Given a number, find the next higher number that uses the same set of digits. For instance, given the number 38276, the next higher number that uses the digits 2 3 6 7 8 is 38627.

While reading the comments, I noticed that some of the C++ solutions used the std::next_permutation function that returns the "lexicographically next greater permutation of elements". Noticing that the math.combinatorics vocabulary lacks a next-permutation word, I thought it would be nice to contribute one to the Factor standard library.

std::next_permutation()

It's a useful place to start to get an overview and example of how the algorithm works. The C++ version is pretty dense and looks like this:

template<typename Iter>
bool next_permutation(Iter first, Iter last)
{
    if (first == last)
        return false;
    Iter i = first;
    ++i;
    if (i == last)
        return false;
    i = last;
    --i;

    for(;;)
    {
        Iter ii = i;
        --i;
        if (*i < *ii)
        {
            Iter j = last;
            while (!(*i < *--j))
            {}
            std::iter_swap(i, j);
            std::reverse(ii, last);
            return true;
        }
        if (i == first)
        {
            std::reverse(first, last);
            return false;
        }
    }
}

Example!

Walking through an example of the algorithm is particularly helpful to understand it:

As an example consider finding the next permutation of:

8342666411

The longest monotonically decreasing tail is 666411, and the corresponding head is 8342.

8342 666411

666411 is, by definition, reverse-ordered, and cannot be increased by permuting its elements. To find the next permutation, we must increase the head; a matter of finding the smallest tail element larger than the head’s final 2.

8342 666411

Walking back from the end of tail, the first element greater than 2 is 4.

8342 666411

Swap the 2 and the 4

8344 666211

Since head has increased, we now have a greater permutation. To reduce to the next permutation, we reverse tail, putting it into increasing order.

8344 112666

Join the head and tail back together. The permutation one greater than 8342666411 is 8344112666.

Implementation

We can use the steps in the example above to help organize our code. First, we find the "cut point" which is the index to the left of the longest monotonic tail:

: cut-point ( seq -- n )
    [ last ] keep [ [ > ] keep swap ] find-last drop nip ;

Next, we want to find the smallest element larger than the element at the "cut point" (searching from the end of the sequence):

: greater-from-last ( n seq -- i )
    [ nip ] [ nth ] 2bi [ > ] curry find-last drop ;

We then need a way to reverse the tail of a sequence (the ! denotes that this modifies the sequence in-place):

: reverse-tail! ( n seq -- seq )
    [ swap 1 + tail-slice reverse! drop ] keep ;

Putting this together gives us a word that looks a bit like our example. I have decided that in the case where the sequence reaches its lexicographic greatest order, we reverse it to its smallest ordering. This allows it to cycle through all possible permutations no matter where you start.

: (next-permutation) ( seq -- seq )
    dup cut-point [
        swap [ greater-from-last ] 2keep
        [ exchange ] [ reverse-tail! nip ] 3bi
    ] [ reverse! ] if* ;

We wrap this with a simple check to make sure the sequence is not empty. Arguably, we could instead check if the length is greater than 1 since a single element sequence has only possible permutation.

: next-permutation ( seq -- seq )
    dup [ ] [ drop (next-permutation) ] if-empty ;

Testing

Factor makes it easy to do unit tests. Here are some of the ones I've used to test this code:

[ "" ] [ "" next-permutation ] unit-test

[ "1" ] [ "1" next-permutation ] unit-test

[ "21" ] [ "12" next-permutation ] unit-test

[ "8344112666" ] [ "8342666411" next-permutation ] unit-test

[ "ABC" "ACB" "BAC" "BCA" "CAB" "CBA" "ABC" ]
[ "ABC" 6 [ dup >string next-permutation ] times ] unit-test

This is available in the latest development branch of Factor.

Bonus

Solving the original programming challenge is as easy as:

IN: scratchpad 38276 number>string next-permutation string>number .
38627