Mathematica Pearls:
Problems and Solutions

Edited by Don T. Piele

Welcome back to Mathematica Pearls, the column devoted to examining interesting solutions to an assortment of small gems. You never know when a simple problem may produce a bushel of unusual solutions. The challenge for this issue – Making Change, solved by Michael Trott and He Yu, is a perfect example of a problem with a bountiful yield. The authors approach the problem in three different ways; direct counting, recursion, and generating functions. Generating functions are broken down further and a solution emerges using partial fractions, which instantly returns a solution for the big bank rolls, such as Bill Gates (18 (10^11))and the Gross National Product (28 (10^14)) . Now that's a Mathematica Pearl.

Making Change

Write a function that will compute how many ways any amount of money m = 100 X + Y
 (with X dollars, and Y cents) can be represented using any combination of US coins and bills up to the 100 dollar bill.The denominations are: penny, nickel, dime, quarter, half dollar, dollar, 2 dollar, 5 dollar, 10 dollar, 20 dollar, 50 dollar, 100 dollar.

Solutions

by Michael Trott (mtrott@wolfram.com) and Yu He (yu@wolfram.com).

1. Direct Counting Method

Given an amount of money m = 100  X + Y (X dollars, Y cents), we wish to evaluate the number of ways in which it can be represented using the denominations ("coins") coins = {(Subscript[c, 1]), (Subscript[c, 2]), …, (Subscript[c, p])}. ((Subscript[c, i]) is expressed in units of cents.) That is, find the number of solution sets {K[1],  K[2],  …,  K[p]} (K[i] integer  ≥  0)  such that m = (∑_(i = 1)^p) K[i] (Subscript[c, i ]) (we use K as the summation variable because it is the default Mathematica summation variable). Although the order of (Subscript[c, i])s is arbitrary in the following programs, for specificity we suppose that (Subscript[c, i]) >  (Subscript[c, i + 1]). Then K[1] can assume the values 0, 1, …,  ⌊m/(Subscript[c, 1])⌋ (where ⌊x⌋ represents the integer part of x) and  given K[1], K[2] can take the values 0, 1,  …, 
      ⌊(m - K[1] (Subscript[c, 1]))/(Subscript[c, 2])⌋, and so on. Therefore, the total number of ways   representing the money m, Change[m, coins],  is simply

(∑_(K[1] = 0
          )^⌊)m/(Subscript[c, 1])⌋(∑_(K[2] = 0
            )^⌊)(m - K[1] (Subscript[c, 1]))/(Subscript[c, 2])⌋… 
          (∑_(K[p - 1] = 0
                )^⌊)(m - (∑_(i = 1)^(p - 1))K[i] (Subscript[c, i]))
                    /(Subscript[c, p - 1])⌋1
                                    
                                            
              (1)

This can be implemented quite easily with Mathematica:

