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!