Hoe maken we de Nederlandse portemonnee lichter?
Dit is een vraag die gesteld werd door Hans Wisbrun en mij erg bezig hield. Na een uitgebreide analyse kwam ik tot interessante ontdekkingen. Voor het tijdschrift Euclides van de Nederlandse Vereniging van Wiskundeleraren heb ik een beknopte, leesbaardere versie geschreven (jaargang 93/4).


Ik heb voor dit onderzoek een aantal scripts gebruikt, die ik via deze pagina deel. De Monte-Carlo berekeningen heb ik via Matlab gedaan. Het gaat hierbij vooral om de structuur en de genomen stappen - deze zijn eenvoudig in andere scripttalen over te nemen. Het bruteforcen van het vinden van muntstelsels met miniale aantal munten als wisselgeld zijn in C geschreven, voor een stelsel van 2 munten, en een variant voor een stelsel van 8 munten, waarin gebruik wordt gemaakt van multithreaded programmeren.


Alle code is volledig vrij te gebruiken naar eigen wens.













BoodschappenMonteCarlo.m

% Monte-Carlo boodschappensimulatie

% Dit zijn prijzen in cent van een sessie boodschappen bij een supermarkt
Prijzen = [102   99  338  278  338  169  259  259  693   57   39  149 1746  298  394 ...
           656  598   78   69  168  299  278   37  398  115  193  165  112  368  539 ...
            62  840  179  182  159   96   96  169  860   79  149  196  125   79  313 ...
           450  124  240  250  398  169  179  267  112  298  269   79  263  185  185 ...
           371  259  229  299  279  259  425   55   65   25  455  175  259  179   55 ...
           209   55  209  169  289  299  289  295  149  129  129  118];

% Aantal trekkingen
NrSamples = 100000;
Samples = zeros(NrSamples, 1);

% Simuleer een rondje boodschappen (trekking)
for j=1:NrSamples
    % Aantal stuks
    NrItems = randi(50, 1);
    
    % Trek NrItems aantal keer uit de mogelijke prijzen die ik heb (uniek)
    idx = randsample(length(Prijzen), NrItems);
    
    % Totaal aantal aan boodschappen voor trekking i
    Samples(j) = sum(Prijzen(idx));
    
    % Neem modulo 500, ofwel rest nadat de rest is betaald met papiergeld
    Samples(j) = mod(Samples(j), 500);
end

% En een plaatje (histogram met 1 cent binning)
figure
hist(Samples, 0:499);
xlim([0 500])
xlabel('wisselgeld (cent)');
ylabel('aantal keer');
title('Verdeling van het wisselgeld');






SlimBoodschappenMonteCarlo.m

% Monte-Carlo portemonnee

% Muntdistributie. Noot: voor 5 cent afronding schalen we naar 1-100
% Coins = [5 10 20 50 100 200]/5;
% MaxAmount = 100;
Coins = [1 2 5 10 20 50 100 200];
MaxAmount = 500;

MaxSamples = 10000;
NumCoins = zeros(1, MaxSamples);

% Portemonnee heeft alleen papiergeld
wallet = [];

for sample = 1:MaxSamples
    % We kopen iets van tussen de 1 en MaxAmount-1 cent
    price = randi(MaxAmount-1, 1);

    % We geven alles wat we in onze portemonnee hebben, en vragen de winkel
    % om het meest effici�nte wisselgeld
    amount = sum(wallet) - price;

    % ... en als het duurder is dan het muntgeld, leggen we paper bij
    amount = mod(amount, MaxAmount);

    % Het wisselgeld
    change = return_change(amount, Coins);

    % En dat gaat direct de portemonnee in
    wallet = change;

    % Aantal munten
    NumCoins(sample) = length(wallet);
end

% Maak een plaatje
figure
hist(NumCoins, min(NumCoins):max(NumCoins));
xlabel('aantal munten');
ylabel('aantal keer');
title(sprintf('Gemiddeld %.1f munten voor stelsel met %d munten.', mean(NumCoins), length(Coins)));






return_change.m

% Change function, returns also the coins used
% Stolen from: http://www.8bitavenue.com/2011/12/dynamic-programming-making-change/
function change = return_change(total, coins)

