Introduction

Today im my post I want to concentrate on a less technical topic (at least to start with 😄).

We will discuss solving the following classic puzzle:

You are given a 1 meter long rod. Someone cuts it in n places. The locations of the cuts are picked independently uniformly at random. As a result of the process you get n+1 pieces. What is the probability that it is possible to create a polygon using them?

There are many approaches that can be used to solve this puzzle. I will want to discuss one that uses only basic probability theory and does not require any tricks.

Next we will validate the solution using Julia 1.9.0-rc1. The solution will lead us to another puzzle, that I think is important to understand by people using random number generation in Julia.

Analytical solution

We start with the observation that if we are given n+1 line segments it is possible to create a polygon using them if and only if every of them has less than 50% of total length of the original rod.

Denote the locations of the cut points on the rod by X(i) for i going from 1 to n. We are measuring location from the left of the rod.

When would be the above condition not met? We have two cases (I am ignoring events having zero probability for simplicity of notation):

all X(i) are greater than 0.5 (then the first segment is overly long), or
for some X(i) it is less than 0.5 and there is no other X(j) such that 0 < X(j) - X(i) < 0.5 (then the segment that begins at point X(i) is overly long). We have n such cases.

The second crucial observation is that the n+1 considered events are disjoint. The reason is that it is impossible that two line segments that are longer than 0.5 at the same time.

Therefore the probability that any of the segments is longer than 0.5 is sum of the probabilities of the individual events.

So what is the probability that all X(i) are greater than 0.5? It is relatively easy. We have n independent events, each having probability 0.5 so the total probability is 1/2^n.

And what is the probability that X(i) is less than 0.5 and there is no other X(j) such that 0 < X(j) - X(i) < 0.5? Again we have n independent events (one for location of X(i) and n-1 for location of X(j) for j other than i) and each of these events has 0.5 probability, so again we have that the probability is 1/2^n.

In total we get that the probability that at least one of the segments is longer than 0.5 is (n+1)/2^n, so the probability that we can create a polygon is 1-(n+1)/2^n.

Checking the solution using simulation

Let us write a simple simulation in Julia to check our result. First write a function that takes n points from the [0, 1] interval and checks if they cut the rod into segments from which one could create a polygon:

check(x) = maximum(diff(sort!([0.0; x; 1.0]))) < 0.5

First add start (0.0) and end (1.0) points of the rod. Next, we arrange the passed points in ascending order in-place using sort!. By running the diff operation, we compute the lengths of the created rod segments, and check if longest of them is less than 0.5.

Now we are ready to write a function that runs this simulation multiple times (1 million in my example) to get an estimate of the probability:

runtest(n) = count(_ -> check(rand(n)), 1:10^6) / 10^6

Let us run the simulation for several values of n and compare the results against analytical result:

julia> using Random

julia> Random.seed!(1234);

julia> [(n=n, theory=1-(n+1)/2^n, simulation=runtest(n)) for n in 1:10]
10-element Vector{NamedTuple{(:n, :theory, :simulation), Tuple{Int64, Float64, Float64}}}:
 (n = 1, theory = 0.0, simulation = 0.0)
 (n = 2, theory = 0.25, simulation = 0.250275)
 (n = 3, theory = 0.5, simulation = 0.500066)
 (n = 4, theory = 0.6875, simulation = 0.687842)
 (n = 5, theory = 0.8125, simulation = 0.812825)
 (n = 6, theory = 0.890625, simulation = 0.890247)
 (n = 7, theory = 0.9375, simulation = 0.9374)
 (n = 8, theory = 0.96484375, simulation = 0.964881)
 (n = 9, theory = 0.98046875, simulation = 0.980502)
 (n = 10, theory = 0.9892578125, simulation = 0.989465)

The results look consistent.

Making the simulation faster

Let us benchmark the speed of our code:

julia> @time [(n=n, theory=1-(n+1)/2^n, simulation=runtest(n)) for n in 1:10];
 10.075038 seconds (100.03 M allocations: 4.621 GiB, 10.55% gc time, 0.69% compilation time)

We could make it faster by avoiding unnecessary allocations:

function runtest_faster(n)
    x = zeros(n+2)
    x[end] = 1.0
    v = @view x[2:end-1]
    return count(1:10^6) do _
        rand!(v)
        sort!(v)
        m = 0.0
        xa = first(x)
        @inbounds for i in 2:n+2
            xb = x[i]
            m = max(m, xb - xa)
            xa = xb
        end
        return m < 0.5
    end / 10^6
end

First let us check the results:

julia> Random.seed!(1234);

