+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: little math/programming problems April 2023 (/thread-19860.html)
little math/programming problems April 2023 - pier4r - 04-28-2023 11:55 AM
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/1651896236457246721?s=20 (contains solutions that spoil the challenge a bit)
RE: little math/programming problems April 2023 - David Hayden - 04-28-2023 10:16 PM
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.08402Dave
RE: little math/programming problems April 2023 - Allen - 04-29-2023 06:50 PM
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.
RE: little math/programming problems April 2023 - Allen - 04-30-2023 01:32 AM
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.05280279510M 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.4931065RE: little math/programming problems April 2023 - David Hayden - 05-02-2023 08:03 PM
(04-29-2023 06:50 PM)Allen Wrote: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?
- 5.490012 expected rolls on average
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
RE: little math/programming problems April 2023 - ThomasF - 05-03-2023 07:10 AM
(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 diceUsing 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 diceEdit: All runs used 5 dice and 10M simulations.
Cheers,
Thomas
RE: little math/programming problems April 2023 - Albert Chan - 05-03-2023 12:18 PM
(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-divisionless-random-integer-generation-on-various-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 ...
RE: little math/programming problems April 2023 - Albert Chan - 05-12-2023 06:14 PM
(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
11111222333334444455666Note 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.
RE: little math/programming problems April 2023 - ttw - 05-14-2023 12:16 PM
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.
RE: little math/programming problems April 2023 - Gilles - 05-15-2023 09:26 PM
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 /
»
»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 ++RE: little math/programming problems April 2023 - pier4r - 05-16-2023 07:47 PM
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.
RE: little math/programming problems April 2023 - Ajaja - 05-16-2023 09:27 PM
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 END100000: 5,49717
1000000: 5,492177
10000000: 5,4922995
RE: little math/programming problems April 2023 - pier4r - 05-17-2023 09:48 PM
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.
RE: little math/programming problems April 2023 - Ajaja - 05-18-2023 07:54 AM
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 ENDSo, standard error of the mean here (SDEV/SQRT(N)):
10000: ~0.05
100000: ~0.016
1000000: ~0.005
10000000: ~0.0016