Mathematica Pearls

Problems and Solutions

by Don Piele

Welcome back to Mathematica Pearls, the column devoted to examining interesting solutions to an assortment of simple yet appealing problems. The problem for this issue appeared in the delightful book Programming Pearls by Jon Bentley,  Addison-Wesley, 1986.  Here he presents this issue's problem as a lesson in finding dramatic performance improvements in computer algorithms.

The Problem

Given a sequence {(Subscript[x, 1]),  (Subscript[x, 2]), ....(Subscript[x, n])}  of integers of length n, write an efficient algorithm (of order n) to find a subsequence {(Subscript[x, L]), (Subscript[x, L + 1]),....(Subscript[x, H])} of consecutive terms with the largest sum. Return the beginning and ending indices L, H and the maximum sum  (∑_(i = L)^H) (Subscript[x, i]). If there is more than one such subsequence, return any one.

Order (n^3) solution

Brute force can be used to try all possibilites. For each lower index L from 1 to n, and each higher index H from L to n, sum the subsequence {(Subscript[x, L]), (Subscript[x, L + 1]),....(Subscript[x, H])}, and see if it's the best so far. If it is, save it  along with the indices L and H.

\(bruteForceOne[X_] := \ 
      Module[
		{MaxSoFar = 0, L = \(H = 0\), Low = \(High = 0\), 
          SubSub = 0, 
		n = Length[X], j}, 
		
        For[L = 1, L <= n, \(L++\), 
			
          For[H = L, H <= n, \(H++\), SubSum = Sum[X[\([j]\)], {j, L, H}]; 

            				If[SubSum >= MaxSoFar, Low = L; High = H]; 
				
            MaxSoFar = Max[MaxSoFar, SubSum]
				]
		]; 

        		{MaxSoFar, {Low, High}}
	]\)

bruteForceOne[winnings] // Timing

For large values of n this solution slows down dramatically since the number of computer operations (hence the time) is proportional to (n^3.)  One n for each loop and one n for the sum computation. To get an estimate for the time it would take to find an answer for a sequence of length 1,000 consider  the following.

Clear[timeInSeconds, timeInHours]solution = 
      Solve[ (timeInSeconds/((1000^3) )) == (( .44 seconds)/(30^3)), 
        timeInSeconds]

timeInSeconds = timeInSeconds  /.  First[solution]

timeInHours =  (timeInSeconds/(3600  (seconds/hours)))

Order (n^2) solution

We can cut it down by one power of n by calculating  the sum  as we move along. Here is the algorithm given by Bentley which runs in time proportional to (n^2.)

\(bruteForceTwo[X_] := \ 
      Module[
		{MaxSoFar = 0, L = \(H = 0\), Low = \(High = 0\), 
          SubSub = 0, 
		n = Length[X], j}, 
		
        For[L = 1, L <= n, \(L++\), 
			SubSum = 0; 
			
          For[H = L, H <= n, \(H++\), SubSum = SubSum + X[\([H]\)]; 
				
            If[SubSum >= MaxSoFar, Low = L; High = H]; 
				
            MaxSoFar = Max[MaxSoFar, SubSum]
				]
		]; 

        		{MaxSoFar, {Low, High}}
	]\)

bruteForceTwo[winnings] // Timing

To get an estimate for the time it would take to find an answer for a sequence of length 1,000 consider  the following.

Clear[timeInSeconds]Solve[ (timeInSeconds/((1000^2) )) == (( .11 seconds)/(30^2)), 
      timeInSeconds];

(timeInSeconds/(60  (seconds/minutes)))  /.  First[%]

The crème de la crème solution.

The following graph shows how timings grow with n for all  algorithms from the least efficient to the most efficient.

Your problem is to find an order n (or crème de la crème) solution that will run in less than one second for n equal to 1,000.

A Previous Problem

Psycho Dice

The challenge that appeared in  Vol 6 No 2 involved the game Psycho Dice which was described in the book, Midnight In The Garden of Good and Evil, by John Berendt. On page 21 it reads, " I've invented a game called Psycho Dice. It's very simple. You take four dice and call out four numbers between one and six – for example,  {4,3,6,6}.Then you throw the dice, and if any of your numbers come up, you leave those dice standing on the board.  You continue to roll the remaining dice until all the dice are sitting on the board, showing your set of numbers. You're eliminated if you roll three times in succession without getting any of the numbers you need."  Compute the probability of winning this game when played with n dice.

Solution
by Robert Dickau (pops@prairienet.org)

Psycho Coins

I started  out by computing the probabilities for Psycho Coins instead of Psycho Dice.  Also,  it seemed easier to figure out the probability of losing instead of winning.  Without any loss in generality, we can assume that we have called all  "heads."

One Coin Game

We know the probability of losing the one-coin game—getting three tails in a row—is ((((1/2)))^3) = (1/8). I'll write the probability of losing the 1-coin game as (Subscript[Pr, 1])(Lose).

Two Coin Game

One way to lose the two-coin game is to get zero heads on any of the first three tries, which happens with probability ((((1/4)))^3) = (1/64).

The other ways to lose are to get one head on the first toss, then throw three tails; or get one head on the second toss, followed by three tails; or get one head on the third toss, followed by three heads. If we denote the probability of getting m heads on the  (t^th) toss (with 2 coins) by (Subscript[Pr, 2])(m, t), the total probability of losing the 2-coin game is