n = length(coins);
M = zeros(total+1, 1);
b = zeros(total, 1);

for j=1:total
    min = Inf;
    b(j) = coins(1);

    for i=1:n
        if (j-coins(i)>=0) && (M(j - coins(i) + 1) < min)
            min = M(j - coins(i) + 1);
            b(j) = coins(i);
        end
    end
    M(j + 1) = 1 + min;
end

change = [];
remainder = total;
while remainder>0
    change = [change b(remainder)];
    remainder = remainder - b(remainder);
end

return;






CoinBruteForce_2c.c

/* Optimal coin set computation by brute force, 2 coin variant */
/* Compile with e.g.:
    gcc -O3 CoinBruteForce_2c.c -o CoinBruteForce_2c
 * or:
    icc -O3 CoinBruteForce_2c.c -o CoinBruteForce_2c
    
 * Run as:
    ./CoinBruteForce_2c
 *   
 
   * Created by Sebastiaan Breedveld *
   http://sebastiaanbreedveld.nl/
 *
 * License: The Unlicense (http://unlicense.org/): do with it whatever you want.
 *
 */

#include <limits.h>
#include <stdio.h>

/* Maximum change is 499 cents */
/* Trick! If you change this to 99 cents, you are able to compute systems which start with 5 cents (and are multiple of 5 cents) */
#define MaxAmount 499

/* Define number of coins in this system */
#define NrCoins 2

/* uint32 should be sufficient and not too limiting */
typedef unsigned int INTTYPE;
#define INTMAX UINT_MAX;

INTTYPE MakeChangeProblem(const INTTYPE Total, const INTTYPE *CoinSet); 

int main(void)
{ 
    /* Inits */
    INTTYPE CoinSet[NrCoins] = {1, 2};
    INTTYPE TotalCoins, c1, MinAmount = INTMAX;
      
    for (c1=2; c1<=MaxAmount; c1++)
    {
        CoinSet[1] = c1;
        TotalCoins = MakeChangeProblem(MaxAmount, CoinSet);
        
        /* Only print sets which are better than found so far */
        if (TotalCoins <= MinAmount)
        {
            MinAmount = TotalCoins;
            printf("1 %3u %10u\n", c1, TotalCoins);
        }
    }
    return 0;
}

INTTYPE MakeChangeProblem(const INTTYPE Total, const INTTYPE *CoinSet)
{
    /* We use a static allocation here, which is faster than dynamic */
    INTTYPE MinCoins[MaxAmount+1];
    INTTYPE i, j, T;
        
    /* This part is taken from the web */
    /* http://codereview.stackexchange.com/questions/92811/find-the-minimum-number-of-coins/92820 */
    MinCoins[0] = 0;
    for (j=1; j<=Total; j++)
    {
        MinCoins[j] = INTMAX;
    }

    for (i = 0; i < NrCoins; i++)
    {
        for (j = CoinSet[i]; j <= Total; j++)
        {
            T = MinCoins[j - CoinSet[i]] + 1;
            if (T < MinCoins[j])
            {
                MinCoins[j] = T;
            }
        }
    }
    
    T = 0;
    for (j=1; j<=Total; j++)
    {
        T += MinCoins[j];
    }
    return T;
}






CoinBruteForce_8c.c

/* Optimal coin set computation by brute force, 8 coin variant */
/* Compile with e.g.:
    gcc -O3 -openmp CoinBruteForce_8c.c -o CoinBruteForce_8c
 * or:
    icc -O3 -openmp CoinBruteForce_8c.c -o CoinBruteForce_8c
    
 * Run as:
    ./CoinBruteForce_8c 
 * or:
    ./CoinBruteForce_8c StartCoin2 EndCoin2 OffsetCoin3
    
    
 * Created by Sebastiaan Breedveld *
   http://sebastiaanbreedveld.nl/
 *
 * License: The Unlicense (http://unlicense.org/): do with it whatever you want.
 *
 */

#include <limits.h>
#include <stdio.h>
#include <omp.h>

/* Maximum change is 499 cents */
#define MaxAmount 499
/* Define number of coins in this system */
#define NrCoins 8

/* uint32 is sufficient and also faster than USHRT */
typedef unsigned int INTTYPE;
#define INTMAX UINT_MAX;

