Re-posted from: https://bkamins.github.io/julialang/2023/10/20/pe.html
Introduction
After several technical posts today I decided to switch back to puzzle-solving mode.
My readers probably know that I like and promote the Project Euler project.
This time I picked a relatively easy puzzle 205 as it nicely shows several
functionalities of Julia that are worth learning.
The post was written under Julia 1.9.2 and StatsBase.jl v0.34.2.
The problem
The problem 205 is stated as follows:
Peter has nine four-sided dice, each with faces numbered from 1 to 4.
Colin has six six-sided dice, each with faces numbered from 1 to 6.
Peter and Colin roll their dice and compare totals: the highest total wins.
The result is a draw if the totals are equal.
What is the probability that Peter beats Colin?
Simulation approach
To get some intuition for our problem let us first try simulating the distribution of the results.
We will draw one million times from Peter’s and Colin’s dice:
julia> using Random
julia> using StatsBase
julia> Random.seed!(1234);
julia> sim_p = [sum(rand(1:4) for _ in 1:9) for _ in 1:10^6]
1000000-element Vector{Int64}:
24
23
20
26
23
⋮
20
27
22
24
23
julia> sim_c = [sum(rand(1:6) for _ in 1:6) for _ in 1:10^6]
1000000-element Vector{Int64}:
21
23
23
24
24
⋮
16
30
15
20
25
Now let us compare the results:
julia> describe(sim_p)
Summary Stats:
Length: 1000000
Missing Count: 0
Mean: 22.498804
Std. Deviation: 3.355036
Minimum: 9.000000
1st Quartile: 20.000000
Median: 22.000000
3rd Quartile: 25.000000
Maximum: 36.000000
Type: Int64
julia> describe(sim_c)
Summary Stats:
Length: 1000000
Missing Count: 0
Mean: 20.996381
Std. Deviation: 4.186906
Minimum: 6.000000
1st Quartile: 18.000000
Median: 21.000000
3rd Quartile: 24.000000
Maximum: 36.000000
Type: Int64
julia> mean(sim_p .> sim_c)
0.5728
We see that Peter’s chances of winning are around 57%.
We also see that both the mean and the median of Peter’s dice are better than Colin’s.
However, the simulation results are only approximate. Let us thus compute the exact result.
Exact approach
To compute the exact probability of Peter’s win first calculate the distribution of
Peter’s and Colin’s dice. With StatsBase.jl it is easy to do using the countmap function:
julia> ex_p = countmap(map(sum, Iterators.product((1:4 for _ in 1:9)...)))
Dict{Int64, Int64} with 28 entries:
16 => 4950
20 => 23607
35 => 9
12 => 165
24 => 27876
28 => 8451
30 => 2598
17 => 8451
23 => 30276
19 => 18351
22 => 30276
32 => 486
11 => 45
36 => 1
9 => 1
31 => 1206
⋮ => ⋮
julia> ex_c = countmap(map(sum, Iterators.product((1:6 for _ in 1:6)...)))
Dict{Int64, Int64} with 31 entries:
16 => 2247
20 => 4221
35 => 6
12 => 456
24 => 3431
28 => 1161
8 => 21
17 => 2856
30 => 456
23 => 3906
19 => 3906
22 => 4221
32 => 126
6 => 1
11 => 252
36 => 1
⋮ => ⋮
As a result we get the number of times out of the total possible outcomes that a given sum on a dice occurs.
Let us check that the total counts of the outcomes are correct. We can do it as we know that on Peter’s dice
we can get 4^9 outcomes and on Colin’s dice 6^6:
julia> sum(values(ex_p)), 4^9
(262144, 262144)
julia> sum(values(ex_c)), 6^6
(46656, 46656)
Indeed the results match.
Using distributions stored in ex_p and ex_c variables we can count the number of times Peter wins with Collin
using the Cartesian product of the distributions (the tosses of the dice are independent) and, in consequence compute the
exact probability that Peter wins:
julia> p_win = sum(pk > ck ? pv*cv : 0 for
(pk, pv) in pairs(ex_p),
(ck, cv) in pairs(ex_c))
7009890480
julia> total = 4^9 * 6^6
12230590464
julia> p_win / total
0.5731440767829801
julia> p_win // total
48679795//84934656
Note that in Julia we can nicely compute both approximate solution (using Float64) and exact solution using rational numbers.
One important aspect that we should have checked when solving this puzzle is if we do not have an integer overflow issue when
computing the result. On 64-bit machine the overflow happens for the value:
julia> typemax(Int)
9223372036854775807
which is much larger than 12230590464. But how can we be sure that we do not get an overflow when computing 4^9 * 6^6?
The easiest check is to take the logarithms of both expressions:
julia> log(typemax(Int))
43.66827237527655
julia> 9 * log(4) + 6 * log(6)
23.227206065447348
Indeed we see that we have a wide safety margin.
Note, however, that on 32-bit machine Julia would use Int32 type to represent integers by default, and hen we have:
julia> typemax(Int32)
2147483647
julia> log(typemax(Int32))
21.487562596892644
And in this case we would have an integer overflow issue, so some care is needed.
Conclusions
I hope you enjoyed the puzzle, the solution, and the code examples I have presented!