\(TraditionalForm
    \`\(Pr\_2\)(Lose) = 
      \(\(Pr\_2\)(1, 1)\) \(\(Pr\_1\)(Lose)\) + 
        \(\(Pr\_2\)(1, 2)\) \(\(Pr\_1\)(Lose)\) + 
        \(\(Pr\_2\)(1, 3)\) \(\(Pr\_1\)(Lose)\) + \(Pr\_2\)(0, 3)\)

Three Coin Game

Similarly, we can lose the three-coin game by getting one head on the first, second, or third toss, followed by losing the two-coin game; or getting two heads in one of the first, second, or third tosses, followed by losing the one-coin game; or throwing all tails on each of the first three throws. The formula would then be:

\(TraditionalForm
    \`\(Pr\_3\)(Lose) = 
      \(\(Pr\_3\)(1, 1)\) \(\(Pr\_2\)(Lose)\) + 
        \(\(Pr\_3\)(1, 2)\) \(\(Pr\_2\)(Lose)\) + 
        \(\(Pr\_3\)(1, 3)\) \(\(Pr\_2\)(Lose)\) + 
        \(\(Pr\_3\)(2, 1)\) \(\(Pr\_1\)(Lose)\) + 
        \(\(Pr\_3\)(2, 2)\) \(\(Pr\_1\)(Lose)\) + 
        \(\(Pr\_3\)(2, 3)\) \(\(Pr\_1\)(Lose)\) + \(Pr\_3\)(0, 3)\)

Putting this together gives this:

Pr[1][ Lose] = ((((1/2)))^3);

Pr[coins_][heads_, try_] := 
      Binomial[coins, heads] (((((((1/2)))^coins)))^try)

Pr[coins_][Lose] := 
      Pr[coins][Lose] = 
        (∑_(i = 1)^(coins - 1))(∑_(j = 1)^3)
                 Pr[coins][i, j] Pr[coins - i][Lose] + 
          Pr[coins][0, 3]

The probability of losing the two-coin game  is:

Pr[2][Lose] // {#, N[#]}&

And the probability of winning is:

1 - %

The probability of winning a three coin game is:

Pr[3][Lose] // {1 - #, N[1 - #]}&

More concisely, the probabilities of winning with n coins are approximately

1.0 - Pr[#][Lose]&/@Range[25]

PsychoDice

Now we turn our Psyco Coins into Psycho Dice. Assume, without any loss in generality,  that we call out all "ones."

One Die Solution

The probability of losing the one-die game is the following.

PrDice[1][ Lose] = ((((5/6)))^3);

PrDice[dice_][ones_, try_] := 
      Binomial[dice, ones] ((((1/6)))^ones) ((((5/6)))^(dice - ones)) 
        (((((((5/6)))^dice)))^(try - 1))

PrDice[dice_][Lose] := 
      PrDice[dice][Lose] = 
        (∑_(i = 1)^(dice - 1))(∑_(j = 1)^3)
                 PrDice[dice][i, j] PrDice[dice - i][Lose] + 
          PrDice[dice][0, 3]

Here are the probabilities of winning Psycho Dice with 1 through 25 dice.

{#, 1.0 - PrDice[#][Lose]}&/@Range[25] // TableForm

The probability converges to .240548.  For the orginal four coin PsychoDice, it's slightly above 1/4 and clearly not a very fair game.

Simulation Solution (editor)

It's always a good idea to check out a solution like this one with a simulation -- since it is so easy and natural to do. The function PsychoDice[x] plays one game with x dice and  returns a 1 for a win or a zero for a loss.

\(\(
PsychoDice[dice_] := 
      Module[{win = 0, \ diceLeft\  = dice, tossesLeft = 3, roll}, 

        \ \ \ \ \ \ \ 
        While[tossesLeft > 0\  && \ win == 0, 
	\ \ \ \ \ 
          roll = Table[Random[Integer, {1, 6}], {diceLeft}]; 
			
          hits = Count[roll, 1]; 
		\ \ \ 
          If[hits > 0, tossesLeft = 3, \(tossesLeft--\)]; 
	\ \ \ \ \ 
          diceLeft = diceLeft - hits; 
	\ \ \ \ \ \ 
          If[diceLeft == 0, win = 1]]; 
\ \ \ \ \ win]\)\)

PsychoDiceSim[x,y] plays x games of PsychoDice with y dice.

PsychoDiceSim[games_, dice_] := 
      N[Plus@@ Table[PsychoDice[dice], {games}]/games]

Here is 1,000 games played with 4 dice.

PsychoDiceSim[1000, 4]

This agrees with  Robert's calculations.

Finally, we play the game 1,000 times for  1 through 10 dice.

{#, PsychoDiceSim[1000, #]}&/@Range[10] // TableForm

Show me the pearls!

Have you got a Mathematica Pearl buried deep in your files?

About The Editor

Don Piele has been interested in creating programming problems since he began the
International Computer Problem Solving Contest in the pages of Creative Computing
 in 1981. In 1992, he organized the USA Computing Olympiad, which selects the top
 four American high school computer programmers to represent the USA at the
annual International Computing Olympiad. He also writes the column, Cowculations,
devoted to computer algorithms using Mathematica in Quantum Magazine. The Web
 address for the USACO is: http://usaco.uwp.edu.

Don Piele 
Mathematics Department
University of Wisconsin-Parkside
Kenosha, WI 53141
piele@cs.uwp.edu