INTTYPE MakeChangeProblem(const INTTYPE *CoinSet); 

int main(int argc, char *argv[])
{ 
    /* Inits */
    INTTYPE CoinSet[NrCoins] = {1, 2, 3, 4, 5, 6, 7, 8};
    INTTYPE TotalCoins, c1, c2, c3, c4, c5, c6, c7;
    INTTYPE cstart=2, cend=MaxAmount, coff2 = 1;
    
    /* We specify a threshold here so it only outputs coin sets which are better than this. You can safely set this to INTMAX */
    INTTYPE MinAmount = 10000;
    
    /* This particular version allows setting an range for the 2nd coin (first coin is always 1) and an offset for the 3rd coin. 
     * It allows running the same executable on multiple machines, each with a different offset */
    if (argc>=2)
    {
      cstart = atoi(argv[1]);
    }
    if (argc>=3)
    {
      cend = atoi(argv[2]);
    }
    if (argc>=4)
    {
      coff2 = atoi(argv[3]);
    }
      
    printf("Going from %u to %u\n", cstart, cend);
    
    for (c1=cstart; c1<=cend; c1++)
    {
        CoinSet[1] = c1;        
        for (c2=c1+coff2; c2<=20; c2++) /* Note that c2 has an upper limit of 20 here - it seems reasonable to assume that this coin will never become larger than 20 cents */
        {
            CoinSet[2] = c2;
            printf("c2 is: %u\n", c2);
            for (c3=c2+1; c3<=50; c3++) /* Ditto */
            {
                for (c4=c3+1; c4<=100; c4++) /* Ditto */
                {
                    /* The next line enables multithreaded execution of the for loop. There is no hard reason
                     * to select this particular for loop for multithreading. If you choose the most inner loop, 
                     * there is an increased overhead of the multithreading, if you choose the top loop, there is
                     * a higher risk of idle threads because some threads are finished sooner than others. */
                    #pragma omp parallel for schedule(dynamic) private(c5, c6, c7, CoinSet, TotalCoins) shared(c1, c2, c3, c4, MinAmount)
                    for (c5=c4+1; c5<=200; c5++)
                    {
                        CoinSet[0] = 1;
                        CoinSet[1] = c1;
                        CoinSet[2] = c2;
                        CoinSet[3] = c3;
                        CoinSet[4] = c4;
                        CoinSet[5] = c5;
                        for (c6=c5+1; c6<=MaxAmount; c6++)
                        {
                            CoinSet[6] = c6;
                            for (c7=c6+1; c7<=MaxAmount; c7++)
                            {
                                CoinSet[7] = c7;
                                TotalCoins = MakeChangeProblem(CoinSet);
                                
                                /* Only print sets which are better than found so far */
                                if (TotalCoins <= MinAmount)
                                {
                                  /* OMP critical: prevent thread collision */
                                  #pragma omp critical
                                  MinAmount = TotalCoins;
                                  printf("1 %3u %3u %3u %3u %3u %3u %3u %6u\n", c1, c2, c3, c4, c5, c6, c7, TotalCoins);
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    return 0;
}

INTTYPE MakeChangeProblem(const INTTYPE *CoinSet)
{
    /* We use a static allocation here, which is faster than dynamic */
    INTTYPE MinCoins[MaxAmount+1];
    INTTYPE i, j, T;
        
    /* This part is taken from the web */
    /* http://codereview.stackexchange.com/questions/92811/find-the-minimum-number-of-coins/92820 */
    MinCoins[0] = 0;
    for (j=1; j<=MaxAmount; j++)
    {
        MinCoins[j] = INTMAX;
    }

    for (i = 0; i < NrCoins; i++)
    {
        for (j = CoinSet[i]; j <= MaxAmount; j++)
        {
            T = MinCoins[j - CoinSet[i]] + 1;
            if (T < MinCoins[j])
            {
                MinCoins[j] = T;
            }
        }
    }
    
    T = 0;
    for (j=1; j<=MaxAmount; j++)
    {
        T += MinCoins[j];
    }
    return T;
}



Thanks to the Online syntax highlighting for the masses! project for code highlighting in HTML.!