Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Need a Faster Solution to Number Theory Problem

• To: mathgroup at smc.vnet.net
• Subject: [mg91402] Re: [mg90990] Need a Faster Solution to Number Theory Problem
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Wed, 20 Aug 2008 04:04:48 -0400 (EDT)
• References: <200808020726.DAA24898@smc.vnet.net>

John Snyder wrote:
> I got hooked on the following little number theory problem which appears in
> the August issue of Discovery magazine.  "Susan has more than one dollar in
> change, but she still can't make exact change for a dollar bill. What is the
> maximum amount of money that Susan could have in pennies, nickels, dimes and
> quarters for this to be true?"
>  
> I solved the problem using the following Mathematica code:
>  
> coins={1,5,10,25};
> dollar=FrobeniusSolve[coins,100];
>  
> Catch[Do[t=Times@@Length/@(Cases[dollar,{p_,n_,d_,q_}/; (p<=#1&&n<=#2 &&
> d<=#3 && q<=#4)]&@@@FrobeniusSolve[coins,a]);
> If[t==0,Throw[a]],{a,130,101,-1}]]
>  
> Starting with a maximum guess of $1.30 this gives the correct answer of
> $1.19 in about 5 seconds.  The problem is that if the problem is made more
> complicated by adding more possible coin denominations (and/or paper
> currency) the execution time become very long indeed.
>  
> Is there a way to speed up the solution to this type of problem so that I
> can handle more complicated cases in a reasonable time?
>  
> Thanks,
>  
> John

Late response, but I didn't get time to look hard at this until now. It 
seems to be a difficult problem in terms of complexity, though I'm not 
certain I have the best possible approach.

First the uninteresting case: if one denomination does not divide the 
avoided amount then the answer is infinity. So we'd want to rule out 
that case, and just handle the rest under the assumption that all 
denominations divide the value we need to avoid.

Here, in brief, is one way to go about this. The set of solutions to the 
Frobenius instance problem carve out a "staircase" in the denominations. 
Any set giving a maximum value must always have an element below the 
"elbows" of this staircase. In effect this gives a domain quite similar 
to the "fundamental domain" of the Frobenius number problem, as 
discussed here:

http://www.integers-ejcnt.org/vol7.html
(look for A15: Frobenius Numbers by Lattice Point Enumeration)

or

http://citeseerx.ist.psu.edu/viewdoc/summary;jsessionid=A56BBBFE608FF49E382BE25E84F2EBC3?doi=10.1.1.88.9057

Given the elbows, we need to find a maximal corner, with l_1 norm (that 
is, Manhattan distance) weighted by the set of denominations. There are 
various ways to do this. I'll give code for one such below (it's adapted 
from code related to the paper at the URLs above). Other ways involve 
ideas from the theory of monomial ideals; see work by Bjarke Roune for 
details on efficient approaches.

Explaining the code below would be an article in itself, so I defer to 
the reference. I simply mention that there might be ways to simplify it 
because the original was written for a slightly different problem, and I 
may have missed modifications that could make it better for the problem 
now at hand. Suffice it to say that the slow step is in tracking through 
all the corners to find the largest. Possibly there is a smart way to 
improve on this.

dominationKernel[X_] :=
  Module[{Y = X[[Ordering[Map[Tr, X]]]], Z, i = 2, len = 1, k},
   Z = Table[{}, {Length[Y]}];
     Z[[1]] = Y[[1]];
     While[i <= Length[Y] && Tr[Y[[i]]] == Tr[Y[[1]]], len++;
       Z[[len]] = Y[[i]]; i++];
     Do[k = 1;
       While[k <= len,
         If[And @@ Thread[Z[[k]] <= Y[[j]]], k = len + 2; Break[],
      k++]];
       If[k == len + 1, len++; Z[[len]] = Y[[j]]], {j, len + 1,
     Length[Y]}];
     Sort[Take[Z, len]]]

ClearNegsAndDeleteZeroVector[vecs_] :=
  If[vecs == {}, {},
   Union[DeleteCases[vecs /. _?Negative :> 0,
     Table[0, {Length[vecs[[1]]]}]]]];

farthest[corns_, A_] :=
  Fold[If[#2.A > #1[[2]], {#2, #2.A}, #1] &, {Table[0, {Length[A]}],
     0}, corns][[1]]

FarthestCorner[A_, elbows_?MatrixQ] :=
   Module[{pts, p1, p2, B, cc, trQ, farvertex, maxes, currfar},
     currfar = farthest[elbows, A].A;
     pts = Select[elbows, #[[1]] > 0 &];
     Scan[({p1, p2} = {First[#], Rest[#]};
        B = Rest /@ Select[elbows, #[[1]] < p1 &];
        B =
         dominationKernel[
          ClearNegsAndDeleteZeroVector[# - p2 & /@ B]];
        maxes = Prepend[Max /@ Transpose[B], 0];
        If[(maxes + #).A > currfar,
         cc = FarthestCornerSub[A, B, currfar, #];
         farvertex = farthest[Prepend[# + p2, p1] & /@ cc, A];
         currfar = Max[currfar, farvertex.A];
         ]) &, pts];
     {currfar - Total[A], farvertex - 1}] /; Length[elbows[[1]]] > 2;

FarthestCorner[A_, elbows_?MatrixQ] :=
  farthest[FarthestCornerSub[A, elbows], A].Rest[A] - Total[A] /;
   Length[elbows[[1]]] == 2

FarthestCornerSub[A_, elbows_, currfar_, backdata_] :=
   Module[{p1, p2, pts, B, maxes},
     pts = Select[elbows, First[#] > 0 &];
     Flatten[Map[(
         {p1, p2} = {#[[1]], Rest[#]};
         B = Rest /@ Select[elbows, First[#] < p1 &];
         B = dominationKernel[
           ClearNegsAndDeleteZeroVector[# - p2 & /@ B]];
         If[B == {}, {},
          maxes = Prepend[Max /@ Transpose[B], 0];
          If[(PadLeft[maxes + #, Length[A]] + backdata).A > currfar,
           Prepend[# + p2, p1] & /@
            FarthestCornerSub[A, B, currfar,
             PadLeft[#, Length[A]] + backdata], {}]
          ]) &, pts], 1]] /; Length[elbows[[1]]] > 2;

FarthestCornerSub[_,
    elbows_, ___] := (Partition[
      Take[Flatten[Reverse /@ Reverse[elbows]], {2, -2}], 2]) /;
    Length[elbows[[1]]] == 2;

maxNoChange[coins_List, amount_Integer /; amount >= 1] /;
   Not[And @@ Map[IntegerQ, amount/coins]] := Infinity

maxNoChange[coins_List, amount_Integer /; amount >= 1] :=
  FarthestCorner[coins, FrobeniusSolve[coins, amount]]

Your example:

In[230]:= maxNoChange[{1, 5, 10, 25}, 100] // Timing
Out[230]= {0.96806, {119, {4, 0, 4, 3}}}

We'll add one coin and up the amount.

In[231]:= maxNoChange[{1, 5, 10, 25, 50}, 150] // Timing
Out[231]= {13.9929, {169, {4, 0, 4, 1, 2}}}

It gets slow quickly, so to speak.

In[235]:= maxNoChange[{1, 5, 10, 25, 50}, 200] // Timing
Out[235]= {82.3371, {219, {4, 0, 4, 1, 3}}}

Daniel Lichtblau
Wolfram Research

• Follow-Ups:
  • Re: Re: Need a Faster Solution to Number Theory
    • From: Daniel Lichtblau <danl@wolfram.com>

• References:
  • Need a Faster Solution to Number Theory Problem
    • From: "John Snyder" <jsnyder@wi.rr.com>