Introduction

During the JuliaCon 2023 conference I got several suggestions for writing some more puzzle-solving posts. Therefore this week I want to present a problem that I have recently learned from Paweł Prałat:

Assume that you have an even number n of cords of equal length lying in a bunch. They are arranged in such a way that they are approximately straight so that you see one end of each cord in one region (call it left) and the other end in another region (call it right). Imagine, for example, that the cords were wrapped around in their middles. However, this unfortunately means, you cannot distinguish which end belongs to which cord. Your task is to tie the ends of the cords in such a way that after removing the middle wrapping they form a single big loop. Assuming that you tie the cords randomly (left and right ends separately) compute the probability that you are going to succeed.

The initial setup of the cords (before we start tying them) is shown on a figure below (I took the image from this source):

Cords

As usual, we are going to solve it analytically and computationally.

The post was written using Julia 1.9.2 and DataFrames.jl 1.6.1.

Analytical solution

Denote by p(n) the probability that we succeed with n cords.

Notice that we can assume that we can first tie n right ends of the cords in any order.

Now let us analyze tying the left ends of the cords. We assume that they are tied randomly. We start tying the left ends by randomly picking cord a and tying it to a random cord b. Let us ask when such a tie is a success and when it is a failure.

If n = 2 we know we succeeded. We have just created a loop using two cords. Thus p(2) = 1.

If n > 2 what we want to avoid is a situation when the cords a and b are already tied on the right side. Why? Because then we would create a loop that would be smaller than the required loop including all cords. The probability that we pick a wrong end of a cord is 1/(n-1) as we have n-1 ends to choose from and one of them is bad. Thus we succeed with probability (n-2)/(n-1).

Now observe that, assuming we succeeded, we have just tied three original cords into a one longer cord. Thus we are left with a situation when we have n-2 cords and the problem has the same structure to what we started with. So we get that for n > 2 we can computep(n) = p(n-2)*(n-2)/(n-1).

This means that we can write (using Julia code notation):

p(n) = prod((i-1)/i for i in n-1:-2:3)

Note that this function works correctly only for even n that is at least 4. We will make its more careful implementation in the computational solution section below.

However, first, let us ask ourselves if the formula for this function can be simplified. Indeed we see that it is equivalent to:

prod(i-1 for i in n-1:-2:3) / prod(i for i in n-1:-2:3)

which in turn can be rewritten as:

prod(i-1 for i in n-1:-2:3)^2 / prod(i for i in n-1:-1:2)

Now observe that the numerator is just 2^(n-2)*factorial(n÷2-1)^2 (remember that n is even) and the denominator is factorial(n-1). Thus the formula further simplifies to:

2^(n-2)*factorial(n÷2-1)^2 / factorial(n-1)

And finally:

2^(n-2) / ((n-1)*binomial(n-2, n÷2-1))

Now from Stirling’s approximation we know that:

binomial(n-2, n÷2-1) ~ 2^(2n-4) / sqrt((n-2)*pi)

so for sufficiently large n:

p(n) ~ sqrt((n/2-1)*pi) / (n-1)

Thus we learn that the probability of getting a full circle is of order O(1/sqrt(n)) for large n.

Let us now check these results computationally.

Computational solution

First start with a more careful implementation of the functions computing p(n):

function connect_exact(n::Integer)
    @assert n > 3 && iseven(n)
    return prod((i-1)/i for i in n-1:-2:3)
end

function connect_approx(n::Integer)
    @assert n > 3 && iseven(n)
    return sqrt(pi * (n / 2 - 1)) / (n - 1)
end

Let us check how close the exact and approximate formulas are. Let us compute percentage deviation of the approximation from the exact result:

julia> [1 - connect_approx(n) / connect_exact(n) for n in 4:2:28]
13-element Vector{Float64}:
11377307454724206
060014397013374854
040631211300167
030689300286045995
024649922854770412
02059439568578203
017683822837349372
015493594528168564
013785863139806342
012417071173843053
011295454766000912
01035962441429672
00956696079055186