julia> [(n=n, theory=1-(n+1)/2^n, simulation=runtest_faster(n)) for n in 1:10]
10-element Vector{NamedTuple{(:n, :theory, :simulation), Tuple{Int64, Float64, Float64}}}:
 (n = 1, theory = 0.0, simulation = 0.0)
 (n = 2, theory = 0.25, simulation = 0.250275)
 (n = 3, theory = 0.5, simulation = 0.500066)
 (n = 4, theory = 0.6875, simulation = 0.687842)
 (n = 5, theory = 0.8125, simulation = 0.812825)
 (n = 6, theory = 0.890625, simulation = 0.890247)
 (n = 7, theory = 0.9375, simulation = 0.9374)
 (n = 8, theory = 0.96484375, simulation = 0.964919)
 (n = 9, theory = 0.98046875, simulation = 0.980717)
 (n = 10, theory = 0.9892578125, simulation = 0.989311)

What is interesting is that up to n=7 we get identical results, but they start to differ for n=8. We will investigate it shortly.

However, first we check the timing:

julia> @time [(n=n, theory=1-(n+1)/2^n, simulation=runtest_faster(n)) for n in 1:10];
  1.792463 seconds (22.75 k allocations: 1.533 MiB, 2.80% compilation time)

As we can see there is a significant reduction in the number of allocations.

Investigating random number generation process

Have a look at this code:

julia> Random.seed!(1234); rand(7)
7-element Vector{Float64}:
32597672886359486
5490511363155669
21858665481883066
8942454282009883
35311164439921205
39425536741585077
9531246272848422

julia> Random.seed!(1234); rand!(zeros(7))
7-element Vector{Float64}:
32597672886359486
5490511363155669
21858665481883066
8942454282009883
35311164439921205
39425536741585077
9531246272848422

julia> Random.seed!(1234); rand!(view(zeros(7), :))
7-element view(::Vector{Float64}, :) with eltype Float64:
32597672886359486
5490511363155669
21858665481883066
8942454282009883
35311164439921205
39425536741585077
9531246272848422

All looks good. The results are identical. Now increase the length of the vector by one:

julia> Random.seed!(1234); rand(8)
8-element Vector{Float64}:
5798621201341324
4112941179498505
9721360824554687
014908849285099945
520354993723718
6395615996802734
8396219340580711
967142768915383

julia> Random.seed!(1234); rand!(zeros(8))
8-element Vector{Float64}:
5798621201341324
4112941179498505
9721360824554687
014908849285099945
520354993723718
6395615996802734
8396219340580711
967142768915383

julia> Random.seed!(1234); rand!(view(zeros(8), :))
8-element view(::Vector{Float64}, :) with eltype Float64:
32597672886359486
5490511363155669
21858665481883066
8942454282009883
35311164439921205
39425536741585077
9531246272848422
7955469475347194

What is going on here? We have just hit a puzzle I have promised that is hidden in our original puzzle.

It turns out that the default random number generator in Julia uses a different method for random number generation for Array than for general AbstractArray in some cases.

The relevant code for our case is as follows. The one used by Julia for Vector{Float64}:

function rand!(rng::Union{TaskLocalRNG, Xoshiro},
               dst::Array{Float64},
               ::SamplerTrivial{CloseOpen01{Float64}})
    GC.@preserve dst xoshiro_bulk(rng,
                                  convert(Ptr{UInt8},
                                  pointer(dst)),
                                  length(dst)*8,
                                  Float64,
                                  xoshiroWidth(),
                                  _bits2float)
    dst
end

and the code used for a view:

function rand!(rng::AbstractRNG, A::AbstractArray{T}, sp::Sampler) where T
    for i in eachindex(A)
        x = rand(rng, sp)
        @inbounds A[i] = x
    end
    A
end

Therefore, what you need to remember that if you change in your code rand into rand! to improve performance the results of your computations might change.

Conclusions

I hope you liked the puzzle. If you wonder if the analytical solution could have been simpler the answer is yes. We can get rid of the special cases in our derivation by the following observation. Instead of analyzing cutting a rod in n places you can imagine that you cut a circle in n+1 places. Then we notice that all cutting points are indistinguishable, and what we need to check if there is a gap of 0.5 length from a given point looking clockwise from it. The probability of such gap is 1/2^n as we have n points that cannot lie in a zone of length 0.5. The rest of the arguments in the derivation are the same as above.

Also, as usual, we have a (hopefully useful) lesson about running stochastic simulations in Julia (especially if you want them to be reproducible): rand and rand! are not guaranteed to produce the same sequences of random numbers.