\(Change1[m_, coins_List] := 
	
      \(Fold[Sum, Sum[1, #1], {##2}]&\)@@
        MapThread[{#1, 0, #2}&, 
          \({K/@#1, 
                ⌊\((
                      m - Reverse[
                          FoldList[#2\ K[#2] + #1&, 0, Reverse[Rest[#1]]]])
                      \)/#1⌋}&\)[Rest[Sort[coins]]]]\)

We can see that Change1 gives (1) by giving a symbolic m and a list of coins.

Change1[m, {1, 5, 25, 10, 50}]

Change1[100, {100, 50, 25, 10, 5, 1}] // Timing

In Mathematica, one can specify the increment ΔK of the summation variable K. So instead of summing K[1] from 0 to ⌊m/(Subscript[c, 1])⌋ we can sum it from 0 to m in steps of ΔK[1] = (Subscript[c, 1]), and so on. This feature of Mathematica is exploited below to speed up and slightly shorten the above implementation.

Change1a[m_, coins_List]:=
Sum[1,##]&@@Drop[MapThread[{#1,0,#2,#3[[1]]}&,{#,m-FoldList[Plus,0,Drop[#,-1]],#}&[K /@ #]]&[Reverse[Sort[coins]]],-1]

Change1a[m, {1, 5, 25, 10, 50}]

Change1a[100, {1, 5, 25, 10, 50, 100}] // Timing

One advantage of this direct counting approach is that it can be used to generate an explicit list of all coin combinations. We simply change Sum to Table and prepend K[p], the number of one cent coins, since it is not explicitly computed in the above programs.

\(Change1Explicit[m_, coins_List] := 

      Module[{s = Sort[coins]}, 

        Flatten[\(Table[
                \(Prepend[K/@#, m - Plus@@\((K/@#)\)]&\)[Rest[s]]/s, ##]&\)@@
            Drop[\(MapThread[{#1, 0, #2, #3[\([1]\)]}&, 
                    \({#, m - FoldList[Plus, 0, Drop[#, \(-1\)]], #}&\)[
                      K
/@\ #]]&\)[Reverse[s]], \(-1\)], Length[s] - 2]]\)

Short[Change1Explicit[100, {1, 5, 25, 10, 50, 100}], 12]

Here the first element of each sublist is the number of 1c coins, and the second the number of 5c coins, and so on. The last entry in each sublist is the number of 100c ($1) coins.

2. Recursive Method

Although the direct counting method given above is straigtforward, an iterative method may be more efficient, especially when repeated calculations are needed. Notice that the sum over K[1] in (1) can be divided into two parts: K[1] = 0 and K[1] ≠ 0. The first part clearly gives the number of ways of representing the same amount of money m using coins = {(Subscript[c, 2]), (Subscript[c, 3, ]) …, (Subscript[c, p])} and the second part is precisely the number of ways of expressing m-(Subscript[c, 1])using all p types of coins. Mathematically, the recursion relation is

Change[m,  {(Subscript[c, 1]),  (Subscript[c, 2]),  (Subscript[c, 3]),  …,  (Subscript[c, p])}]= 
	Change[m - (Subscript[c, 1]),  {(Subscript[c, 1]),  (Subscript[c, 2]),  (Subscript[c, 3]),  …,  (Subscript[c, p])}]+Change[m, {(Subscript[c, 2]),  (Subscript[c, 3]),  …,  (Subscript[c, p])}]

This recursion is implemented below. Note that to speed up the calculation, we define Change2 such that it stores the values calculated in the intermediate steps.

Change2[m_?Negative,  _]  :=  0Change2[m_,   {_}]  :=  1Change2[m_,  coins_] := 
      (Change2[m, #] = 
              Change2[m, Drop[#, -1]] + Change2[m - Last[#], #])&[
        Sort[coins]]

Change2[100, {1, 5, 10, 25, 50, 100}] // Timing

Using this recursive method we can treat much larger amounts of money than we could using Change1. We calculate below the number of ways to change a  $100 bill (i.e., m = 10000).

$RecursionLimit = Infinity;

\(Change2[10000, 

        	\ {1, 5, 10, 25, 50, 100, 500, 1000, 2000, 5000, 10000}] // Timing
      \)

We see that we have computed the solutions to many intermediate values of m.

Length[DownValues[Change2]]

3. Methods Using the Generating Function

In this section we use the standard method of generating functions to solve the problem. We introduce a complex variable z and imagine (Subscript[(z^c), i]) represents one coin of value (Subscript[c, i]) and ((((Subscript[(z^c), i])))^2) = (z^(2  (Subscript[c, i]))), two such coins, and so on. Then, the total number of ways of  representing the money m is simply the coefficient of (z^m) in the following expression (see Pólya (1956)):

\(\((1 + z\^c\_1 + z\^\(2  c\_1\) + z\^\(3  c\_1\) + …)\) 
        \((1 + z\^c\_2 + z\^\(2  c\_2\) + \(z\^\(3  c\_3\)\) …)\) 
        … 
        \((1 + z\^c\_p + z\^\(2  c\_p\) + z\^\(3  c\_p\) + …)\) = 

      	\(1\/\(\((1 - z\^c\_1)\) \((1 - z\^c\_2)\) … 
            \((1 - z\^c\_p)\)\) = 
        ∑\+\(m = 0\)\%∞\( A\_m\) 
            z\^m\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
            \((2)\)\)\)

The left hand of (2) is called the generating function  of the problem since the coefficients of its series expansion give the desired results. In principle, we can of course get the coefficient (Subscript[A, m ])using Series. For example, for coins = {1, 5, 10,  25,  50, 100} and m = 100, we have,

Timing[
      SeriesCoefficient[
        Series[(1/((1 - z) (1 - (z^5)) (1 - (z^10)) 
              (1 - (z^25)) (1 - (z^50)) (1 - (z^100)))), {z, 0, 
            100}], 100]]

However, expanding the generating function to (m^th) power can be slow when m is large. In the following we give other ways of calculating the coefficient (Subscript[A, m]).

3.1 Recursive Method

Multiplying both sides of (2) by (1 - (Subscript[(z^c), 1])) (1 - (Subscript[(z^c), 2])) … (1 - (Subscript[(z^c), p])), we have

(1 - (Subscript[(z^c), 1])) (1 - (Subscript[(z^c), 2]))... (1 - (Subscript[(z^c), p])) 
        (∑_(k = 0)^∞) (Subscript[A, k]) (z^k)  = 1

The coefficient of (z^m)should all be zero except for m = 0 which is 1. By collecting the terms proportional to (z^k) and setting its coefficient to zero, we can write a recursive relation between the  (Subscript[A, k])s. For example, for coins = {1, 5}, we have (1 - z) (1 - (z^5)) = 1 - z - (z^5) + (z^6) and the recursion relation is \(TraditionalForm
      \`A\_\(m\ \)\  = \ A\_\(m - 1\) + A\_\(m - 5\) - \(A\_\(m - 6\).\)\) Here is an implementation of this method.

\(Change3[m_, coins_] := 
	
      \(Module[{
            coef = \(-
                Reverse[
                  Rest[CoefficientList[Expand[Times@@\((1 - z\^coins)\)], z]]]
                  \)}, Nest[Append[Rest[#1], coef.#1]&, 
            Append[Table[0, {Length[coef] - 1}], 1], m]]
          \)[[\(-1\)]]\)

Change3[100, {1, 5, 10, 25, 50, 100}] // Timing

As is done in the earlier recursive method, one can write the program such that it stores the intermediate values if repeated calculation is required.

3.2 Direct Integration

Consider integrating the function

f (z) = 
      ((z^(-(m + 1)))/((1 - (Subscript[(z^c), 1])) (1 - (Subscript[(z^c), 2]))... 
            (1 - (Subscript[(z^c), p])))) = 
        (∑_(k = 0)^∞) (Subscript[A, k]) 
            (z^(k - (m + 1)))              
                     (3)

around a circle of radius r (r < 1). Since

∫(z^n) \[DifferentialD]z=(r^(n + 1 )) 
          (∫_0
              ^(2  π))(e^(I (n + 1) ϕ)) 
              \[DifferentialD]ϕ={({{0         n ≠ -1}, {2  π I   n = -1}})

we see that the integral of (3) will give 2  π i (Subscript[A, m]). So we can simply integrate the left hand side of (3)   with, say, r = 1/2 using NIntegrate and divide the result by 2 π i to get (Subscript[A, m]). For m = 100 and coins = {1, 5, 10, 25, 50, 100},  we have

\(Change4[m_, coins_] := 
	
        Rationalize[
          \(1\/\(2\ π\ I\)\) 
            NIntegrate[
              Evaluate[
                1\/2\ \((
                    z\^\(\(-m\) - 1\)\/Times@@\((1 - z\^coins)\) /. 
                      z -> 1\/2\ Exp[I\ phi])\)\ I\ Exp[I\ phi]], {phi, 
                0, 2\ π}, Method -> Trapezoidal, 
              Compiled -> False, WorkingPrecision -> 50, 
              PrecisionGoal -> 16, MaxRecursion -> 10]]\)

Change4[100, {1, 5, 10, 25, 50, 100}] // Timing

For this example this method seems overly complicated and  slow. However, for more complicated enumeration problems in combinatorics this turns out to be quite a useful method, especially when the generating function is not a rational fucntion, as is the case in the problem under consideration.

3.3 Residue method

In fact, (Subscript[A, m]) can be obtained without having to do the above integral. From the theory of complex functions, the integral is 2  π i times the residue at z = 0. Therefore, (Subscript[A, m]) is precisely the residue of the function f(z) at z = 0. For example, for the same m and coins, we can do

Timing[
      Residue[(1/(z^101)) 
          (1/( (1 - z) (1 - (z^5)) (1 - (z^10)) 
              (1 - (z^25)) (1 - (z^50)) (1 - (z^100)))), {z, 0}]]

Or, in general, we can use the following program:

Change5[m_, coins_] := 
      Residue[((z^(-m - 1))/Times)@@(1 - (z^coins)), {z, 0}]

However, when m is large and there are more coins the different methods described above are usually slow. It is desirable to have a closed-form expression for (Subscript[A, m])s as a function of m, given coins, so that once the expression is calculated, we can just substitute in the value of m and get the result.

Unfortunately, the method of calculating the residue at z = 0 cannot be used for this purpose. This is because Residue cannot be used to get the residue at a pole whose order is not numeric. Again we resort to complex function theory to circumvent this difficulty. Since the integral of f(z) along a circle of radius r ->  ∞ vanishes, the sum of the residues inside the circle is zero. Therefore, the residue at z = 0 is equal to the negative of the sum of the residues of poles on the unit circle. Unlike z = 0, the orders of the poles on the unit circle are independent of m, hence, they can be calculated using Residue.

For example, take coins = {1, 5, 10}. The poles on the unit circle are

\(coins = {1, 5, 10}; 
    poles = Union[
        Flatten[\(Table[E\^\(\(I\ 2\ π\ k\)\/#\), {k, 0, #}]&\)/@coins]]
        \)

These are the residues from the two real poles:

\(f = 
      z\^\(\(-m\) - 1\)\/\(\((1 - z)\)\ \((1 - z\^5)\)\ \((1 - z\^10)\)\); 

    r1 = Plus@@\(\((Residue[f, {z, #1}]&)\)/@{1, \(-1\)}\)\)

Since the residues of a pair of complex conjugate poles are complex conjugate, we only need to calculate the residue of one of the pair and add its complex conjugate. r2 gives the contribution from all complex poles:

r2 = 
      Plus@@(
            (Numerator[#1]/Together[Expand[Denominator[#1]]]&)[
                Together[
                  (#1 + (#1 /. Complex[a_, b_] -> Complex[a, -b])
                            &)[
                    Residue[f, {z, #1}] /. 
                      (((-1))^x)_ -> 
                        Cos[π x] + I Sin[π x]]]]&)/@
          Select[poles, Im[#1] > 0&];

The desired formula is simply the negative of r1 + r2

A[m_, {1, 5, 10}] = Simplify[-r1 - r2]

A[3458472, {1, 5, 10}]

This must be an integer; using RootReduce we can simplify the above expression to an explicit integer.

RootReduce /@ %

The following program gives the formula  A[m,  coins] for general coins. However, when the number of coins is large, the expression is very involved and takes long time to simplify. (Note that given coins, Change6 gives an explicit expression that is a function of m.)

\(Change6[m_, coins_] := 
	
      Module[{f = z\^\(\(-m\) - 1\)\/Times@@\((1 - z\^coins)\), RealPoles, 
          ComplexPoles}, {RealPoles, ComplexPoles} = 
          \(({Cases[#, 1 | \(-1\)], DeleteCases[#, 1 | \(-1\)]}&)\)[
            Union[Flatten[
                \((Table[
                        E\^\(\(I\ 2\ π\ k\)\/#\), {k, 0, 
                          ⌊#\/2⌋}]&)\)/@coins]]]; 
        \(-Plus@@\(\((Together[Residue[f, {z, #}]]&)\)/@RealPoles\)\) - 
          Plus@@\(\((
                \((Numerator[#]\/Together[Expand[Denominator[#]]]&)\)[
                    Together[
                      \((# + \((# /. 
                                  Complex[a_, b_] -> Complex[a, \(-b\)])
                                \)&)\)[
                        Residue[f, {z, #}] /. 
                          \((\(-1\))\)\^x_ -> 
                            Cos[π\ x] + I\ Sin[π\ x]]]]&)\)/@
              ComplexPoles\)]\)

3.4 Partial Fraction Method

In this section we will give yet another method to get a formula for (Subscript[A, m]) as a function of m. We will restrict our discussion to coins = {1, 5, 10, 25, 50, 100}. Using more coins is possible in the approach shown below, but it will dramatically increase the size of some of the expressions as well as the time needed to calculate these expressions.

The generating function (2)

s0 = 
      (1/((1 - z) (1 - (z^5)) (1 - (z^10)) (1 - (z^25)) 
          (1 - (z^50)) (1 - (z^100))));

can be written as

(∑_(i = 1)^p)
        (∑_((Subscript[j, i])))((q 
              ((Subscript[c, i]),  (Subscript[j, i])))/(Subscript[(((1 - (Subscript[(z^c), i])))^j), i]))
                                     
                                                 
                                               (4)

where q((Subscript[c, i]), (Subscript[j, i])) is a polynomial in z. Once the generating function is cast in this form, we can expand the denominator using the binomial formula

(1/(((1 ± (z^k)))^j))=(∑_(n = 0)^∞)(({{-j}, {n}}))(((∓(z^k)))^n)

and then collect the coefficient of (z^m), (Subscript[A, m]). As a first step to obtaining the sum of the above form we apply Apart to the expression s0.

s1 = Apart[s0]

Unfortunately the denominators are not of the form we want. But because the denominators arose from factoring  \(TraditionalForm
      \`\(\((1 - z)\) \((1 - z\^5)\) \((1 - z\^10)\) \((1 - z\^25)\) 
        \((1 - z\^50)\) \((1 - z\^100)\)\ \)\), we can multiply each numerator and denominator by a suitable polynomial to get the desired form. The function toNiceDenominator transforms all fractions to the desired denominator form given in (4).

\(toNiceDenominator[p_\ q_Plus?\((Length[#1] > 2&)\)\^n_Integer?Negative]
         := 
	Module[{
          t = \(Select[{1 - z\^5, 1 - z\^10, 1 - z\^25, 1 - z\^50, 
                  1 - z\^100}, \ \ \ 

                					\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
                \ \ \ \ \ \ \ \ \ \ \ \ \ \ PolynomialQ[Cancel[#/q]]&, 1]
                \)[[1]], f}, 
        f = Cancel[t/q]; FactorTerms[Expand[f\^\(-n\)\ p]]\ t\^n]
\)\(toNiceDenominator[p_\ \((z - 1)\)\^n_Integer?Negative] := 
      p\ \((\(-1\))\)\^n\ \((1 - z)\)\^n
\)toNiceDenominator[p_ ((((z^l)_. + 1))^n)_Integer?Negative] := 
      p (((1 + (z^l)))^n)

Short[s2 = toNiceDenominator/@s1, 8]

Now it is straightforward to convert the fractions into binomials:

\(toBinomial[q_\ z\^k_.\ \((1 + p_.\ z\^m_.)\)\^l_.] = 
      q\ F[\(n - k\)\/m, l, \(n - k\)\/m, p]; 

    toBinomial[q_\ \((1 + p_.\ z\^m_.)\)\^l_.] = q\ F[n\/m, l, n\/m, p]; \)

Short[s3[n_] = toBinomial/@Expand[s2], 8]

The result is not a small expression, it contains 393 summands:

Length[s3[n]]

The explicit form of the function F is given by:

F[x_, y_, z_, p_] := 
      If[IntegerQ[x] && NonNegative[x], 
        (((-1))^x) Binomial[x - y - 1, x] 
          (((-1))^((1/2) z (1 - p))), 0]

Our standard example gives again the answer 293:

s3[100]

But now we can calculate (Subscript[A, m ])for much larger m:

Timing[s3[(10^6)]]

For the Bill Gates order of magnitude

Timing[s3[18 (10^11)]]

and the US GNP order of magnitude:

Timing[s3[28 (10^14)]]

4. The "Super-Rich" (m->∞) Limit

In this section we explore how (Subscript[A, m]) grows in the m ->  ∞ ("super-rich" limit). From the residue method, we know that (Subscript[A, m]) is given by the sum of the residues from the poles of f(z) on the unit circle (with a minus sign). Since the residue from a (k^th) order pole (Subscript[z, 0]) is (1/((k - 1)!)) (Underscript[limit, z -> (Subscript[z, 0])])((d^(k - 1))/(dz^(k - 1)))((f(z) (z - (Subscript[z, 0])))^k)  and since each derivative will bring down a factor of m, the leading term in large m limit comes from the highest order pole (Subscript[z, 0]) = 1. It is straightforward to see from the residue   formula that the leading term is

((m^(p - 1))/((p - 1)! (Subscript[c, 1]) (Subscript[c, 2.]).. (Subscript[c, p])))

For example, for coins = {1, 5, 10} the leading term is

((m^(3 - 1))/((3 - 1)! 5 10))

as is shown in the expression for A[m, {1, 5, 10}] above.

Apart from the leading term, the terms proportional to (m^k) (k < p - 1) have coefficients that are oscillatory (sines and   cosines). It is this oscillatory contribution that gives (Subscript[A, m])=(Subscript[A, m + 1 ]) = 
            …  = (Subscript[A, m + (Subscript[c, p - 1]) - 1]) (where m is divisible by (Subscript[c, p - 1])). Here we plot (Subscript[A, m]) as a function of m for relatively small m, so we can see the oscillatory behavior.  (Note that we plot the real part of A[m, {1, 5, 10}] since it is complex when m is not an integer.)

\(\(Plot[Evaluate[Re[A[m, {1, 5, 10}]]], {m, 1, 80}, 
	\ \ \ \ \ \ \ 
      PlotPoints -> 80, 
      Epilog -> {PointSize[0.005], Hue[0], \ 
          Table[Point[{m, Change2[m, {1, 5, 10}]}], {m, 1, 80}]}]; \)\)

The imaginary part vanishes identically at integer m:

Plot[Evaluate[Im[A[m, {1, 5, 10}]]], {m, 1, 40}];

This is a plot of (Subscript[A, m]) for m = 1 to m = 5000.

Plot[Evaluate[Re[A[m, {1, 5, 10}]]], {m, 1, 5000}];

When m is large, the term proportional to m\(\^\(p - 1\)\) dominates. For m = 5000 the leading term yields

((5000^2)/100)

and the exact expression gives

A[5000, {1, 5, 10}]  // Expand

Here the relative error between the leading asymptotic term and the exact result is shown:

Plot[
      Evaluate[1 - ((m^2)/(100 Re[A[m, {1, 5, 10}]]))], {m, 1, 5000}];

References

The problem of calculating the number of ways to  change a 1$ goes back to Pólya: Pólya, G. American Mathematics Monthly 63,  689 (1956).

The recursice approach to the problem is discussed in: 
R. L. Graham. D. E. Knuth, O. Patashnik. Concrete Mathematics, Addison-Wesley, Reading 1994.

Various methods for the asymptotic evaluation of coefficients of generating functions are discussed in:
A. M. Odlyzko in R. Graham, M. Grötschel, L.   Lovász (eds.). Handbook of Combinatorics, Elsevier, Amsterdam 1995.

Psycho Dice

The challenge for the next time involves the game Psycho Dice which is 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."  Write a function, psychoDice[n]
which will compute the probability of winning this game when played with n dice and n numbers.

SUBMITTING PROBLEMS AND SOLUTIONS

Have you got a Mathematica Pearl buried deep in your  files?  I invite you to share them in this column and add them to our string. If you have the goods, then "Show me the pearls!" The best way to get my attention is to send your work via email to <piele@cs.uwp.edu>. Little pearls slip nicely into an email message you can zip to me from anywhere. You will get immediate feedback. Accompany your work with an explanation of how it was created. That's it, and thanks.

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

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file Pearls62.nb containing Mathematica code for the material described in this article.