We see that approximation is slightly below the exact number and that the percentage deviation decreases as n goes up. With n=28 we are below 1% error.

Let us check some larger value of n:

julia> 1 - connect_approx(10000) / connect_exact(10000)
2.5004688326446534e-5

We see that the values are now really close. If you were afraid that we might be hitting numeric computation issues with connect_approx since we are multiplying a lot of values, we can easily switch to a more precise computation with Julia:

julia> 1 - connect_approx(big(10000)) / connect_exact(big(10000))
2.500468833607760982625749174941669517305288399515417883351990349709823151295288e-05

We see that using normal Float64 was enough for this range of values of n to get enough accuracy.

But what if we were not sure if our derivation of the formula for p(n) was correct? We can use simulation to check it.

Here is the implementation of a simulator:

function connect_sim(n::Integer)
    @assert iseven(n) && n > 3
    left = randperm(n)
    neis2 = zeros(Int, n)
    for i in 1:2:n
        neis2[left[i]] = left[i+1]
        neis2[left[i+1]] = left[i]
    end
    prev = 1
    loc = 2
    visited = 2
    while true
        nei1 = isodd(loc) ? loc+1 : loc-1
        nei2 = neis2[loc]
        loc, prev = (prev == nei1 ? nei2 : nei1), loc
        loc == 1 && return visited == n
        visited += 1
    end
end

The the code we assume that we numbered the cords from 1 to n and that in the right part they are connected 1-2, 3-4, … (note that we can always re-number them to get this).

The neis2 keeps the information about connections on left. To get a random connection pattern we first draw a random n-element permutation and store it in the left variable. Then we assume that the connections are formed by cords left[1]-left[2], left[3]-left[4], … and store these connections in the neis2 vector.

Now we are ready to check if this connection pattern is good, that is, it creates one big loop. To do this we start from cord 1 and assume that we first moved to cord 2. The current location of our travel is kept in variable loc. Then from each cord we move either on right or on left to the next cord. The nei1 variable keeps cords neighbor on right and nei2 on left. We keep track in the prev variable which cord we have visited last. Using this information we know which move we should make next. Notice that since we started from 1 we eventually have to reach it. The number of steps taken to reach 1 is tracked by the visited variable. If when loc == 1 we have that visited == n this means that we have formed a big cycle and we return true. Otherwise we return false.

Let us check if our simulation indeed returns values close to theoretical ones. For this we will record the mean of 100,000 runs of our simulation (and here the power of Julia shines - it is not a problem to run that many samples). We check the results for the values of n we investigated above:

using DataFrames
using Random
using Statistics
connect_sim_mean(n) =
    mean(connect_sim(n) for _ in 1:100_000)
Random.seed!(1234)
df = DataFrame(n=[4:2:28; 10_000])
transform(df, :n .=> ByRow.([connect_exact,
                             connect_approx,
                             connect_sim_mean]))

The results of running this code are given below:

14×4 DataFrame
 Row │ n      n_connect_exact  n_connect_approx  n_connect_sim_mean
     │ Int64  Float64          Float64           Float64
─────┼──────────────────────────────────────────────────────────────
│     4        0.666667          0.590818              0.66662
│     6        0.533333          0.501326              0.53622
│     8        0.457143          0.438569              0.45843
│    10        0.406349          0.393879              0.40671
│    12        0.369408          0.360302              0.36978
│    14        0.340992          0.33397               0.33996
│    16        0.31826           0.312631              0.31743
│    18        0.299538          0.294897              0.29848
│    20        0.283773          0.279861              0.28399
│    22        0.27026           0.266904              0.26879
│    24        0.25851           0.25559               0.25528
│    26        0.248169          0.245598              0.24755
│    28        0.238978          0.236692              0.24052
│ 10000        0.0125335         0.0125331             0.01266

We can see that simulation results match the exact calculations well.

Conclusions

I hope you liked the puzzle and the solution. Next week I plan to present the results of some experiments involving machine learning models in Julia.