little math/programming problems April 2023

[pier4r]

I have a problem someone could probably simulate faster than I can work it out. I have a set of dice with values: 0, 0, 0, 0, 2, 3. You roll the dice, add up the total, then roll that number of dice. And repeat. If I start with five dice, how long will the game last?

Source: https://twitter.com/jamesgrime/status/16...46721?s=20 (contains solutions that spoil the challenge a bit)

[David Hayden]

Simulating the experiment 100,000 times, I get an average of 5.49893 rolls with 5 dice. Here are the results with different numbers of starting dice:

Code:


Initial # dice    Average rolls
1                 2.2266
2                 3.25986
3                 4.10384
4                 4.84957
5                 5.49893
6                 6.08402

Someone smarter than me will have to figure out the exact values.

Dave

[Allen]

Maybe possible with the generating function and some ugly recursions?

\( (z^3 + z^2 + 4)^n \)

where \(n \) is the number of dice.
\( f(5) = z^{15} + 5z^{14} + 10z^{13} + 30z^{12} + 85z^{11} + 121z^{10} + 240z^9 + 500z^8 + 480z^7 + 800z^6 + 1280z^5 + 640z^4 + 1280z^3 + 1280z^2 + 1024 \)

(Interesting it's not possible to roll a sum equal to 1 no matter how many dice you use!)

Since the number of dice changes based on the previous roll, you'd need to recursively add up the probabilities for every possible branch in order to calculate an (more?) exact answer - multiplying the conditional probability of each branch times the total chances of arriving at that. Using the initial state of \( f(5) \) and creating a markov chain 24 layers deep you can cover 99% of the expected cases.

Among 10M games with 5 initial dice I got this estimation:

• 5.490012 expected rolls on average
• mode was 2 rolls
• the longest was 81 rolls!!
• 13% of the games ended on the first roll.
• 99% of the games ended in 24 rolls or less.
  

I also tried with 1000 to 1M dice at a time.. the average number of rolls asymptotically approaches something around \( 4.32 * ln(n) \) with 1000-dice games ending in 32 rolls on average.

[Allen]

I believe the first few terms are now calculated (hopefully correct). Starting with 5 dice, the chances of the game terminating on exactly x rolls is:

Code:


x  approx                exact  
0  0.0                    0
1  0.13168724279835392   32 / 3^5
2  0.1762589048590848    614576800 / 3^20
3  0.14345296021217352   1477716325360817158616745850400 / 3^65
4  0.11119048831439719   1127935854140837598232937023804660987401971125340384 / 3^109
...

The total number of states for depth grows quickly to [1, 15, 375, 17655, 1469850, , ...] After that things get out of hand

Here's the difference between the Monte Carlo sampling and the exact calculation:

Code:


halting count for 10,000,000 rounds sampling Monte Carlo starting with 5 dice
rolls    count    Monte Carlo    Empirical    Abs Error    Percent error
1    1317393    0.1317393      0.131687243    -5.20572E-05    -0.039530938
2    1761181    0.1761181      0.176258905     0.000140805     0.079885246
3    1434023    0.1434023      0.14345296      5.06602E-05     0.03531486
4    1112492    0.1112492      0.111190488    -5.87117E-05    -0.052802795

10M rolls Total Expected Value 5.4931065

Code:


rolls    count    Monte Carlo    PDF    cumulative
1    1317393    0.1317393    0.1317393    0.1317393
2    1761181    0.1761181    0.3522362    0.4839755
3    1434023    0.1434023    0.4302069    0.9141824
4    1112492    0.1112492    0.4449968    1.3591792
5    862577    0.0862577    0.4312885    1.7904677
6    673073    0.0673073    0.4038438    2.1943115
7    531444    0.0531444    0.3720108    2.5663223
8    424292    0.0424292    0.3394336    2.9057559
9    339510    0.033951    0.305559    3.2113149
10    274634    0.0274634    0.274634    3.4859489
11    223553    0.0223553    0.2459083    3.7318572
12    183459    0.0183459    0.2201508    3.952008
13    149193    0.0149193    0.1939509    4.1459589
14    123336    0.0123336    0.1726704    4.3186293
15    100750    0.010075    0.151125    4.4697543
16    83419    0.0083419    0.1334704    4.6032247
17    68726    0.0068726    0.1168342    4.7200589
18    57454    0.0057454    0.1034172    4.8234761
19    47370    0.004737    0.090003    4.9134791
20    39204    0.0039204    0.078408    4.9918871
21    32289    0.0032289    0.0678069    5.059694
22    26840    0.002684    0.059048    5.118742
23    22160    0.002216    0.050968    5.16971
24    18664    0.0018664    0.0447936    5.2145036
25    15493    0.0015493    0.0387325    5.2532361
26    12861    0.0012861    0.0334386    5.2866747
27    10782    0.0010782    0.0291114    5.3157861
28    9111    0.0009111    0.0255108    5.3412969
29    7672    0.0007672    0.0222488    5.3635457
30    6240    0.000624    0.01872    5.3822657
31    5192    0.0005192    0.0160952    5.3983609
32    4240    0.000424    0.013568    5.4119289
33    3544    0.0003544    0.0116952    5.4236241
34    2980    0.000298    0.010132    5.4337561
35    2461    0.0002461    0.0086135    5.4423696
36    2093    0.0002093    0.0075348    5.4499044
37    1686    0.0001686    0.0062382    5.4561426
38    1459    0.0001459    0.0055442    5.4616868
39    1204    0.0001204    0.0046956    5.4663824
40    1034    0.0001034    0.004136    5.4705184
41    809    0.0000809    0.0033169    5.4738353
42    715    0.0000715    0.003003    5.4768383
43    565    0.0000565    0.0024295    5.4792678
44    448    0.0000448    0.0019712    5.481239
45    396    0.0000396    0.001782    5.483021
46    326    0.0000326    0.0014996    5.4845206
47    311    0.0000311    0.0014617    5.4859823
48    235    0.0000235    0.001128    5.4871103
49    172    0.0000172    0.0008428    5.4879531
50    142    0.0000142    0.00071    5.4886631
51    135    0.0000135    0.0006885    5.4893516
52    102    0.0000102    0.0005304    5.489882
53    73    0.0000073    0.0003869    5.4902689
54    78    0.0000078    0.0004212    5.4906901
55    74    0.0000074    0.000407    5.4910971
56    58    0.0000058    0.0003248    5.4914219
57    45    0.0000045    0.0002565    5.4916784
58    40    0.000004    0.000232    5.4919104
59    34    0.0000034    0.0002006    5.492111
60    22    0.0000022    0.000132    5.492243
61    22    0.0000022    0.0001342    5.4923772
62    20    0.000002    0.000124    5.4925012
63    21    0.0000021    0.0001323    5.4926335
64    16    0.0000016    0.0001024    5.4927359
65    5    0.0000005    0.0000325    5.4927684
66    9    0.0000009    0.0000594    5.4928278
67    5    0.0000005    0.0000335    5.4928613
68    6    0.0000006    0.0000408    5.4929021
69    4    0.0000004    0.0000276    5.4929297
70    3    0.0000003    0.000021    5.4929507
71    5    0.0000005    0.0000355    5.4929862
72    2    0.0000002    0.0000144    5.4930006
73    2    0.0000002    0.0000146    5.4930152
74    5    0.0000005    0.000037    5.4930522
76    2    0.0000002    0.0000152    5.4930674
77    2    0.0000002    0.0000154    5.4930828
78    2    0.0000002    0.0000156    5.4930984
81    1    0.0000001    0.0000081    5.4931065

[David Hayden]

(04-29-2023 06:50 PM)Allen Wrote:  

• 5.490012 expected rolls on average
  

I'm surprised that we came up with such different answers. With 10 million trials, I got an average of 5.49603. Here's my C++ code. Does anyone spot a problem? Maybe this is a case of rand() showing it's limitations?

Code:

#include<iostream>                               //Required for cin, cout
#include<cmath>
using namespace std;
int
main()
{
    int sides[6] = {0,0,0,0,2,3};
    int N, D;
    int totalRolls = 0;
    cout << "Enter starting number of dice and number of simulations: ";
    cin >> D >> N;

    for (int i=0; i<N; ++i) {
        int sum = D, count=0;
        int dice;
        for (count=0; (dice = sum); ++count) {  // note that (dice=sum) is assignment, not comparison
            for (int i=sum=0; i<dice; ++i) {
                int side = sides[rand() % 6];
                sum += side;
            }
        }
        totalRolls += count;
    }
    cout << "Average of " << (double)totalRolls/N << " rolls starting with " << D << " dice\n";
}

$ ./foo
Enter starting number of dice and number of simulations: 5 10000000
Average of 5.49603 rolls starting with 5 dice

[ThomasF]

(05-02-2023 08:03 PM)David Hayden Wrote:  Here's my C++ code. Does anyone spot a problem? Maybe this is a case of rand() showing it's limitations?

Hi David!

I assume this is due to the limitations of rand().
With different seeds I get the following result (ie your code but added seed initialization).

Code:

  srand (time(NULL));

Code:

Average of 5.48456 rolls starting with 5 dice
Average of 5.48837 rolls starting with 5 dice
Average of 5.48826 rolls starting with 5 dice

Using another random generator (PCG), I get the following result:

Code:

Average of 5.49419 rolls starting with 5 dice
Average of 5.49255 rolls starting with 5 dice
Average of 5.49275 rolls starting with 5 dice

Edit: All runs used 5 dice and 10M simulations.

Cheers,
Thomas

[Albert Chan]

(05-02-2023 08:03 PM)David Hayden Wrote:  Does anyone spot a problem? Maybe this is a case of rand() showing it's limitations?

rand(), and many other random number generators, have weaker random low bits.
We should use random high bits instead.

< int side = sides[rand() % 6];
> int side = sides[rand() * 6ull / (RANDMAX+1u)];

Bonus: (RAND_MAX+1u) is likely pow-of-2, division should compiled as fast bits shift.
Note: we can speed up a bit more, replacing 6ull by 6u, if it does not cause overflow.

I have ignored slight bias mapping rand() to 0 .. 5
see http://www.azillionmonkeys.com/qed/random.html
see https://lemire.me/blog/2019/06/06/nearly...s-systems/

However, because of long tails, it is hard to get good estimate, even with good random numbers.
For 5 dices, correct expected value = 5.491 819 ...

[Albert Chan]

(05-02-2023 08:03 PM)David Hayden Wrote:  int sides[6] = {0,0,0,0,2,3};

For a rough idea of expected values, we may try a simple "random" generator.
5 dices, with cycled pattern, different initial state:

Code:

000023000023000023000023000023000023 ... Pattern
1111122333
 111112222233444
  11111222223333344555
   1111122222333334444455666
    111112222233333444445555566777
     11111222333334444455666

Note that 6 throws add to sum of only 5, sequence will always terminate.
Expected values = (3 + 4 + 5 + 6 + 7 + 6) / 6 = 31/6 ≈ 5.2

This "random" generator is super lumpy, with maximal 4 0's in a row.
We expected true expected values slightly bigger.

[ttw]

The following pRNG may be helpful if sampling the dice throws is a problem. It's a specially constructed linear congruential generator designed for die rolling.

Letting X stand for the internal state, one has

X = X * 93847949523997033 Mod ( 216172782113784439 )

The die roll is X ( Mod 6 ) or X ( Mod 6 ) + 1

The only constraint is that the first X must be non-zero 0 < Xo < 93847949523997033.

If this version has distributional problems, I got millions more available.

[Gilles]

My (New)RPL version :

Code:

« → Loop 
  « 0 
    1 Loop START
      5 'NbDe' LSTO  
      0 'TotalJet' LSTO 
      WHILE NbDe REPEAT
        1 'TotalJet' STO+  0 'NewDe' LSTO 
        1 NbDe START
          RAND 6 * 
          CASE
            DUP 5 ≥ THEN 3 'NewDe' STO+ END
            DUP 4 ≥ THEN 2 'NewDe' STO+ END
          END
          DROP 
        NEXT
        NewDe 'NbDe' STO 
      END
      TotalJet + 
    NEXT
    Loop / 
  »
»

Size : 296 Bytes

Execution on my laptop and HP50g, NewRPL

Code:

Loops    Time(s)    Result  Time 
         Laptop             HP50g
----------------------------------
1000     0.015      5.438    4.1
10000    0.285      5.4475   39
100000   2.912      5.48428  389 (ratio vs PC ~ 133 )
1E6      30.99      5.49546  ++
10E6     311.79     5.49432  ++

[pier4r]

Finally a calculator version as well, thanks Gilles! (and all the others as well of course)

Admittedly I could produce a calculator version too, but I didn't prioritize the task.

[Ajaja]

My HP-42s/Free42 version:

Code:

00 { 61-Byte Prgm }
01▸LBL "DICE"
02 CLRG
03 STO 11
04 STO 09
05 2
06 STO 04
07 3
08 STO 05
09▸LBL 04
10 0
11 STO 06
12 5
13 STO 07
14▸LBL 03
15 0
16 STO 08
17▸LBL 01
18 RAN
19 6
20 ×
21 RCL IND ST X
22 STO+ 08
23 DSE 07
24 GTO 01
25 1
26 STO+ 06
27 RCL 08
28 STO 07
29 X>0?
30 GTO 03
31 RCL 06
32 STO+ 10
33 DSE 09
34 GTO 04
35 RCL 10
36 RCL 11
37 ÷
38 END

10000: 5,4555
100000: 5,49717
1000000: 5,492177
10000000: 5,4922995

[pier4r]

Interesting the difference between the newRPL and the Free42.

I thought that the math library (at least for basic stuff like RAND) of the HP42 and RPL was similar but likely it is not.

[Ajaja]

I think it just has a large deviation.
I've changed the program a bit to calculate SDEV:

Code:

00 { 57-Byte Prgm }
01▸LBL "DICE"
02 CLRG
03 STO 09
04 2
05 STO 04
06 3
07 STO 05
08▸LBL 04
09 CLX
10 STO 06
11 5
12 STO 07
13▸LBL 03
14 CLX
15 STO 08
16▸LBL 01
17 RAN
18 6
19 ×
20 RCL IND ST X
21 STO+ 08
22 DSE 07
23 GTO 01
24 1
25 STO+ 06
26 RCL 08
27 STO 07
28 X>0?
29 GTO 03
30 RCL 06
31 PRX
32 Σ+
33 DSE 09
34 GTO 04
35 SDEV
36 END

It's ~5
So, standard error of the mean here (SDEV/SQRT(N)):
10000: ~0.05
100000: ~0.016
1000000: ~0.005
10000000: ~0.0016